Next Article in Journal
Role of Colposcopy after Treatment for Cervical Intraepithelial Neoplasia
Next Article in Special Issue
Multiparametric MRI for Prostate Cancer Characterization: Combined Use of Radiomics Model with PI-RADS and Clinical Parameters
Previous Article in Journal
Integrating the Tumor Microenvironment into Cancer Therapy
Previous Article in Special Issue
Automated Classification of Significant Prostate Cancer on MRI: A Systematic Review on the Performance of Machine Learning Applications
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Cancers
Volume 12
Issue 6
10.3390/cancers12061682
cancers-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Multiparametric Analysis of Longitudinal Quantitative MRI Data to Identify Distinct Tumor Habitats in Preclinical Models of Breast Cancer

by
Anum K. Syed
1,
Jennifer G. Whisenant
2,
Stephanie L. Barnes
3,
Anna G. Sorace
4,5,6 and
Thomas E. Yankeelov
1,3,7,8,9,*
1
Department of Biomedical Engineering, The University of Texas at Austin, Austin, TX 78712, USA
2
Department of Medicine, Vanderbilt University Medical Center, Nashville, TN 37232, USA
3
Oden Institute for Computational Engineering and Sciences, The University of Texas at Austin, Austin, TX 78712, USA
4
Department of Biomedical Engineering, The University of Alabama at Birmingham, Birmingham, AL 35294, USA
5
Department of Radiology, The University of Alabama at Birmingham, Birmingham, AL 35294, USA
6
O’Neal Comprehensive Cancer Center, The University of Alabama at Birmingham, Birmingham, AL 35294, USA
7
Department of Diagnostic Medicine, The University of Texas at Austin, Austin, TX 78712, USA
8
Department of Oncology, The University of Texas at Austin, Austin, TX 78712, USA
9
Livestrong Cancer Institutes, The University of Texas at Austin, Austin, TX 78712, USA
*
Author to whom correspondence should be addressed.
Cancers 2020, 12(6), 1682; https://doi.org/10.3390/cancers12061682
Submission received: 19 May 2020 / Revised: 20 June 2020 / Accepted: 22 June 2020 / Published: 24 June 2020
(This article belongs to the Special Issue Radiomics and Cancers)
Browse Figures
Review Reports Versions Notes

Abstract

:
This study identifies physiological tumor habitats from quantitative magnetic resonance imaging (MRI) data and evaluates their alterations in response to therapy. Two models of breast cancer (BT-474 and MDA-MB-231) were imaged longitudinally with diffusion-weighted MRI and dynamic contrast-enhanced MRI to quantify tumor cellularity and vascularity, respectively, during treatment with trastuzumab or albumin-bound paclitaxel. Tumors were stained for anti-CD31, anti-Ki-67, and H&E. Imaging and histology data were clustered to identify tumor habitats and percent tumor volume (MRI) or area (histology) of each habitat was quantified. Histological habitats were correlated with MRI habitats. Clustering of both the MRI and histology data yielded three clusters: high-vascularity high-cellularity (HV-HC), low-vascularity high-cellularity (LV-HC), and low-vascularity low-cellularity (LV-LC). At day 4, BT-474 tumors treated with trastuzumab showed a decrease in LV-HC (p = 0.03) and increase in HV-HC (p = 0.03) percent tumor volume compared to control. MDA-MB-231 tumors treated with low-dose albumin-bound paclitaxel showed a longitudinal decrease in LV-HC percent tumor volume at day 3 (p = 0.01). Positive correlations were found between histological and imaging-derived habitats: HV-HC (BT-474: p = 0.03), LV-HC (MDA-MB-231: p = 0.04), LV-LC (BT-474: p = 0.04; MDA-MB-231: p < 0.01). Physiologically distinct tumor habitats associated with therapeutic response were identified with MRI and histology data in preclinical models of breast cancer.

1. Introduction

Intratumoral heterogeneity plays a fundamental role in tumor progression, treatment resistance, and disease recurrence [1]. Tumors contain phenotypically and genetically diverse cancer cell populations of potentially varying proliferation rates and therapeutic sensitivities [1,2]. Additionally, the tumor microenvironment is remodeled as cancer cells recruit nonmalignant stromal cell populations (e.g., immune cells, fibroblasts, and vascular endothelial cells), thereby creating spatial variation within the tumor [3,4]. The resulting spatial distribution of cancer and stromal cells leads to local alterations in architecture of the microenvironment and the development of physiologically distinct subregions [4,5,6]. For example, high cancer cell proliferation and metabolism results in areas of acidic pH, and disparate vessel recruitment leads to regional variation in oxygen concentration across the tumor [4,5]. These tumor subregions can then exert selective pressures onto the cancer cells, increasing the heterogeneity among the cancer cells as they adapt to their unique local microenvironments [6,7,8,9]. Increased intratumoral heterogeneity has been shown to be associated with poorer patient prognosis and presents a difficult challenge in the treatment of cancer. Specifically in breast cancer, heterogeneity of cancer cell populations, accompanied with non-uniform drug delivery, can yield a diverse cellular response to treatment and result in clonal selection and an altered genomic landscape over the course of therapy [10,11,12,13]. Additionally, heterogeneous vasculature can lead to the emergence of hypoxic niches in regions deprived of oxygen, which have been associated with breast cancer resistance to chemo- and radiation therapies [14,15]. Characterizing intratumoral heterogeneity and its spatiotemporal evolution in response to both targeted and non-targeted therapies could yield better understanding of a patient’s individual tumor physiology and provide guidance for personalized treatment.
Noninvasive imaging enables longitudinal quantification of physiological tissue properties through methods such as diffusion-weighted (DW) magnetic resonance imaging (MRI) and dynamic contrast-enhanced (DCE) MRI. DW-MRI depends on the mobility of water within biological tissue, allowing for noninvasive evaluation of tissue microstructure [16]. The apparent diffusion coefficient (ADC), calculated from DW-MRI data, has been correlated to tumor cellularity and cell membrane permeability [17], and increases in ADC values have been shown to be early markers of breast tumor therapeutic response [16,18]. DCE-MRI provides measures of vascular perfusion and permeability, as well as tissue volume fractions [19]. Pharmacokinetic parameters, Ktrans (i.e., the volume transfer constant) and kep (efflux constant), derived from DCE-MRI data, can be used to quantify vascular alterations in response to treatment and serve as early markers of treatment response [20,21]. These quantitative imaging techniques, quantified on a voxel-wise basis, can characterize and provide biological insights into intratumoral heterogeneity across the entire tumor volume.
Approaches to quantify intratumoral heterogeneity from imaging data include histogram and texture analysis [22]. Histogram analysis involves quantitative assessment of parameter map distributions, demonstrating value over whole-tumor summary statistics [23,24], and allowing for the quantification of longitudinal alterations in intratumoral heterogeneity [25,26]. However, this approach forgoes the spatial information provided in imaging data. To overcome the limitations of one-dimensional analyses, texture analysis provides measures of spatial heterogeneity that have demonstrated utility in the prediction of breast cancer treatment response [27], as well as the classification of benign and malignant lesions [28]. However, this method quantifies spatial variations across the whole tumor, assuming the tumor is heterogeneous, but “well-mixed” [6]. This differs from clinical image data, which has shown tumor heterogeneity presenting regionally as spatially distinct subregions [29,30,31]. In an effort to more comprehensively characterize intratumoral biological heterogeneity, recent work has investigated the segmentation of tumors into distinct subregions by clustering similar voxels from multiparametric image data [32,33,34,35]. This emerging approach has been designated as habitat imaging [36], as the imaging-derived subregions are representative of tumor “habitats”, akin to ecological habitats in which a particular species population resides.
In this contribution, we investigated the ability of multi-parametric, voxel-based characterizations of tumor heterogeneity from DW- and DCE-MRI data to spatially resolve physiologically distinct tumor habitats in preclinical models of breast cancer. Two xenograft models of breast cancer were used, one a model of human epidermal growth factor receptor 2 positive (HER2+) disease (BT-474), and the other of triple negative disease (MDA-MB-231). We then utilized the identified habitats to quantify the longitudinal alterations of the tumor microenvironment in response to clinically relevant therapeutics. Finally, using end-point histology data, we biologically validated the imaging-derived habitats, through correlation with the habitats identified in the histological images.

2. Results

2.1. Discovery of MRI Tumor Habitats

Figure 1 presents the study design for the discovery of tumor habitats from quantitative MRI data and the analysis of longitudinal tumor composition defined by these habitats. The clustering analysis of quantitative imaging data yielded three primary clusters, or habitats, for both the BT-474 and MDA-MB-231 tumor cohort. The resultant dendrogram from agglomerative clustering (Figure S1), along with gap statistic plots were used to select a cluster number (k) of 3 (among two to 10 clusters) for both tumor cohorts.
The mean parameter values (Ktrans, kep, ADC, and ve) for each cluster were then evaluated to determine the physiological habitat each cluster represented. Statistically significant differences were found between the three identified clusters for all parameter values, in both the BT-474 (Figure 2a) and MDA-MB-231 (Figure 2b) tumor cohort, indicating that the agglomerative clustering algorithm was robust in identifying physiologically-distinct tumor habitats across the two preclinical models of breast cancer.
To identify the physiology each cluster represented, the quantitative measures of vascularity, Ktrans and kep, and cellularity, ADC and ve, were considered. For the BT-474 cohort (Figure 2a), the first identified cluster had a mean Ktrans of 0.18 min−1 and a mean kep of 0.84 min−1, and was thus labeled as a “high-vascularity” (HV) habitat. In comparison, the second and third clusters had mean Ktrans values of 0.05 and 0.04 min−1, respectively, and mean kep values of 0.31 and 0.08 min−1, respectively, and were labeled as “low-vascularity” (LV) habitats. Similar differences between clusters were observed in the MDA-MB-231 cohort (Figure 2b). The first cluster, labeled a “high-vascularity” (HV) habitat, had a mean Ktrans of 0.13 and a mean kep of 0.77, while the second and third clusters had mean Ktrans values of 0.04 and 0.01, respectively, and mean kep values of 0.32 and 0.01, respectively, and correspondingly were labeled as “low-vascularity” (LV) habitats.
Next, quantitative measures of cellularity, ADC and ve, were considered. For the BT-474 cohort (Figure 2a), the third identified cluster had a mean ADC of 1.2 × 10−3 mm2/s and a mean ve of 0.78 and was thus labeled as a “low-cellularity” (LC) habitat. In comparison, the first and second clusters had mean ADC values of 6.9 × 10−4 and 6.2 × 10−4 mm2/s, respectively, and mean ve values of 0.27 and 0.18, respectively, and were labeled as “high-cellularity” (HC) habitats. Similar differences between cluster values were observed in the MDA-MB-231 cohort (Figure 2b) for the ve parameter. The third cluster had a mean ve of 0.99, while the first and second clusters had mean ve values of 0.27 and 0.17, respectively. While each cluster showed statistically significant differences in mean ADC for the MDA-MB-231 cohort, all three clusters had ADC values within the “high-cellularity” range (i.e., less than 1 × 10−3 mm2/s, See Materials & Methods) with values of 6.9 × 10−4, 5.1 × 10−4, and 5.7 × 10−4 mm2/s for the first, second, and third cluster, respectively. Taken with ve measures of cellularity, the third cluster was labeled as a “low-cellularity” (LC) habitat, while the first and second clusters were labeled as “high-cellularity” (HC) habitats. Altogether, the three identified clusters from both tumor cohorts represented the following physiological habitats: high-vascularity high-cellularity (HV-HC), low-vascularity high-cellularity (LV-HC), and low-vascularity low-cellularity (LV-LC).
Figure 2c shows a representative habitat map and the corresponding quantitative parameter maps for a BT-474 tumor at day 0. Although no spatial information was used in the clustering process, the identified clusters were observed to colocalize spatially, demonstrating that we were able to spatially resolve physiologically distinct tumor habitats. These qualitative observations of spatial colocalization were confirmed quantitatively through multiregional spatial interaction (MSI) matrix analysis (Figure S2d). Habitat voxels mapped to their original spatial index showed statistically higher MSI values (p < 0.01) than if mapped to a random index, confirming that the identified habitats colocalize spatially.

