Breast cancer is the second leading cause of cancer death among all cancers occurring in American women (1). The survival rate and prognosis of a patient with breast cancer is dependent on the stage of cancer at diagnosis. Locally advanced breast cancers (generally with tumor size >2 cm) are often treated with neoadjuvant chemotherapy (NACT) before surgery to reduce the tumor size for breast-conserving surgery (2, 3). A pathological complete response (pCR) to NACT is considered a surrogate marker for overall and long-term disease-free survival (4). However, the pCR rate is only 6%–45% depending on breast cancer subtypes and treatment regimen (5, 6, 7, 8). It is therefore important to identify the nonresponders at an early stage so that their treatment regimen can be modified, sparing them the long- and short-term toxicities from ineffective chemotherapies. Currently, in standard of care, the response to NACT is evaluated based on the histological examination of a surgical specimen taken after the completion of NACT. Noninvasive or minimally invasive methods that can predict therapy response at the early stages of NACT can potentially play an important role in the emerging era of precision medicine to help guide regimen de-escalation/alteration in NACT treatment of breast cancer (9).
A significant change in the microenvironment of the tumor, such as perfusion/permeability and metabolism, usually precedes a reduction in tumor size as response to chemotherapy (10–13). As a noninvasive imaging method for assessment of microvascular perfusion/permeability, dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) is increasingly used in research and early-phase clinical trial settings to predict and evaluate cancer response to treatment (9, 10). Several studies (9, 14–22) have shown that changes in quantitative parameters estimated from pharmacokinetic modeling of DCR-MRI data can be useful markers for early prediction of breast cancer response to NACT.
When compared to normal tissue vasculature, tumor vasculature exhibits greater spatiotemporal heterogeneity. The heterogeneity of the tumor vasculature reflects the tumor stage and the disease progression (23). The aforementioned DCE-MRI studies (9, 14–22) generally reported changes in mean parameter values of the entire breast tumor, masking the potential changes in spatial heterogeneity of the microvasculature in response to NACT. Image texture features that can capture the heterogeneity of tumor vasculature from DCE-MRI images or voxel-based parametric maps could be highly useful in assessing tumor response to therapy. Several texture analysis methods such as gray-level co-occurrence matrix (GLCM) and gray-level run length matrix (RLM) have been frequently used in DCE-MRI analysis (24, 25). They were initially used on DCE-MRI images directly. Teruel et al. (24) analyzed T1-weighted DCE-MRI images using GLCM features to predict breast cancer response to NACT. They extracted 16 textural features at each time point of a DCE-MRI acquisition, and the most significant feature yielded a receiver operating curve (ROC) area under the curve (AUC) of 0.77 for prediction of pCR versus stable disease. Similarly, Golden et al. (25) used GLCM features from pre- and post-NACT 2-dimensional (2D) DCE image slices to evaluate NACT response. The pre-NACT features were able to predict pCR with an AUC of 0.68. Although the post-NACT features showed more favorable performances in predicting pCR, these were obtained after the completion of NACT and were not useful for early prediction of pCR. Several studies have performed the same texture analysis on voxel-based maps of pharmacokinetic parameters estimated from pharmacokinetic modeling of DCE-MRI data. Banerjee et al. (26) extracted a combination of intensity, texture, shape, and edge-based features from 2D maps of pharmacokinetic parameters before and after NACT to assess treatment response. Their best model obtained an AUC of 0.83, using a concatenation of Riesz and first-order statistical features. However, the use of only pre- and post-NACT data limits the utility of this model for early prediction of NACT response. In our previous study, we (27) have extracted multiple statistical texture features from 3-dimensional (3D) pharmacokinetic parametric maps before and after 1 cycle of NACT, and found that 3D GLCM features were most effective for early prediction of NACT response through correlation with index values of residual cancer burden (RCB) using a regression model.
In all the analysis methods described above, texture has been studied on a statistical level, by analyzing the spatial distribution of the gray-level values. Textures can also be characterized by fractals, which describe irregular structures that show self-similarity at various scales. Fractal-based texture analysis correlates texture heterogeneity to fractal dimension (FD), which is a mathematical descriptor of a structure's geometrical complexity, based on the concept of spatial pattern self-similarity. Rose et al. (28) showed that fractal analysis could be used to quantify spatial heterogeneity in DCE-MRI parametric maps and differentiate between low- and high-grade tumors. Several other studies (29, 30) have used fractal analysis of breast DCE-MRI images to classify benign versus malignant tumors.
Another important aspect while considering textures is the scale. It has been shown that human visual system processes information in a multiscale approach (different cells in the visual cortex respond to different frequencies and orientations) (31). Owing to the highly heterogeneous nature of the tumor vasculature, analyzing images at a single resolution may not be able to capture the entire complexity of the tumor vasculature. A multiresolution approach can decompose an image into different levels of resolution, giving an opportunity to extract informative features at each level. Lower resolution levels best represent large structures or high contrast, while higher resolutions describe small size or low-contrast objects (32). Multiresolution analysis gives the advantage of analyzing both small- and large-object characteristics in a single image at several resolutions and therefore may be better suited to describe the highly heterogeneous tumor vasculature structure. Multiresolution methods, such as the wavelet analysis, transform images into a representation containing both frequency and spatial information (33). The mean and entropy values extracted from the subimages resulting from wavelet decomposition of DCE-MRI images have been used to classify malignant and benign breast tumors (34, 35). Braman et al. (36) used Gabor wavelet, co-occurrence measures and energy measures to generate 1980 features from DCE images to predict breast cancer response to NACT. A feature selection step was carried out to select top 10 features for final classification. Al-Kadi et al. combined wavelet analysis with fractal analysis and used multiresolution fractal descriptors on ultrasonography images to characterize the tissue and showed that tumor heterogeneity described by this feature improved prediction of response to therapy and disease characterization (37). To the best of our knowledge, fractal analysis at multiple resolutions has not been conducted on breast MRI images for prediction of response to NACT. In this preliminary study, we evaluated the potential of multiresolution fractal analysis of volumetric DCE-MRI pharmacokinetic parametric maps for early prediction of breast cancer response to NACT, and compared it with the conventional methods of GLCM, RLM and single-resolution fractal analysis.
Materials and Methods
Patient Cohort and Study Schema
In total, 55 patients diagnosed with locally advanced breast cancer received standard-of-care NACT. They were consented to participate in a longitudinal research DCE-MRI study approved by the local IRB. The NACT regimen typically consists of 4 cycles of doxorubicin–cyclophosphamide administration every 2 weeks followed by 4 cycles of taxane every 2 weeks, or 6 cycles of the combination of all 3 drugs every 3 weeks (9, 27). The targeted agent trastuzumab was added to the regimen for tumors with positive HER2 (human epidermal growth factor receptor 2) receptor status. A full NACT course therefore would normally last 4–5 months.
In total, 4 DCE-MRI examinations were performed before, during, and after the NACT course: pre-NACT (visit-1), after the first NACT cycle (visit-2), at NACT midpoint (visit-3; usually after 3 or 4 cycles of NACT, or before the change of NACT agents), and after the completion of NACT but before surgery (visit-4). Except for the visit-1 examination, all examinations were performed at least a week after administering the latest cycle of NACT agents to allow time for the drugs to take effect. Pathological analysis of the post-NACT surgical specimens was performed to determine the status of pathologic response to NACT. The values of cross-sectional size of the tumor in 2D, tumor cell density, number of lymph nodes involved, and the greatest dimension in the largest involved node were measured and used in the equation given by Symmans et al. (38) to compute the RCB. A pCR is defined as the absence of residual invasive tumor, indicated by RCB = 0. Non-pCR includes all cases with RCB > 0.
In this preliminary study, only data from the visit-1 and visit-2 DCE-MRI studies were used for image feature analysis and correlations with response endpoint of pCR versus non-pCR to assess the capability for early prediction of breast cancer response to NACT.
DCE-MRI Data Acquisition and Analysis
DCE-MRI data acquisition was performed using a Siemens 3 T system (Siemens, Erlangen, Germany) with the body coil as the transmitter and a 4-channel bilateral phased-array breast coil as the receiver. During each MRI session, following pilot scans and precontrast axial T1- and T2-weighted MRI acquisitions, axial bilateral DCE-MRI images with full breast coverage were acquired using a 3D gradient echo-based Time-resolved angiography With Stochastic Trajectories (TWIST) sequence (9). DCE-MRI acquisition parameters included the following: flip angle = 10°, echo time/repetition time = 2.9/6.2 milliseconds, parallel imaging acceleration factor of 2, field of view = 30 to 34 cm, in-plane matrix size = 320 × 320, and slice thickness = 1.4 mm. About 32–34 image volume sets of 104–128 slices each were acquired over a period of about 10 minutes with a temporal resolution of 14–20 seconds. The contrast agent gadolinium (HP-DO3A) was injected intravenously (0.1 mmol/kg at 2 mL/s) using a programmable power injector after acquisition of 2 baseline image volumes, followed by a 20-mL saline flush at the same injection rate.
Three experienced breast radiologists manually delineated the tumor region of interest (ROI) on postcontrast (90–120 seconds after the injection of the contrast agent) DCE-MRI image slices that contained the contrast-enhanced tumor. To minimize interobserver variability in tumor ROI drawing for the same patient, 1 radiologist drew ROIs for the entire longitudinal study of a single patient. With only 2 patients having multifocal disease, ROIs were drawn for the primary breast tumors only. Figure 1 shows an example of postcontrast DCE images from a pCR patient, drawn tumor ROIs, and ROI mean signal intensity ratio time-courses at visit-1 and visit-2. For pharmacokinetic analysis, precontrast tissue T1 value, T10, was determined using a proton density method (9) by acquiring proton density images just before DCE-MRI that were spatially coregistered with the DCE images. The DCE time-course data from the voxels within the tumor ROI was fitted with a 2-compartment—3-parameter shutter-speed model (9, 39), using a population-averaged arterial input function from the axillary artery (9). This pharmacokinetic analysis yielded the following 4 parameters: Ktrans (volume transfer rate constant), ve (volume fraction of extravascular and extracellular space), kep (=Ktrans/ve, efflux rate constant), and τi (mean intracellular water lifetime). Figure 2 shows examples of voxel-based parametric maps of these 4 parameters for a pCR (Figure 2A) and a non-pCR (Figure 2B) tumor at visit-1 and visit-2. The parametric maps from the visit-1 and visit-2 studies of all patients were subjected to multiresolution fractal analysis described in detail below, as well as the traditional texture analysis methods of GLCM, RLM and single-resolution fractal analysis.
Multiresolution Fractal Analysis
Each of the parametric maps was decomposed into a multiresolution representation using wavelet analysis, and subsequently, FDs were calculated at each resolution level.
Wavelet Analysis for Multiresolution Decomposition.
Wavelet analysis is used to decompose the parametric maps into a set of frequency sub-bands based on small basis functions of varying frequency and limited time duration called wavelets, enabling the characterization of texture at appropriate frequency levels. The wavelet is scaled and translated to cover the time-frequency domain. The discrete wavelet transform for a function f(x, y, z) of size (M, N, K) can be represented as:40). The wavelet transform depends mainly on the scaling (φ) and wavelet (ψ) functions, but it is not necessary to define their explicit form. Instead, a low-pass and high-pass filter that characterize the interaction of these functions are used. The process of decomposing the parametric maps can be viewed as passing them through a series of low-pass and high-pass filters and down-sampling successively. The 3D volume is first filtered along the columns resulting in a low-pass-filtered subvolume and a high-pass-filtered subvolume. These resulting subvolumes are further filtered along rows and slices resulting in 8 decomposed subvolumes. We have used Daubechies wavelets, as these filters have been designed to account for signal discontinuities and self-similarity, which make them the most suitable wavelet for describing signals exhibiting fractal patterns (41). Unlike Haar wavelet, they use overlapping windows that help capture changes in high frequency, and they also demonstrate better recognition of fine characteristic structures (40). One level of decomposition results in 8 subvolumes. The FD for each of these subvolumes is calculated.
Multiresolution Fractal Analysis.
The FD is calculated based on the power spectrum analysis of the 3D Fourier transformation of the subvolumes (42). The 3D discrete Fourier transform is defined as:
The frequency space is evenly divided into 12 zenith and 24 azimuth directions, and 30 points are uniformly sampled along the radial component in each of these directions. The power spectral density is plotted against sampled radial frequency in a log-log plot. The slope β of a least-squares regression line of the log-log plot is related to the FD as:
In standard wavelet analysis, the energy of the subvolume is used to guide further decomposition, but this value is highly dependent on the intensity values of the subvolume. In this work, instead of using energy, the subvolume with the highest FD was selected for further decomposition.
Each parametric map was decomposed down to 4 levels, using FD to guide the sub-band tree structure. Finally, we concatenated the highest and the lowest FD at each level of decomposition to form a feature vector. Therefore, for each parametric map an 8-dimensional feature vector is generated from multiresolution FD analysis.
Conventional Texture Feature Analysis
We compared the performance of multiresolution fractal analysis features with that of GLCM, RLM, and single-resolution fractal analysis. GLCM is a second-order statistical method, which estimates the joint probability P(i, j | d θ), where 2 voxels with intensity i and j are separated by distance d and direction θ. A GLCM matrix was constructed by averaging the matrices obtained over 13 directional offsets at distance d = 1 (27). Twelve Harlick features (43) were derived from this GLCM matrix. RLM P(i, r | θ) is defined as the number of pixels with gray-level i and run-length r, for a given direction θ. RLM was computed by adding all possible run lengths in the 13 directions of the 3D space and 13 statistical features were derived from this matrix (44). Fractal analysis describes the roughness or smoothness of the texture through the FD measure. Here, single-resolution fractal analysis refers to the estimation of FD of the tumor ROI from 3D parametric maps directly (39).
Evaluation of Predictive Performance for NACT Response
For each of the features obtained from the GLCM, RLM, multiresolution, and single-resolution fractal analysis, the percentage change in the feature values was calculated between the visit-1 and visit-2 DCE-MRI studies. These percentage changes were given as input to support vector machine (45), a robust classifier, to generate a predictive model for classification of pCR versus non-pCR. The performances of the models were evaluated using the ROC AUC, sensitivity and specificity analysis. Sensitivity here refers to the proportion of pCRs correctly identified as pCRs, while specificity refers to the proportion of non-pCRs correctly identified as non-pCRs.
The support vector machine classification performance was evaluated by calculating the average over 10 random partitions of the data for training and testing. For each partition, pCRs and non-pCRs were randomly divided into training and testing data sets as described below. The mean and standard deviation of AUC, sensitivity, and specificity values obtained over the 10 partitions of training and testing data sets are reported. The predictive performance was assessed for the features extracted from each of the 4 parametric maps as well as those constructed by concatenating the texture features from all 4 parametric maps of Ktrans, kep, ve, and τi, designated as “All.”
The ROC AUC values of the multiresolution fractal features were compared with those of the conventional features by calculating the critical ratio according to the Hanley and McNeil formula (46). The statistical significance was set at P < .05.
Among the 55 patients in the study cohort, 14 achieved pCR to NACT, while the other 41 patients were non-pCRs based on pathological analysis of the surgical specimens. Table 1 shows the clinicopathological characteristics of the pCR and non-pCR groups. Nine pCRs/31 non-pCRs and 5 pCRs/10 non-pCRs were selected randomly to form the training and testing sets, respectively. Figure 3 shows the ROC AUC values for classification of pCR versus non-pCR using the GLCM, RLM, single-resolution fractal, and multiresolution fractal features from parametric maps of different DCE-MRI parameters considered individually and the concatenated feature from all 4 parametric maps. GLCM and RLM features seemed to overfit on the training data, as they had high training AUCs but low testing AUCs. For example, the AUC values were 0.76 and 0.75 from the training Ktrans maps, and 0.47 and 0.40 from the testing Ktrans maps for the GLCM and RLM methods, respectively. Overall, the multiresolution fractal features from each DCE-MRI parametric map and the concatenated features performed the best in prediction of pCR versus non-pCR in both the training and testing data sets with AUC = 0.85, 0.86, 0.87, 0.86, 0.91 (for Ktrans, kep, ve, τi, and All, respectively), and 0.80, 0.63, 0.74, 0.70, 0.78 for the training and testing data sets, respectively. The only exception was from the kep maps in the testing data sets, where the single-resolution fractal analysis provided the highest AUC of 0.71 among the 4 feature analysis methods. Within the testing or training sets, the predictive performances of multiresolution fractal features were significantly better than the GLCM features from the Ktrans map (AUC = 0.47, P = .022) in the testing set, RLM features from the τi map (AUC = 0.71, P = .012) in the training set, RLM features from the Ktrans (AUC = 0.40, P = .013), and “All” (AUC = 0.47, P = .049) maps in the testing set, and single-resolution fractal features from the ve (AUC = 0.71, P = .012) and τi (AUC = 0.67, P = .037) maps in the training set.
We evaluated the specificities of the classification models at 2 levels of sensitivities for the testing data sets: 60% (3 out of 5 pCRs were classified correctly) and at 80% (4 out of 5 pCRs were classified correctly), as shown in Table 2. At both sensitivity levels, with a few exceptions, fractal features presented higher specificities than the GLCM and RLM features, with the multiresolution method generally outperforming the single-resolution method (except when they were applied to the kep map).
Figure 4 shows the ROC AUC values (for the training and testing data sets) for each level of decomposition of the concatenated feature vectors obtained by combining multiresolution fractal features extracted from the Ktrans, kep, ve, and τi maps. The combination of features from all 4 levels and all 4 parametric maps is represented by the black bar “All.” It can be observed that the first level of decomposition alone achieved a predictive performance (AUC = 0.87 and 0.82 for the training and testing sets, respectively) comparable to that of the combined features from all levels (AUC = 0.91 and 0.78 for the training and testing sets, respectively), while features from levels 2, 3, and 4 individually had rather poor performances in the test set with AUC = 0.57, 0.45, and 0.54, respectively.
This preliminary study shows that multiresolution fractal analysis has the potential to better capture the heterogeneity in the breast tumor vasculature as measured by DCE-MRI, and the extracted features from voxel-based DCE-MRI parametric maps are good early predictors of breast cancer response to NACT. In general, the concatenated features extracted from parametric maps of all the DCE-MRI parameters provide the best predictive performance. Multiresolution analysis filters out irrelevant features and noise at different resolutions, rendering more emphasis on distinct features, and fractal analysis at each level appears to be able to capture these distinct features. The GLCM and RLM features reflect the overall correlation between adjacent voxels in terms of second-order and higher-order statistical features, respectively (27). For the small data set used in this study, the generally higher AUC values from the multiresolution fractal analysis when compared to GLCM and RLM methods suggest that decomposing the texture may give further insights into the heterogeneity of the tumor microvasculature shown on DCE-MRI parametric maps and help capture the subtle variations in the texture which cannot be assessed by the single-resolution approach. However, this observation needs to be validated with a larger patient cohort. Consistent with the studies reporting mean parameter changes (9, 14–22), the results from this study provide further proof that changes in vascular perfusion/permeability represented by DCE-MRI imaging biomarkers are important features in identifying responders and nonresponders at the early stage of NACT.
On inspecting the AUC values from Figure 3, it can be observed that single-resolution fractal features performed consistently well in prediction of response for both the training and testing sets, although not as well as the multiresolution approach. The higher AUCs for fractal-based features suggest that they provide a richer representation of the heterogeneity in the tumor when compared to GLCM and RLM methods. The low dimensionality (d = 1) of single-resolution fractal feature is less likely to cause overfitting for the small sample size of our data set and therefore could lead to a good discriminative model. This could be one of the reasons that contributed to its effectiveness. On the other hand, in spite of increased dimensionality (d = 32), multiresolution fractal features exhibited better predictive performance, suggesting that analyzing heterogeneity at multiple resolutions provides a more comprehensive measure of the texture and thus increases the discriminative power of the feature. At each level of decomposition, the approximation coefficient [Wφ from equation (1)] represents the low-frequency component, which characterizes the coarse structure of the data, and the detail coefficients [Wψi from equation (2)] represent the high-frequency components, which capture the discontinuities and singularities in the data. Therefore, combination of features from different scales and frequencies gives a richer representation of the overall underlying texture. The advantages of multiresolution fractals can be expected to be even more significant when the data set is large enough to offset their high dimensionality.
The tumor heterogeneity appears to be captured well at the first level of decomposition as shown by Figure 4. Each decomposition level analyzes the signal at a particular band of frequencies. Higher decomposition levels have better frequency resolution. The first level of decomposition encompasses the entire frequency band of the input data in its subvolumes. Thereafter, we select the subvolume with the highest FD and perform multiresolution fractal analysis on that sub-band alone. By doing this we are effectively looking at finer frequency resolutions of the selected sub-band frequencies alone. Considering features from these finer frequency resolutions in isolation do not appear to have as much discriminative power as the first level of decomposition, but combining the finer frequency resolutions features with the features from first level appears to enrich the representation and provide incremental improvement.
As shown in Table 2 for the test data set, at fixed sensitivity, the higher AUC values from the multiresolution fractal features generally resulted in higher specificity values compared to those from other features. It is important to have high sensitivities so that most pCR patients will be correctly identified and continue with the original or de-escalated NACT regimen. At 80% sensitivity, the >60% specificity (except for the kep features) of the multiresolution fractal features implies that were this method used in clinical care, more than half of the non-pCRs would be correctly classified after the first NACT cycle, potentially enabling alteration of treatment plans for these nonresponders at the early stage of NACT to receive more personalized care.
This study has several limitations. The first being the small size of the data set used. The preliminary results obtained need to be evaluated on a larger patient cohort. Also due to the small size of the data set, dimensionality increase in feature vectors impedes the performance of the classifier. Larger data set can enable the choice of a richer feature vector from different levels in the multiresolution fractal decomposition, which might consistently outperform the other features. Finally, the DCE-MRI parametric maps used for feature analysis were obtained with the shutter-speed model, which is not commonly used in pharmacokinetic analysis of DCE-MRI data. In future studies, parametric maps obtained with the widely used standard Tofts model (47, 48), which generates only the Ktrans, ve, and kep parameters and thus results in reduced dimensionality of the feature vector, will be used for feature extractions and the results will be compared with those presented here.
In this preliminary study, we have demonstrated that multiresolution fractal analysis of voxel-based DCE-MRI parametric maps could be a promising tool for early prediction of breast cancer response to NACT. The multiresolution fractal features generally have better predictive performances than those extracted with the more conventional methods of GLCM, RLM, and single-resolution fractal analysis. Furthermore, compared to features extracted from individual DCE-MRI parametric maps, the use of concatenated features from all DCE-MRI parameters generally further improves prediction of NACT response.