2.2. Longitudinal Alterations in Tumor Composition in Response to Therapy

Figure 3a shows representative tumor habitat maps for a BT-474 control and trastuzumab-treated tumor, on days 0, 1, and 4. Prior to treatment, at day 0, no significant differences (p > 0.05) were observed in tumor composition between treated and control tumors, with mean percent tumor volumes (and 95 % confidence interval) of 43.5 ± 10.9% and 29.3 ± 17.1%, for the HV-HC habitat, 32.1 ± 7.7% and 41.8 ± 11.1%, for the LV-HC habitat, and 24.4 ± 6.3% and 28.9 ± 13.6% for the LV-LC habitat, respectively (Figure 3b). At day 4, trastuzumab-treated tumors showed a significant increase in percent tumor volume of the HV-HC habitat compared to baseline, with a mean value of 66.9 ± 24.3% (p = 0.02), and a significant decrease in percent tumor volume of the LV-HC habitat compared to baseline, with a mean value of 9.6 ± 8.6% (p < 0.01). For control tumors, no significant changes from baseline in percent tumor volume of any habitat were observed over the course of the study (p > 0.05). At day 4, trastuzumab-treated tumors were observed to have significantly higher HV-HC percent tumor volume compared to control tumors (33.3 ± 19.3%, p = 0.03), and significantly lower LV-HC percent tumor volume compared to control (47.8 ± 29.4%, p = 0.03). No significant differences were observed for the LV-LC habitat between treated and control tumors (p > 0.05).
Figure 3c shows representative tumor habitat maps for an MDA-MB-231 control and low-dose albumin-bound paclitaxel (ABP) treated tumor, at day 0, 1, and 3. Prior to treatment, at day 0, no significant differences (p > 0.05) were observed in tumor composition between low-dose ABP, high-dose ABP, and control tumors, with mean percent tumor volumes of 19.0 ± 6.4%, 25.0 ± 10.5%, and 26.1 ± 9.7% for the HV-HC habitat, 60.5 ± 5.7%, 54.3 ± 7.1%, and 57.1 ± 9.6% for the LV-HC habitat, and 20.5 ± 4.4%, 20.6 ± 4.6%, and 16.8 ± 5.7% for the LV-LC habitat, respectively (Figure 3d). Low-dose ABP-treated tumors showed a significant decrease from baseline in mean percent tumor volume of the LV-HC habitat at day 1 (45.2 ± 8.32%, p < 0.01) and at day 3 (49.5 ± 7.1%, p = 0.01). A corresponding nonsignificant increase in the percent tumor volume of the HV-HC habitat was observed at day 1 (32.8 ± 14.6%, p = 0.12) and at day 3 (27.4 ± 11.0%, p = 0.15). Control and high-dose ABP-treated tumors did not show any significant changes from baseline in percent tumor volume of any habitat over the course of the study (p > 0.05). Additionally, no significant differences were observed between treatment groups for any habitat over the course of the study (p > 0.05).

2.3. Discovery of Histological Tumor Habitats

Figure 1 shows the overall study design for the discovery of tumor habitats from histological data, through the use of stain density maps, and the biological validation of MRI-derived habitats using histological habitats. The clustering analysis of histological stain density maps yielded three primary clusters, or habitats, for both the BT-474 and MDA-MB-231 tumor cohort. For each cohort, a k = 3 showed the optimal gap statistic among 2 to 10 clusters.
The mean vessel count and nuclei count for each cluster was then evaluated to determine the physiological habitat each cluster represented. Statistically significant differences were found between the three identified clusters for all stain density parameter values, in the BT-474 (Figure 4a) and MDA-MB-231 (Figure 4b) tumor cohort. For both tumor cohorts, the first identified cluster was labeled as “high-vascularity” (HV), with mean vessel counts/unit (in this context, “unit” refers to the 150 μm × 150 μm neighborhood around a pixel) of 5.5 for BT-474 cohort and 3.9 for MDA-MB-231 cohort. The second and third clusters showed mean vessel values of 2.2 and 0.7 counts/unit for the BT-474 cohort, and 0.58 and 0.28 counts/unit for the MDA-MB-231 cohort, respectively, and were labeled as “low-vascularity” (LV) habitats. For both tumor cohorts, the third identified cluster was labeled as “low-cellularity” (LC), with mean nuclei counts/unit of 17.8 for the BT-474 cohort and 60.1 for the MDA-MB-231 cohort. The first and second clusters showed mean nuclei values of 142.3 and 109.5 counts/unit, for the BT-474 cohort, and 249.9 and 220.8 counts/unit for the MDA-MB-231 cohort, respectively, and were labeled as “high-cellularity” (HC) habitats. Altogether, the three identified clusters from both tumor cohorts represented the following histological habitats: high-vascularity high-cellularity (HV-HC), low-vascularity high-cellularity (LV-HC), and low-vascularity low-cellularity (LV-LC).
Figure 4c shows a representative habitat map and the corresponding stain density maps for an example MDA-MB-231 tumor, from the low-dose ABP treatment group. As with the discovery of MRI habitats, while no spatial information was used in the clustering process to identify histological habitats, the identified clusters were observed to colocalize spatially, which was confirmed quantitatively through MSI matrix analysis (Figure S3).
To further investigate the underlying physiologies that the identified histological habitats represent, we evaluated the necrosis and Ki-67 expression levels for each habitat. The LV-LC histological habitat showed high levels of necrosis (Figure 4), with a mean necrosis value of 48.1 percent area/unit for the BT-474 cohort and 21.1 percent area/unit for the MDA-MB-231 cohort. In comparison, the high cellularity habitats, HV-HC and LV-HC, showed low levels of necrosis with mean necrosis values of 9.4 and 10.9 percent area/unit for the BT-474 cohort, and 4.3 and 3.6 percent area/unit for the MDA-MB-231 cohort, respectively. The HV-HC histological habitat showed high levels of Ki-67 expression (Figure 4), with mean Ki-67+ nuclei values of 31.0 counts/unit for the BT-474 cohort and 23.7 counts/unit for the MDA-MB-231 cohort. In comparison, the low vascularity habitats, LV-HC and LV-LC, showed low levels of Ki-67 expression with mean Ki-67+ nuclei values of 6.6 and 0.5 counts/unit for the BT-474 cohort, and 8.5 and 0.9 counts/unit for the MDA-MB-231 cohort, respectively.

2.4. Correlation between MRI-Derived and Histology-Derived Habitats

A positive linear correlation was observed between MRI percent tumor volume and histological percent tumor area of HV-HC habitats (Figure 5) for the BT-474 tumor cohort, with an r2 = 0.74 (p = 0.03), and a nonsignificant positive correlation for the MDA-MB-231 cohort, with an r2 = 0.23 (p > 0.05). A positive linear correlation was also observed between MRI percent tumor volume and histological percent tumor area of LV-LC habitats (Figure 5) for both the BT-474 tumor cohort, with an r2 = 0.70 (p = 0.04), and the MDA-MB-231 cohort, with an r2 = 0.73 (p < 0.01). Correlation analysis of the LV-HC habitats showed a nonsignificant positive correlation with an r2 = 0.55 (p = 0.09), for the BT-474 cohort, and a significant positive correlation with an r2 = 0.53 (p = 0.04), for the MDA-MB-231 cohort.

3. Discussion

In this contribution, we utilized quantitative imaging parameters derived from DW- and DCE-MRI to spatially resolve tumor habitats in two preclinical models of breast cancer. Three habitats were identified and labeled based on imaging measures of cellularity and vascularity as HV-HC, LV-HC, and LV-LC. Although no spatial information was used in the clustering process, identified habitats were found to colocalize spatially (Figure S2d), demonstrating that the clustering of quantitative imaging parameters alone was robust to identifying tumor habitats.
Interestingly, the HV-HC subregions were often observed to localize towards the outer edge of the tumor (Figure 2c), characteristic of the “enhancing rim” phenotype often empirically noted in tumors. Conversely, the LV-LC subregions were observed to localize towards the center of the tumor, possibly representing core necrosis frequently found in tumors. These habitat physiologies identified with MRI were biologically validated by comparing to clusters identified in corresponding histology data (Figure 5). Independent clustering of stain density maps yielded three histological habitats analogous to those identified from MRI data (HV-HC, LV-HC, and LV-LC), lending support to the selection of a k = 3 in the MRI habitat analysis. In both the BT-474 and MDA-MB-231 cohorts, the percent tumor area for all three histological habitats showed a positive linear correlation with the percent tumor volume of their corresponding MRI-derived habitats, (Figure 5b). This correlation indicates that tumor habitats derived from MRI data represent biologically distinct tumor subregions, defined in terms of cellularity and vascularity. Further analysis of histological habitats showed higher counts of Ki-67 positive nuclei in the HV-HC habitat compared to LV-HC and LV-LC habitats (Figure 4), suggesting that this habitat may represent viable tumor regions. These results align with findings by others of higher Ki-67 expression in regions of vascular perfusion at the tumor periphery [37,38]. Additionally, the LV-LC histological habitat showed significantly higher levels of necrosis compared to the HV-HC and LV-HC habitats (Figure 4), indicating this habitat may represent regions of necrosis, an observation supported by the literature [37]. It is noteworthy that habitat trends in mean parameter values (e.g., Ktrans of HV-HC > Ktrans of LV-HC > Ktrans of LV-LC) were consistent between the two tumor models, with the exception of ADC for the LV-LC habitat. The LV-LC habitat showed increased mean ADC, typical of necrosis [16,17], compared to other habitats for the BT-474 cohort (Figure 2a). However, the same was not observed for the MDA-MB-231 cohort (Figure 2b), a finding which was surprising since the imaged LV-LC habitat showed a significant correlation to the histological LV-LC (necrosis) habitat (Figure 5b). The LV-LC mean ADC for the MDA-MB-231 cohort (5.7 x 10−4 mm2/s) was similar to that observed for regions of necrosis in other studies of the same tumor model [39,40]. Differences in BT-474 and MDA-MB-231 necrosis morphology may attribute to differences observed in ADC values of the LV-LC habitat.
The longitudinal evaluation of tumor habitat composition was facilitated by the use of quantitative imaging parameters, allowing the pooling of voxel data across imaging time points for each tumor cohort prior to clustering analysis. Clustering cohort voxel data, as opposed to individual tumor data, enabled the habitat analysis to be robust to inter-tumoral variations, avoiding the assumption that all habitats are present in each tumor, and ensure the habitat definitions were consistent across time. Prior to treatment, no difference in tumor composition was observed between treatment groups in either cohort (Figure 3). However, tumor composition at baseline differed between the HER2+ and triple negative tumors, with the BT-474 tumor cohort showing a higher percent tumor volume of the HV-HC habitat and lower percent tumor volume of the LV-HC habitat compared to MDA-MB231 (Figure 3b,d and Figure S4). The differences between subtypes in baseline tumor composition may be due to HER2’s involvement in VEGF expression and vascular recruitment [41], leading to potentially different vascular patterns between these two tumor models. Several radiomics studies have demonstrated the utility of clinical image features to predict molecular subtypes of breast cancer [42,43,44,45]. The observed differences in tumor composition at baseline indicate the potential to discriminate breast cancer subtypes using identified habitats.
In this study, habitats of therapeutic response were identified for both targeted and chemotherapies. At day 4 post-treatment with the anti-HER2 targeted therapy, trastuzumab, BT-474 tumors showed a significant decrease in the percent tumor volume of the LV-HC habitat (p < 0.01) and corresponding increase in the percent tumor volume of the HV-HC habitat (p = 0.02) from baseline. The LV-HC habitat represents a subregion with decreased vascular perfusion and higher cellularity, potentially corresponding to hypoxic regions with continued cell proliferation. The decrease in percent tumor volume of the LV-HC habitat upon trastuzumab treatment is concordant with findings from previous work, where trastuzumab therapy led to a longitudinal decrease in hypoxia in BT-474 tumors, quantified by 18F-FMISO positron emission tomography (PET) and pimonidazole staining of histological sections [46]. Observed increases in percent tumor volume of the HV-HC habitat and decreases in the LV-HC habitat upon trastuzumab treatment are supported by other studies demonstrating improved vascular perfusion resulting from trastuzumab therapy [47,48,49,50,51]. MDA-MB-231 tumors were treated with a low or high dose of albumin-bound paclitaxel (ABP), an established taxane cell-cycle inhibitor [52]. At day 3, the low dose of ABP therapy showed significant decreases in percent tumor volume of the LV-HC habitat (p = 0.01) compared to baseline values and corresponding increases (p > 0.05) in percent tumor volume of the HV-HC habitat. Paclitaxel has been shown to have anti-angiogenic effects [53,54], and its albumin-bound form has been shown to induce vascular remodeling [55]. Decreases in the LV-HC habitat upon low-dose ABP therapy may indicate improved vascular perfusion from the vascular remodeling effects of ABP therapy.
Regions of low vascular perfusion and high cellularity, as with the LV-HC habitat, may potentially represent habitats where cells have adapted to low oxygen and nutrient conditions, resisting apoptosis and continuing to proliferate [6]. Thus, identifying these cell-dense regions of hypoxia may be of clinical interest as they could represent regions where underlying cell populations are more likely to be resistant to therapeutic interventions [5,15]. The LV-HC habitat may represent regions deprived of oxygen, although hypoxia was not directly measured in this study. However, several groups have evaluated the relationship between MRI-derived habitats and hypoxia, and found regions of moderate perfusion to often correspond to hypoxia [37,39,56]. Using spatially coregistered histology and multiparameteric MRI data, Jardim-Perassi et al. found MRI habitats corresponding to viable-hypoxic tissue showed decreased area under the curve values calculated from DCE-MRI, compared to viable-normoxic habitats [39]. Others have clustered DCE-MRI signal intensity time courses to identify subregions of different contrast agent uptake patterns and related identified subregions to tissue oxygenation. Stoyanova et al. found regions with delayed contrast enhancement and wash out to correspond to hypoxia, as measured by pimonidazole staining and 18F-FMISO-PET imaging [57]. Regions of similar enhancement patterns were identified by Featherstone et al., who found such regions to correspond with decreased oxygen enhancement as measured by oxygen-enhanced MRI [33].
Often quantitative imaging-based evaluations of treatment response involve summarizing an entire tumor volume down to a single statistic, largely ignoring the spatial heterogeneity within a tumor. In this study, we demonstrated the utility of habitat imaging to evaluate longitudinal alterations in tumor heterogeneity through the identification of tumor habitats of response to both targeted and non-targeted treatments. Other groups have also investigated the application of quantitative MRI parameter maps to identify tumor subregions of response. Quantitative ADC and T2 maps have been used to segment tumor subregions of viable tissue and necrosis, with research showing significant changes in viable tumor fractions in response to anti-angiogenic [58] and radiation [59] therapy in preclinical tumor models. Longo et al. found analysis of tumor subregions derived from quantitative DCE-MRI parameters outperformed whole tumor analysis, when evaluating response to antiangiogenic therapy in a murine model of breast cancer [60]. These studies, accompanying our own, demonstrate the potential application of quantitative habitat imaging to measure the spatiotemporal evolution of intratumoral heterogeneity in response to therapy.
One limitation to our study was the lack of spatial registration of ex vivo histological habitats to those identified from in vivo imaging data. While such registration could provide a more direct biological validation to our in vivo findings, tissue section alignment with imaged slices, along with tissue shrinkage or deformation post-tumor excision, presents a difficult challenge for spatial registration of histology data to MRI [61]. While we were unable to spatially correlate imaged habitats with those derived from histology, we found significant correlations between tumor composition of histology and central-slice MRI data to lend biological support for MRI-derived habitats. Another important limitation is that identified habitats are constrained to the imaging parameters used to identify them, as well as the spatial resolution of the MRI acquisition. If additional imaging measures were employed, such as those sensitive to pH, hypoxia, or mechanical stiffness, it is possible that additional physiologically distinct habitats may be discovered. The spatial resolution of MRI also influences the extent of partial volume effects, as multiple habitats may be averaged within a tumor voxel [39]. Employing soft clustering techniques that provide a probabilistic cluster assignment to each voxel, such as gaussian mixture modeling, may be more robust to multiple habitats within a voxel and address these limitations. Finally, the imaging parameter thresholds used to label the MRI habitats as high or low cellularity or vascularity were selected manually, and future efforts will involve more systematic methods to determine these thresholds.
To apply these methods in the clinical setting, several challenges must be identified and then overcome. These include the limited signal-to-noise ratio and spatial and temporal resolution, as well as the presence of image artifacts (due to, for example, patient motion or biopsy clips), and the ability to obtain high-quality quantitative DW- and DCE-MRI parametric maps. The limited spatial resolution available in clinical MRI, may particularly confound this habitat imaging approach, as it is predicated upon being able to resolve and quantify intratumoral heterogeneity. Fortunately, clinical studies have shown that DW- and DCE-MRI metrics can be acquired with repeatability and reproducibility in the clinical setting [62,63,64], providing evidence of feasibility for the future application of our technology to clinical studies. In particular, the individual quantitative parameters derived from DW- and DCE-MRI have been identified as potential biomarkers of (for example) therapeutic response for breast cancer patients undergoing neoadjuvant therapy [16,19,20,65], thereby providing support for potential clinical translation of our habitat imaging approach. Furthermore, habitat imaging has already made significant contributions in other disease locations, such as in the brain [66,67] and prostate [68,69]; thus, it is a reasonable hypothesis that these methods may be applied for the predicting response in breast cancer patients.

4. Materials and Methods

4.1. Cell Culture & Animal Model

Two xenograft tumor models were employed using the human-derived BT-474 (HER2+ breast cancer) and MDA-MB-231 (triple negative breast cancer) breast cancer cell lines (ATCC, Manassas, VA, USA). BT-474 cells were cultured in improved minimal essential medium (Invitrogen, Carlsbad, CA, USA) supplemented with 10% fetal bovine serum and 1% insulin at 37 °C in 5% CO2. MDA-MB-231 cells were cultured in Dulbecco’s modified essential medium supplemented with 10% fetal bovine serum at 37 °C in 5% CO2.
Data from two independent cohorts of mice as part of previous studies were used in this study [50,70]. Cohort one of female athymic nude mice were subcutaneously implanted with a 0.72 mg, 60-day release, 17β-estradiol pellet (Innovative Research of America, Sarasota, FL, USA) 24 h prior to subcutaneous injection of BT-474 cells. The second cohort of female athymic nude mice were subcutaneously injected with MDA-MB-231 cells. Tumors were grown until they reached approximately 275 mm3 in volume, when they were then randomly assigned to treatment groups. The BT-474 xenograft tumor cohort (N = 20) were either treated with trastuzumab (10 mg/kg) or saline control. The MDA-MB-231 xenograft tumor cohort (N = 37) were either treated with 15 mg/kg (low-dose) albumin-bound paclitaxel (ABP), 25 mg/kg (high-dose) ABP, or saline control. Treatment and imaging schedules for both cohorts of mice are shown in Figure 6. A flow chart of sample number for each treatment group, over the course of the study, is shown in Figure S5.

4.2. Magnetic Resonance Imaging and Analysis

While details are provided elsewhere [50,51,70], here we briefly summarize the salient MRI procedures. Mice were imaged using a 7T small animal MRI scanner (Agilent Technologies, Palo Alto, CA, USA) to acquire high resolution T2-weighted anatomical images, DW- and DCE-MRI. The day prior to the baseline imaging, a jugular catheter was implanted into mice for exogenous delivery of the contrast agent for DCE-MRI. All MRI data were acquired with transaxial slices, 1 mm slice thickness, and a 28 × 28 × 15 mm field of view. High resolution T2-weighted anatomical images were collected for tumor segmentation using a 128 × 128 × 15 matrix and TR/TE/NEX = 5500 ms/ 35.65 ms/2. Tumor regions-of-interest (ROIs) were manually drawn around the tumor boundary for all slices containing the tumor using the high resolution T2-weighted images and custom-built MATLAB scripts (MATLAB, MathWorks, Natick, MA, USA). ROIs were modified manually to exclude regions where no T1-weighted or T2-weighted signal was observed, suggestive of calcifications [71]; see Figure S6.
Diffusion-weighted data was collected using a standard pulse gradient spin echo sequence with three b-values (150, 500, and 800 s/mm2) and the following scan parameters: 64 × 64 matrix, TR/TE/NEX = 2000 ms/30ms/2; gradient duration and intervals of δ = 3 ms and Δ = 20 ms, respectively. The signal intensity data from acquired DW-MRI was fit voxel-wise (MATLAB) to Equation (1) to calculate ADC maps via:
S b = S 0 exp − A D C × b ,
where S(b) is the signal intensity after gradient application, S0 is the signal intensity prior to gradient application, and b is the strength of the diffusion gradient. Voxels within the ROI were removed from analysis if the signal intensities at the three b-values did not decrease monotonically; i.e., we require that S(b = 800) < S(b = 500) < S(b = 150). The average (and 95% confidence interval) percent of tumor voxels removed was 7.31% (1.30%) per DWI scan. Slices of the DW-MRI data showing motion artifacts were removed from subsequent analysis if greater than 20% of the slice’s voxel were removed via the aforementioned criteria. DW-MRI data for a tumor on a given imaging time point was discarded entirely if more than half of the slices were removed from analysis due to the aforementioned criteria (BT-474 cohort: 6 out of 56 scans, MDA-MB-231 cohort: 3 out of 110 scans, Figure S5).
After DW-MRI acquisition, pre-contrast T1 maps were acquired using a multislice inversion recovery snapshot fast low angle shot gradient echo sequence. Seven inversion times ranging from 250 to 10,000 ms were used in addition to the following scan parameters: 64 × 64 matrix, TR/TE/NEX = 12,000 ms/2.1 ms/2, and a flip angle of α = 15 °. T1 maps were calculated by fitting voxel data (MATLAB) to Equation (2):
S T I = S 0 1 − 2 exp − T I T 1 + exp − T R T I ,
where S(TI) is the signal intensity at inversion time TI and S0 is the baseline signal intensity.
DCE-MRI acquisition involved the collection of dynamic T1-weighted images using a spoiled gradient echo sequence with a temporal resolution of 12.8 s for 20 minutes, with the following scan parameters: 64 × 64 matrix, TR/TE/NEX = 100 ms/2.1 ms/2, and a flip angle of α = 25 °. Pre-contrast images were acquired for approximately two minutes before a bolus of 0.05 mmol/kg gadolinium-diethylenetriaminepentaacetic acid (Magnevist, Bayer, Whippany, NJ) was delivered via a jugular catheter. DCE-MRI data were voxel-wise fit to the standard Kety-Tofts model (Equation (3) [72,73]:
C t t = K t r a n s ∫ 0 t C p u exp − K t r a n s v e t − u d u ,
using a population-derived [74] arterial input function (AIF) [74,75], Cp(t), to estimate Ktrans, and ve (i.e., the extravascular extracellular volume fraction). To reduce systematic errors in the model fits, the AIF was calibrated for each DCE-MRI scan to segmented skeletal muscle tissue [76] prior to fitting for pharmacokinetic parameters. Additionally, each tumor voxel time-course was fit to a half-logistic function to determine the bolus arrival-time, which was used to shift the calibrated AIF on a voxel-wise basis to align the enhancement time of Cp(t) with the bolus arrival-time [77]. Voxel-wise fits to the standard Kety-Tofts model were bound within the physiological range for each parameter (e.g., 0 to 5 min−1 for Ktrans, 0 to 1 for ve). The efflux constant, kep, was calculated as Ktrans/ve, with values beyond the physiological range (i.e., if 0 min−1 < kep < 5 min−1) removed from analysis. The average (and 95% confidence interval) percent of tumor voxels removed was 0.76% (0.15%) per DCE scan. DCE-MRI data was not used if contrast agent delivery was unsuccessful (BT-474 cohort: four out of 56 scans, MDA-MB-231 cohort: 12 out of 110 scans, Figure S5).

4.3. Discovery of MRI Tumor Habitats

Multiparametric voxel data was extracted for tumor ROIs yielding a four-dimensional vector whose entries consist of Ktrans, ve, kep, and ADC for each voxel. For each tumor cohort, voxel data was pooled across all imaging time points and each parameter distribution was scaled to have a mean of zero and standard deviation of one, to allow all parameters to contribute equally in the clustering process [78]. Tumor habitats were then identified through automated clustering of voxel data using an agglomerative clustering [79] algorithm with the Ward linkage and Euclidean distance measure (RStudio, “Rclusterpp” package). Details of agglomerative clustering are described in the supplemental methods. In this study, both the dendrogram branch height and gap statistic [80] were used to determine the number of clusters (RStudio, “cluster” package). No spatial information was used in the clustering process.
Cluster-labeled voxels were then mapped back to their original spatial index and a multiregional spatial interaction (MSI) matrix (adapted from Wu et al. [81]) was used to quantify spatial colocalization of cluster-labeled voxels. Details of MSI matrix analysis are described in supplemental methods. To identify what type of physiology each cluster (habitat) represented, the mean value for each MRI parameter was calculated to label each identified cluster in terms of high or low “vascularity” or “cellularity”. Clusters with mean Ktrans > 0.1 min−1 or kep values > 0.5 min−1 were designated as high-vascularity (HV) habitats, and otherwise as low-vascularity (LV) habitats. Clusters with mean ADC < 1 × 10−3 mm2/s or mean ve < 0.3 were designated as high-cellularity (HC) habitats, and otherwise as low-cellularity (LC) habitats.

4.4. Quantifying Longitudinal Alterations in Tumor Composition

Longitudinal alterations in tumor heterogeneity were quantified by changes in the percent of tumor volume comprised by each habitat to determine if habitat tumor composition changes with respect to a therapeutic intervention. The percent tumor volume of each habitat was averaged within each treatment group and evaluated at every time point for each tumor model cohort. Tumors that were larger than 500 mm3 at day 0 were excluded from treatment response analysis (BT-474: N = 1, MDA-MB-231: N = 2, Figure S5).

4.5. Immunohistochemistry Staining and Image Processing

On the final imaging time point, (i.e., day 3 for the MDA-MB-231 cohort, and day 4 for the BT-474 cohort) each animal was sacrificed and the tumor excised for histological analysis (BT-474 cohort: N = 16, MDA-MB-231 cohort: N = 19, Figure S5). Tumors were sliced at the largest cross-section corresponding (as best as possible) to the central in vivo imaging plane and fixed in 10% formalin for 48 h, then stored in 70% ethanol for processing. Upon processing, tissue samples were embedded in paraffin and sliced into 5 μm-thick sections. Sections were mounted onto glass slides and stained with hemotoxylin and eosin (H&E), anti-CD31, or anti-Ki-67 as previously described [51,70]. A Leica SCN400 (Leica Biosystems, Nussloch, Germany) slide scanner was used to digitally scan slides in high-resolution (20×) brightfield. Images were processed using Bio-Formats [82] (http://openmicroscopy.org) and custom-built MATLAB scripts.
Images of histology sections were then automatically segmented to yield stain-positive masks across the entire tumor slice (Figure S7a,b). A k-means clustering algorithm was used for automated segmentation of the H&E stained sections to generate masks for nuclei and necrosis across the entire tumor slice. Otsu thresholding [83] was used for automated segmentation of proliferative cell nuclei from anti-Ki-67 stained sections and microvessels from anti-CD31 stained sections. Further details of the automated segmentation techniques are described in previous work [26]. Masks of whole tissue for each stained tissue slice were semi-automatically generated using gray (non-tissue) area masks and manually drawn ROIs (Figure S7c). These whole tissue masks were used to spatially align different histological stains from separate slides to the corresponding H&E stain for each tumor sample, through a rigid registration algorithm (Figure S7d).
To quantify regional densities in cellularity and vascularity, a neighborhood analysis was performed to generate stain density maps for each tissue section. A 150 μm × 150 μm [84] neighborhood around each pixel was evaluated to facilitate comparison to the larger MRI voxels (437.5 μm × 437.5 μm) without losing the ability to assess regional variations of the tumor microenvironment. To quantify CD31+ microvessel density, H&E nuclear density, and Ki-67+ nuclear density, individual connected components within a 150 μm × 150 μm neighborhood around each pixel were counted across each tissue slice. Each neighborhood count was then assigned as that pixel’s nuclei (H&E, Ki-67) or microvessel (CD31) density value, thus yielding a neighborhood heatmap of stain density (Figure S7e). A similar neighborhood analysis was performed for H&E necrosis density, where the percent stain-positive area within a 150 μm × 150 μm neighborhood around each pixel was calculated and assigned as that pixel’s necrosis density value (Figure S7e). Stain density maps are presented in counts per unit (counts/unit) or percent area per unit (percent area/unit), where a “unit” represents the 150 μm × 150 μm neighborhood around each pixel.

4.6. Discovery of Histological Tumor Habitats

To reduce redundancy in the clustering analysis of stain density maps, each stain density map was downsampled by a factor of 0.01, using nearest neighbor interpolation (Figure S7f). Multiparametric pixel data from spatially-aligned stain density maps were extracted for the tissue-masked areas of each tumor, yielding a four-dimensional vector for each pixel (nuclei, necrosis, Ki-67+ nuclei, and vessel density). Pixel data was pooled within each tumor cohort (BT-474 or MDA-MB-231) and each parameter distribution was scaled to have a mean of zero and standard deviation of one, so that all parameters contribute equally in the clustering process [78]. Pixel data was then clustered (Figure S7g) through agglomerative clustering, in the same manner as MRI-derived habitats (described above), to identify physiologically-distinct histological subregions within tumor cross-sections. No spatial information was used in the clustering process. Tumors were included in the analysis only if all three (H&E, CD31, and Ki-67) stained sections were available (BT-474 cohort: N = 12, MDA-MB-231 cohort: N = 13, Figure S5).
Cluster-labeled pixels were then mapped back to their original spatial index and spatial colocalization of cluster-labeled pixels was measured through MSI matrix analysis, in the same manner as MRI-derived habitats. To identify what type of physiology each cluster (habitat) represented, the mean value for vessel and nuclei count was calculated to label each identified cluster in terms of high or low “vascularity” or “cellularity”, respectively. Clusters with mean vessel count/unit values greater than three were designated as high-vascularity (HV) habitats, and otherwise as low-vascularity (LV) habitats. Clusters with mean nuclei count/unit values greater than 100 were designated as high-cellularity (HC) habitats, and otherwise as low-cellularity (LC) habitats.

4.7. Correlations between MRI and Histological Habitats

To biologically validate MRI-derived habitats, a correlation analysis was performed between endpoint MRI and corresponding histology data. Tumors with both endpoint MRI data and histology data taken from the central tumor cross-section, were included in the analysis (BT-474 cohort N = 6, MDA-MB-231 cohort: N = 8, Figure S5). To identify central slices of the tumor ROI, high-resolution T2-weighted MRI and b = 800 DW-MRI data were visualized alongside corresponding H&E histology sections. Each of the MRI slices were stepped through and the anatomical similarity of the shape of the tumor slice as well as the presence and location of regions of necrosis were assessed to determine the closest possible match to the histology slice (Figure S8). The percent tumor volume of each habitat was calculated for the central slice of the tumor ROI at the study end point (day 4 for BT-474 cohort, day 3 for MDA-MB-231 cohort). These values were then correlated with the percent tumor area of each habitat across the histological section.

4.8. Statistical Analysis

To test if differences in mean parameter values between habitats (MRI-derived or histology-derived) were significant, a one-way analysis of variance (ANOVA) followed by Tukey’s honest significant difference test was used. Differences from baseline within treatment groups were compared using a nonparametric Wilcoxon rank sum test. Differences between treatment groups were compared using a nonparametric Wilcoxon rank sum test (BT-474) or one-way ANOVA followed by Tukey’s honest significant difference test (MDA-MB-231). Correlations were tested using Pearson’s product-moment correlation. In all statistical tests a p-value less than 0.05 was considered significant.

5. Conclusions

This study employed quantitative MRI data to spatially resolve physiologically-distinct tumor habitats of varying cellularity and vascular properties in two preclinical models of breast cancer. We identified tumor habitats associated with therapeutic response to targeted and chemotherapies. We also provided biological validation of imaged habitats through correlation analysis with histological habitats identified from stain density maps. Application of habitat imaging techniques utilizing physiological parameters derived from quantitative MRI could provide greater understanding of a patient’s individual tumor physiology and provide guidance for personalized cancer treatment.

Supplementary Materials

The following are available online at https://www.mdpi.com/2072-6694/12/6/1682/s1, Figure S1: Resultant dendrograms from clustering of voxel data, Figure S2: Quantifying spatial colocalization of MRI habitats, Figure S3: Quantifying spatial colocalization of histological habitats, Figure S4: Comparison of baseline tumor composition between BT-474 and MDA-MB-231 cohorts, Figure S5: Sample numbers across study, Figure S6: Example ROI segmentation, Figure S7: Image processing steps for histology data, Figure S8: Example central tumor ROI slice.

Author Contributions

Conceptualization, A.K.S., A.G.S., T.E.Y.; methodology, A.K.S., J.G.W., S.L.B., A.G.S., T.E.Y.; software, A.K.S.; validation, A.K.S., A.G.S., T.E.Y.; formal analysis, A.K.S.; investigation J.G.W., S.L.B., A.G.S.; resources, A.K.S., J.G.W., S.L.B., A.G.S., T.E.Y.; data curation, A.K.S.; writing—original draft preparation, A.K.S.; writing—review and editing, A.K.S., J.G.W., S.L.B., A.G.S., T.E.Y.; visualization, A.K.S.; supervision, A.G.S., T.E.Y.; project administration, A.G.S., T.E.Y.; funding acquisition, A.G.S., T.E.Y.; All authors have read and agreed to the published version of the manuscript.

Funding

We thank the Cancer Prevention and Research Institute of Texas for support through RR160005, and the American Cancer Society for support through RSG-18-006-01-CCE. We thank the National Cancer Institute for funding through R01CA138599, U01CA17470, 1U24 CA226110, and R01CA240589. We thank the National Institute of Biomedical Imaging and Bioengineering for support through T32EB007507.

Acknowledgments

T.E.Y. is a CPRIT Scholar in Cancer Research. We thank Alfiia Galimzianova, Imon Banerjee, Ryan Woodall, and Chengyue Wu for several informative discussions.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Marusyk, A.; Polyak, K. Tumor heterogeneity: Causes and consequences. Biochim. Biophys. Acta BBA Rev. Cancer 2010, 1805, 105–117. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  2. Greaves, M.; Maley, C.C. Clonal evolution in cancer. Nature 2012, 481, 306–313. [Google Scholar] [CrossRef] [PubMed]
  3. Junttila, M.R.; de Sauvage, F.J. Influence of tumour micro-environment heterogeneity on therapeutic response. Nature 2013, 501, 346–354. [Google Scholar] [CrossRef] [PubMed]
  4. Hanahan, D.; Coussens, L.M. Accessories to the crime: Functions of cells recruited to the tumor microenvironment. Cancer Cell 2012, 21, 309–322. [Google Scholar] [CrossRef] [Green Version]
  5. Koumenis, C.; Hammond, E.; Giaccia, A. (Eds.) Tumor Microenvironment and Cellular Stress; Advances in Experimental Medicine and Biology Series; Springer: New York, NY, USA, 2014; Volume 772, ISBN 978-1-4614-5914-9. [Google Scholar]
  6. Gatenby, R.A.; Grove, O.; Gillies, R.J. Quantitative imaging in cancer evolution and ecology. Radiology 2013, 269, 8–14. [Google Scholar] [CrossRef] [Green Version]
  7. Polyak, K.; Haviv, I.; Campbell, I.G. Co-evolution of tumor cells and their microenvironment. Trends Genet. 2009, 25, 30–38. [Google Scholar] [CrossRef]
  8. Tabassum, D.P.; Polyak, K. Tumorigenesis: It takes a village. Nat. Rev. Cancer 2015, 15, 473–483. [Google Scholar] [CrossRef]
  9. Merlo, L.M.F.; Pepper, J.W.; Reid, B.J.; Maley, C.C. Cancer as an evolutionary and ecological process. Nat. Rev. Cancer 2006, 6, 924–935. [Google Scholar] [CrossRef]
  10. Balko, J.M.; Giltnane, J.M.; Wang, K.; Schwarz, L.J.; Young, C.D.; Cook, R.S.; Owens, P.; Sanders, M.E.; Kuba, M.G.; Sánchez, V.; et al. Molecular profiling of the residual disease of triple-negative breast cancers after neoadjuvant chemotherapy identifies actionable therapeutic targets. Cancer Discov. 2014, 4, 232–245. [Google Scholar] [CrossRef] [Green Version]
  11. Sato, F.; Saji, S.; Toi, M. Genomic tumor evolution of breast cancer. Breast Cancer 2016, 23, 4–11. [Google Scholar] [CrossRef]
  12. Zardavas, D.; Irrthum, A.; Swanton, C.; Piccart, M. Clinical management of breast cancer heterogeneity. Nat. Rev. Clin. Oncol. 2015, 12, 381–394. [Google Scholar] [CrossRef]
  13. Koren, S.; Bentires-Alj, M. Breast tumor heterogeneity: Source of fitness, hurdle for therapy. Mol. Cell 2015, 60, 537–546. [Google Scholar] [CrossRef] [PubMed]
  14. Kunkler, I.H.; Ward, C.; Speers, C.; Pierce, L.; Feng, F.; Meehan, J.; Langdon, S.P. Personalisation of radiotherapy for breast cancer. In Personalized Treatment of Breast Cancer; Toi, M., Winer, E., Benson, J., Klimberg, S., Eds.; Springer Japan: Tokyo, Japan, 2016; pp. 131–149. ISBN 978-4-431-55551-3. [Google Scholar]
  15. Gilkes, D.M. Hypoxia and Cancer Metastasis; Springer: Cham, Switzerland, 2019; ISBN 978-3-030-12734-3. [Google Scholar]
  16. Padhani, A.R.; Liu, G.; Mu-Koh, D.; Chenevert, T.L.; Thoeny, H.C.; Takahara, T.; Dzik-Jurasz, A.; Ross, B.D.; Van Cauteren, M.; Collins, D.; et al. Diffusion-weighted magnetic resonance imaging as a cancer biomarker: Consensus and recommendations. Neoplasia 2009, 11, 102–125. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  17. Arlinghaus, L.R.; Yankeelov, T.E. Diffusion-weighted MRI. In Quantitative MRI in Cancer; Imaging in Medical Diagnosis and Therapy; CRC Press: Boca Raton, FL, USA, 2011; ISBN 978-1-4398-2058-2. [Google Scholar]
  18. Galbán, C.J.; Ma, B.; Malyarenko, D.; Pickles, M.D.; Heist, K.; Henry, N.L.; Schott, A.F.; Neal, C.H.; Hylton, N.M.; Rehemtulla, A.; et al. Multi-site clinical evaluation of DW-MRI as a treatment response metric for breast cancer patients undergoing neoadjuvant chemotherapy. PLoS ONE 2015, 10. [Google Scholar] [CrossRef] [Green Version]
  19. Yankeelov, T.; Gore, J. Dynamic contrast enhanced magnetic resonance imaging in oncology: Theory, data acquisition, analysis, and examples. Curr. Med. Imaging Rev. 2007, 3, 91–107. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  20. Li, X.; Arlinghaus, L.R.; Ayers, G.D.; Chakravarthy, A.B.; Abramson, R.G.; Abramson, V.G.; Atuegwu, N.; Farley, J.; Mayer, I.A.; Kelley, M.C.; et al. DCE-MRI analysis methods for predicting the response of breast cancer to neoadjuvant chemotherapy: Pilot study findings: DCE-MRI to predict breast cancer treatment response. Magn. Reson. Med. 2014, 71, 1592–1602. [Google Scholar] [CrossRef] [Green Version]
  21. Tudorica, A.; Oh, K.Y.; Chui, S.Y.-C.; Roy, N.; Troxell, M.L.; Naik, A.; Kemmer, K.A.; Chen, Y.; Holtorf, M.L.; Afzal, A.; et al. Early prediction and evaluation of breast cancer response to neoadjuvant chemotherapy using quantitative DCE-MRI. Transl. Oncol. 2016, 9, 8–17. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  22. O’Connor, J.P.B.; Rose, C.J.; Waterton, J.C.; Carano, R.A.D.; Parker, G.J.M.; Jackson, A. Imaging intratumor heterogeneity: Role in therapy response, resistance, and clinical outcome. Clin. Cancer Res. 2015, 21, 249–257. [Google Scholar] [CrossRef] [Green Version]
  23. Just, N. Improving tumour heterogeneity MRI assessment with histograms. Br. J. Cancer 2014, 111, 2205–2213. [Google Scholar] [CrossRef] [Green Version]
  24. O’Connor, E.; Fieller, N.; Holmes, A.; Waterton, J. How to analyse dynamic MRI in oncology; advanced histogram analysis gives better statistical power & insight than simple averaging. Proc. Int. Soc. Magn. Reason. Med. 2005, 13, 2101. [Google Scholar]
  25. Aliu, S.O.; Wilmes, L.J.; Moasser, M.M.; Hann, B.C.; Li, K.-L.; Wang, D.; Hylton, N.M. MRI methods for evaluating the effects of tyrosine kinase inhibitor administration used to enhance chemotherapy efficiency in a breast tumor xenograft model. J. Magn. Reson. Imaging 2009, 29, 1071–1079. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  26. Syed, A.K.; Woodall, R.; Whisenant, J.G.; Yankeelov, T.E.; Sorace, A.G. Characterizing trastuzumab-induced alterations in intratumoral heterogeneity with quantitative imaging and immunohistochemistry in HER2+ breast cancer. Neoplasia 2019, 21, 17–29. [Google Scholar] [CrossRef] [PubMed]
  27. Teruel, J.R.; Heldahl, M.G.; Goa, P.E.; Pickles, M.; Lundgren, S.; Bathen, T.F.; Gibbs, P. Dynamic contrast-enhanced MRI texture analysis for pretreatment prediction of clinical and pathological response to neoadjuvant chemotherapy in patients with locally advanced breast cancer. NMR Biomed. 2014, 27, 887–896. [Google Scholar] [CrossRef]
  28. Chen, W.; Giger, M.L.; Li, H.; Bick, U.; Newstead, G.M. Volumetric texture analysis of breast lesions on contrast-enhanced magnetic resonance images. Magn. Reson. Med. 2007, 58, 562–571. [Google Scholar] [CrossRef]
  29. Chaudhury, B.; Zhou, M.; Goldgof, D.B.; Hall, L.O.; Gatenby, R.A.; Gillies, R.J.; Patel, B.K.; Weinfurtner, R.J.; Drukteinis, J.S. Heterogeneity in intratumoral regions with rapid gadolinium washout correlates with estrogen receptor status and nodal metastasis: ER status and nodal metastasis. J. Magn. Reson. Imaging 2015, 42, 1421–1430. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  30. Méndez, C.A.; Pizzorni Ferrarese, F.; Summers, P.; Petralia, G.; Menegaz, G. DCE-MRI and DWI integration for breast lesions assessment and heterogeneity quantification. Int. J. Biomed. Imaging 2012, 2012, 1–8. [Google Scholar] [CrossRef]
  31. O’Connor, J.P.B. Cancer heterogeneity and imaging. Semin. Cell Dev. Biol. 2017, 64, 48–57. [Google Scholar] [CrossRef]
  32. Wu, J.; Gong, G.; Cui, Y.; Li, R. Intratumor partitioning and texture analysis of dynamic contrast-enhanced (DCE)-MRI identifies relevant tumor subregions to predict pathological response of breast cancer to neoadjuvant chemotherapy. J. Magn. Reson. Imaging 2016. [Google Scholar] [CrossRef]
  33. Featherstone, A.K.; O’Connor, J.P.B.; Little, R.A.; Watson, Y.; Cheung, S.; Babur, M.; Williams, K.J.; Matthews, J.C.; Parker, G.J.M. Data-driven mapping of hypoxia-related tumor heterogeneity using DCE-MRI and OE-MRI: Data-driven identification of hypoxia-related tumor regions. Magn. Reson. Med. 2017. [Google Scholar] [CrossRef]
  34. Chang, Y.-C.C.; Ackerstaff, E.; Tschudi, Y.; Jimenez, B.; Foltz, W.; Fisher, C.; Lilge, L.; Cho, H.; Carlin, S.; Gillies, R.J.; et al. Delineation of tumor habitats based on dynamic contrast enhanced MRI. Sci. Rep. 2017, 7. [Google Scholar] [CrossRef] [Green Version]
  35. Zhou, M.; Hall, L.; Goldgof, D.; Russo, R.; Balagurunathan, Y.; Gillies, R.; Gatenby, R. Radiologically defined ecological dynamics and clinical outcomes in glioblastoma multiforme: Preliminary results. Transl. Oncol. 2014, 7, 5–13. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  36. Jardim-Perassi, B.V.; Martinez, G.; Gillies, R. Habitat imaging of tumor evolution by magnetic resonance imaging (MRI). In Radiomics and Radiogenomics; Li, R., Xing, L., Napel, S., Rubin, D.L., Eds.; Chapman and Hall/CRCCRC Press, Taylor & Francis Group: Boca Raton, FL, USA, 2019; pp. 115–120. ISBN 978-1-351-20827-7. [Google Scholar]
  37. Katiyar, P.; Divine, M.R.; Kohlhofer, U.; Quintanilla-Martinez, L.; Schölkopf, B.; Pichler, B.J.; Disselhorst, J.A. A novel unsupervised segmentation approach quantifies tumor tissue populations using multiparametric MRI: First results with histological validation. Mol. Imaging Biol. 2017, 19, 391–397. [Google Scholar] [CrossRef] [Green Version]
  38. Chaudhury, B.; Zhou, M.; Farhidzadeh, H.; Goldgof, D.B.; Hall, L.O.; Gatenby, R.A.; Gillies, R.J.; Weinfurtner, R.J.; Drukteinis, J.S. Predicting Ki67% expression from DCE-MR images of breast tumors using textural kinetic features in tumor habitats. Proc. SPIE 2016, 9785, 97850T. [Google Scholar]
  39. Jardim-Perassi, B.V.; Huang, S.; Dominguez-Viqueira, W.; Poleszczuk, J.; Budzevich, M.M.; Abdalah, M.A.; Pillai, S.R.; Ruiz, E.; Bui, M.M.; Zuccari, D.A.P.C.; et al. Multiparametric MRI and coregistered histology identify tumor habitats in breast cancer mouse models. Cancer Res. 2019, 79, 3952–3964. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  40. Kakkad, S.; Zhang, J.; Akhbardeh, A.; Jacob, D.; Krishnamachary, B.; Solaiyappan, M.; Jacobs, M.A.; Raman, V.; Leibfritz, D.; Glunde, K.; et al. Collagen fibers mediate MRI-detected water diffusion and anisotropy in breast cancers. Neoplasia 2016, 18, 585–593. [Google Scholar] [CrossRef] [Green Version]
  41. Dean-Colomb, W.; Esteva, F.J. Her2-positive breast cancer: Herceptin and beyond. Eur. J. Cancer 2008, 44, 2806–2812. [Google Scholar] [CrossRef]
  42. Youk, J.H.; Son, E.J.; Chung, J.; Kim, J.-A.; Kim, E. Triple-negative invasive breast cancer on dynamic contrast-enhanced and diffusion-weighted MR imaging: Comparison with other breast cancer subtypes. Eur. Radiol. 2012, 22, 1724–1734. [Google Scholar] [CrossRef]
  43. Agner, S.C.; Rosen, M.A.; Englander, S.; Tomaszewski, J.E.; Feldman, M.D.; Zhang, P.; Mies, C.; Schnall, M.D.; Madabhushi, A. Computerized image analysis for identifying triple-negative breast cancers and differentiating them from other molecular subtypes of breast cancer on dynamic contrast-enhanced MR images: A feasibility study. Radiology 2014, 272, 91–99. [Google Scholar] [CrossRef] [Green Version]
  44. Li, H.; Zhu, Y.; Burnside, E.S.; Huang, E.; Drukker, K.; Hoadley, K.A.; Fan, C.; Conzen, S.D.; Zuley, M.; Net, J.M.; et al. Quantitative MRI radiomics in the prediction of molecular classifications of breast cancer subtypes in the TCGA/TCIA data set. NPJ Breast Cancer 2016, 2, 16012. [Google Scholar] [CrossRef]
  45. Fan, M.; Li, H.; Wang, S.; Zheng, B.; Zhang, J.; Li, L. Radiomic analysis reveals DCE-MRI features for prediction of molecular subtypes of breast cancer. PLoS ONE 2017, 12, e0171683. [Google Scholar] [CrossRef]
  46. Sorace, A.G.; Syed, A.K.; Barnes, S.L.; Quarles, C.C.; Sanchez, V.; Kang, H.; Yankeelov, T.E. Quantitative [18F]FMISO PET imaging shows reduction of hypoxia following trastuzumab in a murine model of HER2+ breast cancer. Mol. Imaging Biol. 2017, 19, 130–137. [Google Scholar] [CrossRef] [Green Version]
  47. Izumi, Y.; Xu, L.; Di Tomaso, E.; Fukumura, D.; Jain, R.K. Jain Herceptin acts as an anti-angiogenic cocktail. Nature 2002, 416, 279–280. [Google Scholar] [CrossRef] [PubMed]
  48. McCormack, D.R.; Walsh, A.J.; Sit, W.; Arteaga, C.L.; Chen, J.; Cook, R.S.; Skala, M.C. In vivo hyperspectral imaging of microvessel response to trastuzumab treatment in breast cancer xenografts. Biomed. Opt. Express 2014, 5, 2247–2261. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  49. Hardee, M.E.; Eapen, R.J.; Rabbani, Z.N.; Dreher, M.R.; Marks, J.; Blackwell, K.L.; Dewhirst, M.W. Her2/neu signaling blockade improves tumor oxygenation in a multifactorial fashion in Her2/neu+ tumors. Cancer Chemother. Pharmacol. 2009, 63, 219–228. [Google Scholar] [CrossRef] [Green Version]
  50. Sorace, A.G.; Quarles, C.C.; Whisenant, J.G.; Hanker, A.B.; McIntyre, J.O.; Sanchez, V.M.; Yankeelov, T.E. Trastuzumab improves tumor perfusion and vascular delivery of cytotoxic therapy in a murine model of HER2+ breast cancer: Preliminary results. Breast Cancer Res. Treat. 2016, 155, 273–284. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  51. Whisenant, J.G.; Sorace, A.G.; McIntyre, J.O.; Hakmook, K.; Sánchez, V.; Loveless, M.E.; Yankeelov, T.E. Evaluating treatment response using DW-MRI and DCE-MRI in trastuzumab responsive and resistant HER2-overexpressing human breast cancer xenografts. Transl. Oncol. 2014, 7, 768–779. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  52. Yardley, D.A. nab-Paclitaxel mechanisms of action and delivery. J. Controlled Release 2013, 170, 365–372. [Google Scholar] [CrossRef]
  53. Bocci, G.; Di Paolo, A.; Danesi, R. The pharmacological bases of the antiangiogenic activity of paclitaxel. Angiogenesis 2013, 16, 481–492. [Google Scholar] [CrossRef] [Green Version]
  54. Belotti, D.; Vergani, V.; Drudis, T.; Borsotti, P.; Pitelli, M.R.; Viale, G.; Giavazzi, R.; Taraboletti, G. The microtubule-affecting drug paclitaxel has antiangiogenic activity. Clin. Cancer Res. 1996, 2, 1843–1849. [Google Scholar] [PubMed]
  55. Sun, X.; Yang, L.; Yan, X.; Sun, Y.; Zhao, D.; Ji, Y.; Wang, K.; Chen, X.; Shen, B. DCE-MRI-derived parameters in evaluating abraxane-induced early vascular response and the effectiveness of its synergistic interaction with cisplatin. PLoS ONE 2016, 11, e0162601. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  56. Egeland, T.A.M.; Gulliksrud, K.; Gaustad, J.-V.; Mathiesen, B.; Rofstad, E.K. Dynamic contrast-enhanced-MRI of tumor hypoxia. Magn. Reson. Med. 2012, 67, 519–530. [Google Scholar] [CrossRef]
  57. Stoyanova, R.; Huang, K.; Sandler, K.; Cho, H.; Carlin, S.; Zanzonico, P.B.; Koutcher, J.A.; Ackerstaff, E. Mapping tumor hypoxia in vivo using pattern recognition of dynamic contrast-enhanced MRI data. Transl. Oncol. 2012, 5, 437. [Google Scholar] [CrossRef] [PubMed]
  58. Berry, L.R.; Barck, K.H.; Go, M.A.; Ross, J.; Wu, X.; Williams, S.P.; Gogineni, A.; Cole, M.J.; Van Bruggen, N.; Fuh, G.; et al. Quantification of viable tumor microvascular characteristics by multispectral analysis. Magn. Reson. Med. 2008, 60, 64–72. [Google Scholar] [CrossRef] [PubMed]
  59. Henning, E.C.; Azuma, C.; Sotak, C.H.; Helmer, K.G. Multispectral tissue characterization in a RIF-1 tumor model: Monitoring the ADC and T2 responses to single-dose radiotherapy. Part II. Magn. Reson. Med. 2007, 57, 513–519. [Google Scholar] [CrossRef] [PubMed]
  60. Longo, D.L.; Dastrù, W.; Consolino, L.; Espak, M.; Arigoni, M.; Cavallo, F.; Aime, S. Cluster analysis of quantitative parametric maps from DCE-MRI: Application in evaluating heterogeneity of tumor response to antiangiogenic treatment. Magn. Reson. Imaging 2015, 33, 725–736. [Google Scholar] [CrossRef] [PubMed]
  61. Alic, L.; Haeck, J.C.; Bol, K.; Klein, S.; Van Tiel, S.T.; Wielepolski, P.A.; de Jong, M.; Niessen, W.J.; Bernsen, M.; Veenland, J.F. Facilitating tumor functional assessment by spatially relating 3D tumor histology and in vivo MRI: Image registration approach. PLoS ONE 2011, 6, e22835. [Google Scholar] [CrossRef] [Green Version]
  62. Sorace, A.G.; Wu, C.; Barnes, S.L.; Jarrett, A.M.; Avery, S.; Patt, D.; Goodgame, B.; Luci, J.J.; Kang, H.; Abramson, R.G.; et al. Repeatability, reproducibility, and accuracy of quantitative MRI of the breast in the community radiology setting. J. Magn. Reson. Imaging 2018, 48, 695–707. [Google Scholar] [CrossRef]
  63. Galbraith, S.M.; Lodge, M.A.; Taylor, N.J.; Rustin, G.J.S.; Bentzen, S.; Stirling, J.J.; Padhani, A.R. Reproducibility of dynamic contrast-enhanced MRI in human muscle and tumours: Comparison of quantitative and semi-quantitative analysis. NMR Biomed. 2002, 15, 132–142. [Google Scholar] [CrossRef]
  64. Spick, C.; Bickel, H.; Pinker, K.; Bernathova, M.; Kapetas, P.; Woitek, R.; Clauser, P.; Polanec, S.H.; Rudas, M.; Bartsch, R.; et al. Diffusion-weighted MRI of breast lesions: A prospective clinical investigation of the quantitative imaging biomarker characteristics of reproducibility, repeatability, and diagnostic accuracy. NMR Biomed. 2016, 29, 1445–1453. [Google Scholar] [CrossRef]
  65. Partridge, S.C.; Amornsiripanitch, N. The role of DWI in the assessment of breast lesions. Top. Magn. Reson. Imaging TMRI 2017, 26, 201–209. [Google Scholar] [CrossRef]
  66. Stringfield, O.; Arrington, J.A.; Johnston, S.K.; Rognin, N.G.; Peeri, N.C.; Balagurunathan, Y.; Jackson, P.R.; Clark-Swanson, K.R.; Swanson, K.R.; Egan, K.M.; et al. Multiparameter MRI predictors of long-term survival in glioblastoma multiforme. Tomography 2019, 5, 135–144. [Google Scholar] [CrossRef] [PubMed]
  67. Lee, J.; Narang, S.; Martinez, J.; Rao, G.; Rao, A. Spatial habitat features derived from multiparametric magnetic resonance imaging data are associated with molecular subtype and 12-month survival status in glioblastoma multiforme. PLoS ONE 2015, 10. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  68. Parra, N.A.; Lu, H.; Choi, J.; Gage, K.; Pow-Sang, J.; Gillies, R.J.; Balagurunathan, Y. Habitats in DCE-MRI to predict clinically significant prostate cancers. Tomography 2019, 5, 68–76. [Google Scholar] [CrossRef] [PubMed]
  69. Stoyanova, R.; Pollack, A.; Takhar, M.; Lynne, C.; Parra, N.; Lam, L.L.; Alshalalfa, M.; Buerki, C.; Castillo, R.; Jorda, M.; et al. Association of multiparametric MRI quantitative imaging features with prostate cancer gene expression in MRI-targeted prostate biopsies. Oncotarget 2016, 7, 53362. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  70. Barnes, S.L.; Sorace, A.G.; Whisenant, J.G.; McIntyre, J.O.; Kang, H.; Yankeelov, T.E. DCE- and DW-MRI as early imaging biomarkers of treatment response in a preclinical model of triple negative breast cancer. NMR Biomed. 2017, 30, e3799. [Google Scholar] [CrossRef]
  71. Wilson, G.H.; Gore, J.C.; Yankeelov, T.E.; Barnes, S.; Peterson, T.E.; True, J.M.; Shokouhi, S.; McIntyre, J.O.; Sanders, M.; Abramson, V.; et al. An approach to breast cancer diagnosis via PET imaging of microcalcifications using 18F-NaF. J. Nucl. Med. 2014, 55, 1138–1143. [Google Scholar] [CrossRef] [Green Version]
  72. Kety, S.S. The theory and applications of the exchange of inert gas at the lungs and tissues. Pharmacol. Rev. 1951, 3, 1–41. [Google Scholar]
  73. Tofts, P.S.; Brix, G.; Buckley, D.L.; Evelhoch, J.L.; Henderson, E.; Knopp, M.V.; Larsson, H.B.W.; Lee, T.-Y.; Mayr, N.A.; Parker, G.J.M.; et al. Estimating kinetic parameters from dynamic contrast-enhanced t1-weighted MRI of a diffusable tracer: Standardized quantities and symbols. J. Magn. Reson. Imaging 1999, 10, 223–232. [Google Scholar] [CrossRef]
  74. Loveless, M.E.; Halliday, J.; Liess, C.; Xu, L.; Dortch, R.D.; Whisenant, J.; Waterton, J.C.; Gore, J.C.; Yankeelov, T.E. A quantitative comparison of the influence of individual versus population-derived vascular input functions on dynamic contrast enhanced-MRI in small animals. Magn. Reson. Med. 2012, 67, 226–236. [Google Scholar] [CrossRef] [Green Version]
  75. Loveless, M.E.; Lawson, D.; Collins, M.; Nadella, M.V.P.; Reimer, C.; Huszar, D.; Halliday, J.; Waterton, J.C.; Gore, J.C.; Yankeelov, T.E. Comparisons of the efficacy of a Jak1/2 inhibitor (AZD1480) with a VEGF signaling inhibitor (cediranib) and sham treatments in mouse tumors using DCE-MRI, DW-MRI, and histology. Neoplasia 2012, 14, 54–64. [Google Scholar] [CrossRef]
  76. Li, X.; Rooney, W.D.; Várallyay, C.G.; Gahramanov, S.; Muldoon, L.L.; Goodman, J.A.; Tagge, I.J.; Selzer, A.H.; Pike, M.M.; Neuwelt, E.A.; et al. Dynamic contrast enhanced-MRI with extravasating contrast reagent: Rat cerebral glioma blood volume determination. J. Magn. Reson. 2010, 206, 190–199. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  77. Wu, C.; Pineda, F.; Hormuth, D.A.; Karczmar, G.S.; Yankeelov, T.E. Quantitative analysis of vascular properties derived from ultrafast DCE-MRI to discriminate malignant and benign breast tumors. Magn. Reson. Med. 2019, 81, 2147–2160. [Google Scholar] [CrossRef] [PubMed]
  78. Han, J.; Kamber, M.; Pei, J. Data preprocessing. In Data Mining: Concepts and Techniques; Elsevier Science & Technology: St. Louis, MO, USA, 2011; ISBN 978-0-12-381480-7. [Google Scholar]
  79. Nielsen, F. Hierarchical clustering. In Introduction to HPC with MPI for Data Science; Undergraduate Topics in Computer Science; Springer International Publishing: Cham, Switzerland, 2016; pp. 195–211. ISBN 978-3-319-21902-8. [Google Scholar]
  80. Tibshirani, R.; Walther, G.; Hastie, T. Estimating the number of clusters in a data set via the gap statistic. J. R. Stat. Soc. Ser. B Stat. Methodol. 2001, 63, 411–423. [Google Scholar] [CrossRef]
  81. Wu, J.; Cao, G.; Sun, X.; Lee, J.; Rubin, D.L.; Napel, S.; Kurian, A.W.; Daniel, B.L.; Li, R. Intratumoral spatial heterogeneity at perfusion MR imaging predicts recurrence-free survival in locally advanced breast cancer treated with neoadjuvant chemotherapy. Radiology 2018, 172462. [Google Scholar] [CrossRef]
  82. Linkert, M.; Rueden, C.T.; Allan, C.; Burel, J.-M.; Moore, W.; Patterson, A.; Loranger, B.; Moore, J.; Neves, C.; MacDonald, D.; et al. Metadata matters: Access to image data in the real world. J. Cell Biol. 2010, 189, 777–782. [Google Scholar] [CrossRef] [Green Version]
  83. Otsu, N. A threshold selection method from gray-level histograms. IEEE Trans. Syst. Man Cybern. 1979, 9, 62–66. [Google Scholar] [CrossRef] [Green Version]
  84. Maley, C.C.; Koelble, K.; Natrajan, R.; Aktipis, A.; Yuan, Y. An ecological measure of immune-cancer colocalization as a prognostic factor for breast cancer. Breast Cancer Res. 2015, 17. [Google Scholar] [CrossRef] [Green Version]
Cancers 12 01682 g001 550
Figure 1. Pipeline for the habitat imaging analysis. Longitudinal DCE- and DW-MRI data was acquired for each animal and processed to generate individual multi-parametric maps of the tumor region (top row, left). Each tumor voxel was then defined as a vector where each element corresponded to one of the parameters returned from the quantitative MRI analysis. Voxel data within each individual tumor cohort (BT-474 or MDA-MB-231) were pooled together. The pooled voxel data is then clustered to identify tumor habitats (middle row, left). Tumor composition over time was quantified in terms of the percent of each tumor’s volume having a particular habitat (bottom row, left). End point histology data was acquired and processed to generate stain density maps (top row, right). Each pixel was then defined as a vector of its stain density values and pixel data within each tumor cohort were pooled together. The pooled pixel data is clustered to identify histological tumor habitats (middle row, right). Histological tumor composition was then quantified in terms of the percent of each tumor’s area having a particular habitat (bottom row, right). To biologically validate imaging-derived habitats, a correlation analysis was performed between habitat percent tumor volume from the imaging-derived habitats, and habitat percent tumor area from the histology-derived habitats.
Figure 1. Pipeline for the habitat imaging analysis. Longitudinal DCE- and DW-MRI data was acquired for each animal and processed to generate individual multi-parametric maps of the tumor region (top row, left). Each tumor voxel was then defined as a vector where each element corresponded to one of the parameters returned from the quantitative MRI analysis. Voxel data within each individual tumor cohort (BT-474 or MDA-MB-231) were pooled together. The pooled voxel data is then clustered to identify tumor habitats (middle row, left). Tumor composition over time was quantified in terms of the percent of each tumor’s volume having a particular habitat (bottom row, left). End point histology data was acquired and processed to generate stain density maps (top row, right). Each pixel was then defined as a vector of its stain density values and pixel data within each tumor cohort were pooled together. The pooled pixel data is clustered to identify histological tumor habitats (middle row, right). Histological tumor composition was then quantified in terms of the percent of each tumor’s area having a particular habitat (bottom row, right). To biologically validate imaging-derived habitats, a correlation analysis was performed between habitat percent tumor volume from the imaging-derived habitats, and habitat percent tumor area from the histology-derived habitats.
Cancers 12 01682 g001
Cancers 12 01682 g002 550
Figure 2. The mean parameter values for each of the identified habitats are shown for the (a) BT-474 and (b) MDA-MB-231 cohorts. The mean parameter values were used to identify the physiological habitat each cluster represented in terms of high or low “vascularity” or “cellularity”. The three identified habitats were as follows: high-vascularity high-cellularity (HV-HC), low-vascularity high-cellularity (LV-HC), low-vascularity low-cellularity (LV-LC). All three habitats were found in both the BT-474 and MDA-MB-231 cohorts. Error bars show standard deviation and ** indicates a p < 0.01. Panel (c) shows representative parameter maps and a corresponding habitat map for an example BT-474 tumor at baseline, where the HV-HC habitat is shown in red, LV-HC in green, and LV-LC in blue.
Figure 2. The mean parameter values for each of the identified habitats are shown for the (a) BT-474 and (b) MDA-MB-231 cohorts. The mean parameter values were used to identify the physiological habitat each cluster represented in terms of high or low “vascularity” or “cellularity”. The three identified habitats were as follows: high-vascularity high-cellularity (HV-HC), low-vascularity high-cellularity (LV-HC), low-vascularity low-cellularity (LV-LC). All three habitats were found in both the BT-474 and MDA-MB-231 cohorts. Error bars show standard deviation and ** indicates a p < 0.01. Panel (c) shows representative parameter maps and a corresponding habitat map for an example BT-474 tumor at baseline, where the HV-HC habitat is shown in red, LV-HC in green, and LV-LC in blue.
Cancers 12 01682 g002
Cancers 12 01682 g003 550
Figure 3. Panel (a) shows representative tumor habitat maps for a BT-474 control (CTL) and trastuzumab-treated (TRA) tumor on days 0 (pre-treatment), 1, and 4. Trastuzumab-treated tumors showed longitudinal increases in percent tumor volume (panel (b)) of the HV-HC habitat at day 4 (p = 0.02), and significantly higher percent tumor volumes of the HV-HC habitat compared to control (p = 0.03). Trastuzumab-treated tumors also showed longitudinal decreases in percent tumor volume (panel (b)) of the LV-HC habitat at day 4 (p < 0.01) and significantly lower percent tumor volumes of the LV-HC habitat compared to control (p = 0.03). Panel (c) shows representative tumor habitat maps for an MDA-MB-231 control (CTL) and low-dose albumin-bound paclitaxel (ABP-LD) tumor on days 0 (pre-treatment), 1, and 3. Low-dose ABP-treated tumors showed significant decreases in percent tumor volume (panel (d)) of the LV-HC habitat at day 1 (p < 0.01) and day 3 (p = 0.01). No significant changes were observed for control and high-dose ABP (ABP-HD) tumors over the course of the study. No significant differences were observed between treatment groups for the MDA-MB-231 cohort. Statistical significance is denoted with a * for p < 0.05 and ** for p < 0.01, error bars represent standard error.
Figure 3. Panel (a) shows representative tumor habitat maps for a BT-474 control (CTL) and trastuzumab-treated (TRA) tumor on days 0 (pre-treatment), 1, and 4. Trastuzumab-treated tumors showed longitudinal increases in percent tumor volume (panel (b)) of the HV-HC habitat at day 4 (p = 0.02), and significantly higher percent tumor volumes of the HV-HC habitat compared to control (p = 0.03). Trastuzumab-treated tumors also showed longitudinal decreases in percent tumor volume (panel (b)) of the LV-HC habitat at day 4 (p < 0.01) and significantly lower percent tumor volumes of the LV-HC habitat compared to control (p = 0.03). Panel (c) shows representative tumor habitat maps for an MDA-MB-231 control (CTL) and low-dose albumin-bound paclitaxel (ABP-LD) tumor on days 0 (pre-treatment), 1, and 3. Low-dose ABP-treated tumors showed significant decreases in percent tumor volume (panel (d)) of the LV-HC habitat at day 1 (p < 0.01) and day 3 (p = 0.01). No significant changes were observed for control and high-dose ABP (ABP-HD) tumors over the course of the study. No significant differences were observed between treatment groups for the MDA-MB-231 cohort. Statistical significance is denoted with a * for p < 0.05 and ** for p < 0.01, error bars represent standard error.
Cancers 12 01682 g003
Cancers 12 01682 g004 550
Figure 4. The mean stain density values for each of the histological habitats are shown for the (a) BT-474 and (b) MDA-MB-231 cohorts. The mean vessel count and nuclei count values were used to identify the physiological habitat each cluster represented in terms of high or low “vascularity” or “cellularity”, respectively. The three identified habitats were as follows: high-vascularity high-cellularity (HV-HC), low-vascularity high-cellularity (LV-HC), low-vascularity low-cellularity (LV-LC). All three habitats were found in both the BT-474 and MDA-MB-231 cohorts. Error bars show standard deviation and ** indicates a p < 0.01. Panel (c) shows representative stain density maps and a corresponding habitat map for an example MDA-MB-231 tumor from the low-dose ABP treatment group, where the HV-HC habitat is shown in red, LV-HC in green, and LV-LC in blue.
Figure 4. The mean stain density values for each of the histological habitats are shown for the (a) BT-474 and (b) MDA-MB-231 cohorts. The mean vessel count and nuclei count values were used to identify the physiological habitat each cluster represented in terms of high or low “vascularity” or “cellularity”, respectively. The three identified habitats were as follows: high-vascularity high-cellularity (HV-HC), low-vascularity high-cellularity (LV-HC), low-vascularity low-cellularity (LV-LC). All three habitats were found in both the BT-474 and MDA-MB-231 cohorts. Error bars show standard deviation and ** indicates a p < 0.01. Panel (c) shows representative stain density maps and a corresponding habitat map for an example MDA-MB-231 tumor from the low-dose ABP treatment group, where the HV-HC habitat is shown in red, LV-HC in green, and LV-LC in blue.
Cancers 12 01682 g004
Cancers 12 01682 g005 550
Figure 5. Correlation of MRI-derived tumor habitats with histological habitats. Stain density maps were used to generate histological habitat maps. The percent tumor area of each habitat was calculated for each histological section (panel (a)). Central imaging slices were selected from study endpoint (day 4 for BT-474 cohort, day 3 for MDA-MB-231 cohort) MRI data (panel (a)) and percent tumor volume of each habitat was calculated for each tumor’s central slice. The correlation was calculated between percent tumor volume (MRI) and percent tumor area (histology) for each habitat. For the BT-474 cohort (panel (b), left column), percent tumor volume of MRI-derived habitats showed a significant positive correlation with percent tumor area of histology-derived habitats for the HV-HC (r2 = 0.74, p = 0.03) and LV-LC (r2 = 0.70, p = 0.04) habitats. For the MDA-MB-231 cohort (panel (b), right column), percent tumor volume of MRI-derived habitats showed a significant positive correlation with percent tumor area of histology-derived habitats for the LV-HC (r2 = 0.70, p = 0.04) and LV-LC habitats (r2 = 0.73, p < 0.01).
Figure 5. Correlation of MRI-derived tumor habitats with histological habitats. Stain density maps were used to generate histological habitat maps. The percent tumor area of each habitat was calculated for each histological section (panel (a)). Central imaging slices were selected from study endpoint (day 4 for BT-474 cohort, day 3 for MDA-MB-231 cohort) MRI data (panel (a)) and percent tumor volume of each habitat was calculated for each tumor’s central slice. The correlation was calculated between percent tumor volume (MRI) and percent tumor area (histology) for each habitat. For the BT-474 cohort (panel (b), left column), percent tumor volume of MRI-derived habitats showed a significant positive correlation with percent tumor area of histology-derived habitats for the HV-HC (r2 = 0.74, p = 0.03) and LV-LC (r2 = 0.70, p = 0.04) habitats. For the MDA-MB-231 cohort (panel (b), right column), percent tumor volume of MRI-derived habitats showed a significant positive correlation with percent tumor area of histology-derived habitats for the LV-HC (r2 = 0.70, p = 0.04) and LV-LC habitats (r2 = 0.73, p < 0.01).
Cancers 12 01682 g005
Cancers 12 01682 g006 550
Figure 6. Tumor implantation, imaging, and treatment schedule for the tumor cohorts. MDA-MB-231 tumors were imaged on days 0, 1, and 3, while BT-474 tumors were imaged on days 0, 1, 4. MDA-MB-231 tumors were treated on day 0 (post-imaging) and day 2, with either saline (control), 15 mg/kg albumin-bound paclitaxel (ABP), or 25 mg/kg ABP treatment. BT-474 tumors were treated on days 0 (post-imaging) and 3, with either saline (control) of 10 mg/kg trastuzumab treatment.
Figure 6. Tumor implantation, imaging, and treatment schedule for the tumor cohorts. MDA-MB-231 tumors were imaged on days 0, 1, and 3, while BT-474 tumors were imaged on days 0, 1, 4. MDA-MB-231 tumors were treated on day 0 (post-imaging) and day 2, with either saline (control), 15 mg/kg albumin-bound paclitaxel (ABP), or 25 mg/kg ABP treatment. BT-474 tumors were treated on days 0 (post-imaging) and 3, with either saline (control) of 10 mg/kg trastuzumab treatment.
Cancers 12 01682 g006

Share and Cite

MDPI and ACS Style

Syed, A.K.; Whisenant, J.G.; Barnes, S.L.; Sorace, A.G.; Yankeelov, T.E. Multiparametric Analysis of Longitudinal Quantitative MRI Data to Identify Distinct Tumor Habitats in Preclinical Models of Breast Cancer. Cancers 2020, 12, 1682. https://doi.org/10.3390/cancers12061682

AMA Style

Syed AK, Whisenant JG, Barnes SL, Sorace AG, Yankeelov TE. Multiparametric Analysis of Longitudinal Quantitative MRI Data to Identify Distinct Tumor Habitats in Preclinical Models of Breast Cancer. Cancers. 2020; 12(6):1682. https://doi.org/10.3390/cancers12061682

Chicago/Turabian Style

Syed, Anum K., Jennifer G. Whisenant, Stephanie L. Barnes, Anna G. Sorace, and Thomas E. Yankeelov. 2020. "Multiparametric Analysis of Longitudinal Quantitative MRI Data to Identify Distinct Tumor Habitats in Preclinical Models of Breast Cancer" Cancers 12, no. 6: 1682. https://doi.org/10.3390/cancers12061682

APA Style

Syed, A. K., Whisenant, J. G., Barnes, S. L., Sorace, A. G., & Yankeelov, T. E. (2020). Multiparametric Analysis of Longitudinal Quantitative MRI Data to Identify Distinct Tumor Habitats in Preclinical Models of Breast Cancer. Cancers, 12(6), 1682. https://doi.org/10.3390/cancers12061682

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop