Abstract

Objective: To determine if tumor habitat perfusion and diffusion characteristics are related to treatment response in soft tissue sarcoma patients.

Methods: Eight patients (58.6 ± 10.2 y/o) with soft tissue sarcomas underwent pre-treatment Dynamic Contrast Enhanced- (DCE-) and Diffusion Weighted-MRI and post-treatment resection allowing pathology to determine treatment response. Tumors were manually segmented from T1-weighted post-contrast images. Tumor habitats, classified as well-perfused, hypoxic and necrotic, were determined from DCE-MRI using a pattern recognition technique. Each habitat was characterized with pharmacokinetic parameters, k_ep and K^trans, calculated from the DCE-MRI sequence and the apparent diffusion coefficients from the co-registered Diffusion Weighted-MRI. The quantitative imaging features were examined for associations with treatment response. Patients were classified as responders and nonresponders based on histopathology.

Results: Using unsupervised clustering on the imaging features, the patients were divided accurately between responders and nonresponders. Out of all determined features, k_ep (p = 0.04), K^trans (p < 0.01), and percent volume (p = 0.02) for well-perfused habitats were significantly lower in nonresponders, whereas volume (p = 0.04) and percent volume (p < 0.01) for necrotic habitats were significantly higher in responders.

Conclusion: Prediction of treatment response in soft tissue sarcoma patients yielded promising results when utilizing differences in Dynamic Contrast Enhanced- and Diffusion Weighted-MRI features between unique tumor habitats.

Keywords

Pattern Recognition; DCE-MRI; DW-MRI; mpMRI; Soft Tissue Sarcoma

Introduction

Soft tissue sarcomas are a rare but diverse form of cancer with over 50 different subtypes [1]. Due to the heterogeneity of sarcomas, patients often require personalized treatment plans involving radiation, chemotherapy, and surgery, but treatment often is ineffective with a five year survival rate of 64% [2]. Methods to evaluate and monitor soft tissue sarcomas are needed to improve management of the disease. Among potential prognostic factors, the amount of necrosis in soft tissue sarcoma tumors has been found to be predictive of local outcome and recurrence based on resected tumor samples [3]. An in vivo technique to identify tumor composition would potentially allow prediction of patient response to treatment. In imaging, the concept of “habitats” was introduced to map tumor heterogeneity and quantify amounts of each different microenvironment [4–6].

Several studies have explored the use of Dynamic Contrast Enhanced-MRI (DCE-MRI) for prediction of treatment response in soft tissue sarcoma patients [7–10]. DCE-MRI measures perfusion by observing the flow of contrast through tissue over time, producing signal-versus-time curves. Pharmacokinetic fitting of signal-versus-time curves in tumor voxels revealed parameters that allowed prediction of treatment response with promising results [7, 9, 10]. Another functional technique, Diffusion Weighted Imaging (DWI), may provide additional predictive parameters. DWI offers a measure of cellularity by observing water diffusion patterns in tissue, producing Apparent Diffusion Coefficient (ADC) values. DWI has never been demonstrated on pre-treatment prediction of response in sarcoma patients but has shown promise in treatment monitoring [11–15]. However, analyses of both DCE-MRI and DWI datasets in these studies are still performed on the whole tumor ROI, considering the tumor as a homogenous region. The presence of heterogeneous habitats could potentially be masked by an analysis averaging over the entire area of the tumor.

Consideration of signal-versus-time curves and ADC values from different tumor habitats would offer additional features for prediction while providing a more accurate depiction of tumor state than grouping curves from different habitats in the same analyses. Using pattern recognition techniques, a method to determine the location and distribution of tumor habitats from DCE-MRI datasets has been previously demonstrated [5]. Well-perfused, hypoxic, and necrotic habitats were determined in pre-clinical prostate and brain tumor models, as well as in clinical sarcoma and prostate datasets. Habitat quantitative features, extracted from DCE-MRI and DWI sequences, could potentially be associated with the treatment outcomes of soft tissue sarcomas through a radiomics approach. Radiomics data have proven capable for building descriptive and predictive models relating image features to outcomes, including in sarcomas [16].

In this paper, we utilize knowledge of the location and volumes of tumor habitats in soft tissue sarcomas to generate habitat-specific features for use in prediction of treatment response. Features include pharmacokinetic parameters from modeling of signal-versus-time curves from DCE-MRI and ADC values. The use of multiple features and features native to individual habitats is explored to provide more accurate prediction that may facilitate effective treatment management for soft tissue sarcoma patients.

Methods

Study Protocol

All procedures performed in studies involving human participants were in accordance with the ethical standards of the institutional and/or national research committee and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards. Eight patients (58.6 ± 10.2 y/o) diagnosed with soft tissue sarcomas were enrolled in an Institutional Review Board-approved study. For this retrospective study, informed consent was waved. DCE- and DWI sequences were acquired for each patient, prior to receiving neoadjuvant chemo/radiation therapy. Pathology was performed upon resection of tumors, revealing percent necrosis of sample by which response was determined and by which results from MRI analysis could be compared. Patients’ clinical information and treatment response are described in Table 1.

Table 1. Patient clinical information and treatment response.

MRI Protocol

Imaging parameters for DW and DCE-MRI scans for each patient are listed in detail in Table 2. Briefly, were carried out either on a 3T Skyra or 1.5T Symphony MR Scanner (Siemens, Erlangen, Germany).

Table 2. Patient imaging parameters.

Abbreviations: TR = Repetition time; TE = Echo time.

DCE-MRI analysis

Tumor Volumes of Interest (VOIs) for each patient were contoured manually in MIM (MIM, Cleveland, Ohio) by an experienced skeletomuscular radiologist (T.S.). DCE-MRI signal-versus-time curves from VOIs were analyzed using Constrained Non-Negative Matrix Factorization (cNMF), as detailed in previous publications [5, 17], to reveal the distribution of habitats. The method uncovers the main patterns underlying the signal-versus-time curves and the weights of these constituent patterns in each voxel. The number of these patterns, associated with different tumor habitats was set to three to account for well-perfused, hypoxic, and necrotic habitats.

To assign voxels to specific habitats, a threshold of 60% was applied to the cNMF weights, whereby voxels containing a given habitat at a greater fraction than the threshold were assigned to that habitat. For instance, if for a given voxel the fraction of hypoxic, necrotic and well-perfused habitats were 10%, 20% and 70%, the voxel was considered part of the well-perfused habitat. Conversely, if the contribution of the habitats were 30%, 30% and 40%, the voxel was not assigned to any habitat. The signal-versus-time curves for the voxels corresponding to a given habitat were averaged together, producing a representative curve for each habitat. Pharmacokinetic modeling was performed on these representative curves using a Tofts model [18, 19]. Using the synthetic Parker fixed population average Arterial Input Function (AIF) [20] the volume transfer constant between plasma and Extracellular Extravascular Space (EES), K^trans (related to perfusion and permeability per unit volume of tissue), the fractional volume v_e of the EES and k_ep, the ratio of K^trans and v_e, also known as the efflux rate constant, were quantified. In addition, the volume of each habitat was determined and also divided against the overall tumor volume to obtain percent volume of habitats.

ADC analysis

DCE-MRI and ADC maps were co-registered in MIM. The same tumor VOIs and tumor habitats determined from the DCE-MRI analysis were applied to the ADC. Average ADC values for voxels corresponding to well-perfused, hypoxic, and necrotic habitats were determined.

Prediction of patient response

The association between (i) the imaging features for each habitat: K^trans, k_ep, volume, percent from tumor volume, and ADC and (ii) patient response to neoadjuvant therapy was investigated using hierarchical clustering of patient features. Features involving habitats that were not present in a given tumor were set to 0. The final two clusters at the top level of the clustergram were determined to be the groups of patients corresponding to patients predicted to be responders and nonresponders.

Results

Quantitative measures, including habitat-specific K^trans, k_ep, volume, percent volume, and mean ADC for well-perfused, hypoxic, and necrotic habitats are listed in Table 3. Three habitats were identified in 7 patients and 2 habitats in one patient. Roughly, similar fractions of necrotic habitats were determined between the present approach and pathology. An example of habitat locations and signal-versus-time curves in a patient is shown in Figure 1. K^trans showed an increasing trend from necrotic to hypoxic to well-perfused habitats, which generally contained voxels with curves with the highest K^trans. Similarly, k_ep showed a decreasing trend from necrotic to hypoxic to well-perfused habitats, which generally contained voxels with curves with the lowest k_ep. Volume, percent volume, and mean ADC did not have an apparent trend with habitat type.

Table 3. Quantitative measures for patients based on analysis of DCE-MRI data.

Abbreviations: P = Well-perfused; H = Hypoxic; N = Necrotic; Vol = Volume; Per = Percent Volume.

Figure 1. (A) Axial DCE-MRI of responding 56 year-old extraskeletal osteosarcoma female patient post-treatment. Well-perfused, hypoxic, and necrotic habitats determined in the patient are overlaid with red, green, and blue respectively. Early enhancing image from DCE-MRI (left) and ADC (right) are shown in inset from left to right. (B) Axial DCE-MRI of nonresponding 53 year-old synovial sarcoma female patient post-treatment. Well-perfused, hypoxic, and necrotic habitats determined in the patient are overlaid with red, green, and blue respectively. Early enhancing image from DCE-MRI (left) and ADC (right) are shown in inset from left to right. (C) Corresponding average signal-versus-time curves for each habitat are shown below for Fig. 1a. (D) Corresponding average signal-versus-time curves for each habitat are shown below for Fig. 1b.

Hierarchical clustering was performed on subject data using the quantitative measures described above (Figure 2). Two clusters were readily apparent, separating responding from non-responding patients. Parameters which were significantly correlated with response after correction for multiple comparison with Tukey’s honest significant difference test were k_ep well-perfused (p = 0.04), K^trans well-perfused (p < 0.01), volume necrotic (p = 0.04), percent well-perfused (p = 0.02), and percent necrotic (p < 0.01). The well-perfused habitat’s characteristics: k_ep, K^trans, and percent volume were lower in responders whereas volume necrotic and percent necrotic were higher in responders.

Figure 2. Heatmap showing hierachical clustering of 8 patients using quantitative imaging features. Nonresponders, patients 5, 7, and 8, are shown to be in a separate cluster from responders, patients 1, 2, 3, 4, and 6.

Discussion

In this feasibility study, a radiomics approach utilizing perfusion, diffusion and volumetric information from unique tumor habitats was applied to the prediction of treatment response in soft tissue sarcoma patients. Patients from a small sample set were correctly classified into responders and nonresponders, and habitat-related features were identified that had differentiating potential for the two populations. The success of the approach suggests that consideration of tumor heterogeneity is important for treatment prediction and that the state of the tumor may not be accurately represented by analyses of the whole tumor.

This importance of analysis of unique habitats to prognosis has primarily been recognized in breast cancer tumors. DCE-MRI studies of breast cancer have shown tumors with a larger proportion of voxels with fast washin and a steady plateau of washout, corresponding to the shape of curves in well-perfused habitats, are more likely to be malignant [21]. Numerous studies on breast cancer tumors have focused analyses on these ‘hot spots’ in tumors [22–27], determined by looking for contiguous ROIs of high enhancement in DCE-MRI datasets. Similar analyses have been performed on soft tissue sarcoma tumors [9, 28]. However, although ‘hot spot’ methods begin to consider the heterogeneity in tumors, they only focus on well-perfused habitats, whereas identification and analysis of additional tumor habitats can provide additional discriminatory features for treatment prediction as in the present study and offer a broader understanding of circumstances in the tumor. In addition, common problems with analysis of ‘hot spots’ include the selected ROIs are often manually contoured [22–26, 28] although semi-automated [29, 30] and automated methods [27] have been proposed. Manual contours produce significant variability in measurements derived from selected ROIs [30]. Analysis of ‘hot spots’ also often precludes the inclusion of anatomical features as only the highest enhancing voxels are selected, which may not extend to the whole extent of the habitat. The present approach has the benefit of addressing the two issues by automating delineation of tumor habitats, avoiding variability with manual selection of ROIs. The approach also allows detection of habitats at a subpixel resolution [5], leading to more accurate determination of habitat locations and allowing the analysis of both habitat-specific pharmacokinetic and anatomical parameters. Specifically, in sarcoma, a similar pattern recognition approach on DCE-MRI data has been applied with good agreement between determined habitats and histopathology [31].

Significant differences were seen in perfusion, diffusion and volumetric parameters from necrotic and well-perfused habitats between responders and nonresponders. The relation of treatment response with features in necrotic habitats may be expected due to the correlation of response to necrosis in pathological samples [3]. Percent volume and volume of necrosis directly relate to the amount of necrosis in tumors, whereas k_ep and K^trans of necrotic habitats may represent how ‘necrotic’ those areas are, with lower values representing greater extent of cell death. Similarly, features associated with well-perfused regions may be expected to correlate with response based on past DCE-MRI findings. As mentioned, in breast cancer tumors, a greater proportion of well-perfused habitats have correlated with a greater chance of malignancy [21]. Since vascularity is necessary for growing tumors, well-perfused habitats would offer a window into the ‘health’ of tumors. In the case of sarcomas, lower k_ep and K^trans in well-perfused habitats indicate slower diffusion from vessels into the extracellular space, suggesting a tumor lagging in growth and potentially less active. In the present study, no significant differences in parameters from DW-MRI, namely ADC, were found between responders and responders. Given prior work linking ADC to treatment monitoring in soft tissue sarcoma patients [11–15], the lack of discriminatory ability for ADC was unexpected. A possible explanation may the use of multiple b-values across multiple MR systems in DW-MRI acquisition. Different patients would potentially have datasets with different diffusion weighting and signal-to-noise confounding ADC values between responders and nonresponders. Standardized imaging parameters would allow more definitive assessment of the value of ADC to prediction of treatment response.

The study has several limitations that are typical for such pilot investigations. It is retrospective and includes a small number of patients. The analyzed patients had different types of sarcoma and varying treatment plans and imaging exam parameters. However, despite the heterogeneity in the dataset, results are promising on the limited samples. The nature of the approach provides quantitative measures robust to the variety of cases likely to be seen in sarcoma patients. Both ADC and pharmacokinetic parameters are derived to measure true physiological characteristics of the tumor (diffusion, contrast-to-tissue exchange rates) and thus can be assumed to be relatively independent from the imaging sequences.

The appeal of the proposed approach is the non-invasive way it can identify tumor habitats in vivo. For this purpose, a pattern recognition method was utilized that was previously validated on pre-clinical [5, 32] and clinical data [5]. The method assumed the presence of three habitats in tumors which corresponded well with results seen in the pre-clinical model [5, 32] as well as with literature where three shapes of curves have been observed in tumors [21]. Finally, patients enrolled in the study had heterogeneous presentations of sarcomas, including location, tissue types, and stage of sarcomas. Although differences in the patient population could provide a confounding factor in prediction of treatment response, these differences should serve to negatively impact prediction as prediction would be more difficult on a diverse set of patients.

In conclusion, in this pilot study the feasibility for prediction of treatment response in soft tissue sarcoma patients was demonstrated based on identification and analysis of tumor heterogeneity from mpMRI. Pharmacokinetic and anatomical information extracted from unique tumor habitats provided features that could separate responders and nonresponders, the combination of which provided more discriminatory ability than a single feature. A more extensive study will be needed in the future to validate the performance of this approach. For soft tissue sarcomas which encompass a broad range of tumors, additional information derived from the state and characteristics of tumor habitats may lead to a more accurate prediction of treatment response in patients.

Funding: NEI F30-EY027162 (sponsor had no role in the study design; in the collection, analysis and interpretation of the data; in the writing of the report; and in the decision to submit the paper for publication).

References

1. In GK, Hu JS, Tseng WW (2017) Treatment of advanced, metastatic soft tissue sarcoma: latest evidence and clinical considerations. Ther Adv Med Oncol 9: 533–550.
2. Society AC (2017) Cancer facts & figures. 2017.
3. Eilber FC, Rosen G, Eckardt J, Forscher C, et al. (2001) Treatment-induced pathologic necrosis: a predictor of local recurrence and survival in patients receiving neoadjuvant therapy for high-grade extremity soft tissue sarcomas. J Clin Oncol 19: 3203–3209.
4. Gatenby RA1, Grove O, Gillies RJ (2013) Quantitative imaging in cancer evolution and ecology. Radiology 269: 8–15. [crossref]
5. Chang YC, Ackerstaff E, Tschudi Y, Jimenez B, et al. (2017) Delineation of Tumor Habitats based on Dynamic Contrast Enhanced MRI. Sci Rep 7: 9746.
6. Tschudi Y, Pollack A, Punnen S, Ford JC, et al. (2018) Automatic Detection of Prostate Tumor Habitats using Diffusion MRI. Scientific Reports 8: 16801.
7. Huang W, Beckett BR, Tudorica A, Meyer JM, Afzal A, et al. (2016) Evaluation of Soft Tissue Sarcoma Response to Preoperative Chemoradiotherapy Using Dynamic Contrast-Enhanced Magnetic Resonance Imaging. Tomography 2: 308–316. [crossref]
8. Xia W, Yan Z, Gao X (2017) Volume fractions of DCE-MRI parameter as early predictor of histologic response in soft tissue sarcoma: A feasibility study. Eur J Radiol 95: 228–235.
9. Alic L, van Vliet M, van Dijke CF, Eggermont AM, Veenland JF, et al. (2011) Heterogeneity in DCE-MRI parametric maps: a biomarker for treatment response? Phys Med Biol 56: 1601–1616. [crossref]
10. Guo J, Reddick WE, Glass JO, Ji Q, et al. (2012) Dynamic contrast-enhanced magnetic resonance imaging as a prognostic factor in predicting event-free and overall survival in pediatric patients with osteosarcoma. Cancer 118: 3776–3785.
11. Baunin C, Schmidt G, Baumstarck K, Bouvier C, et al. (2012) Value of diffusion-weighted images in differentiating mid-course responders to chemotherapy for osteosarcoma compared to the histological response: preliminary results. Skeletal Radiol 41: 1141–1149.
12. Oka K, Yakushiji T, Sato H, Hirai T, et al (2010) The value of diffusion-weighted imaging for monitoring the chemotherapeutic response of osteosarcoma: a comparison between average apparent diffusion coefficient and minimum apparent diffusion coefficient. Skeletal Radiol 39: 141–146.
13. Dudeck O, Zeile M, Pink D, Pech M, et al. (2008) Diffusion-weighted magnetic resonance imaging allows monitoring of anticancer treatment effects in patients with soft-tissue sarcomas. J Magn Reson Imaging 27: 1109–1113.
14. Hayashida Y, Yakushiji T, Awai K, Katahira K, et al. (2006) Monitoring therapeutic responses of primary bone tumors by diffusion-weighted image: Initial results. Eur Radiol 16: 2637–2643.
15. Thoeny HC, De Keyzer F, Chen F, Vandecaveye V, et al. (2005) Diffusion-weighted magnetic resonance imaging allows noninvasive in vivo monitoring of the effects of combretastatin a-4 phosphate after repeated administration. Neoplasia 7: 779–787.
16. Crombé A, Périer C, Kind M, De Senneville BD, et al. (2018) T2 -based MRI Delta-radiomics improve response prediction in soft-tissue sarcomas treated by neoadjuvant chemotherapy. J Magn Reson Imaging. [crossref]
17. Stoyanova R, Huang K, Sandler K, Cho H, et al. (2012) Mapping Tumor Hypoxia In Vivo Using Pattern Recognition of Dynamic Contrast-enhanced MRI Data. Transl Oncol 5: 437-U114.
18. Tofts PS (1997) Modeling tracer kinetics in dynamic Gd-DTPA MR imaging. J Magn Reson Imaging 7: 91–101. [crossref]
19. Tofts PS, Brix G, Buckley DL, Evelhoch JL, et al. (1999) Estimating kinetic parameters from dynamic contrast-enhanced T(1)-weighted MRI of a diffusable tracer: standardized quantities and symbols. J Magn Reson Imaging 10: 223–232.
20. Parker GJ, Roberts C, Macdonald A, Buonaccorsi GA, et al. (2006) Experimentally-derived functional form for a population-averaged high-temporal-resolution arterial input function for dynamic contrast-enhanced MRI. Magn Reson Med 56: 993–1000.
21. Kuhl CK, Mielcareck P, Klaschik S, Leutner C, et al. (1999) Dynamic breast MR imaging: are signal intensity time course data useful for differential diagnosis of enhancing lesions? Radiology 211: 101–110.
22. Dongfeng H, Daqing M, Erhu J (2012) Dynamic breast magnetic resonance imaging: pretreatment prediction of tumor response to neoadjuvant chemotherapy. Clin Breast Cancer 12: 94–101.
23. Su MY, Cheung YC, Fruehauf JP, Yu H, et al. (2003) Correlation of dynamic contrast enhancement MRI parameters with microvessel density and VEGF for assessment of angiogenesis in breast cancer. J Magn Reson Imaging 18: 467–477.
24. Pickles MD, Manton DJ, Lowry M, Turnbull LW (2009) Prognostic value of pre-treatment DCE-MRI parameters in predicting disease free and overall survival for breast cancer patients undergoing neoadjuvant chemotherapy. Eur J Radiol 71: 498–505.
25. Nagashima T, Sakakibara M, Nakamura R, Arai M, et al. (2006) Dynamic enhanced MRI predicts chemosensitivity in breast cancer patients. Eur J Radiol 60: 270–274.
26. Hayes CP, A.R.; Leach, M.O. (2002) Assessing changes in tumour vascular function using dynamic contrast-enhanced magnetic resonance imaging. NMR Biomed 15: 154–163.
27. Chen W, Giger ML, Bick U, Newstead GM (2006) Automatic identification and classification of characteristic kinetic curves of breast lesions on DCE-MRI. Med Phys 33: 2878–2887.
28. van Rijswijk CS, Geirnaerdt MJ, Hogendoorn PC, Peterse JL, et al. (2003) Dynamic contrast-enhanced MR imaging in monitoring response to isolated limb perfusion in high-grade soft tissue sarcoma: initial results. Eur Radiol 13: 1849–1858.
29. Liney GP, Sreenivas M, Gibbs P, Garcia-Alvarez R, et al. (2006) Breast lesion analysis of shape technique: Semiautomated vs. manual morphological description. Journal of Magnetic Resonance Imaging 23: 493–498.
30. Mussurakis S, Buckley DL, Horsman A (1997) Dynamic MRI of invasive breast cancer: assessment of three region-of-interest analysis methods. J Comput Assist Tomogr 21: 431–438.
31. Venianaki M, Salvetti O, de Bree E, Maris T, et al. (2018) Pattern recognition and pharmacokinetic methods on DCE-MRI data for tumor hypoxia mapping in sarcoma. Multimedia Tools and Applications 77: 9417–9439.
32. Stoyanova R, Huang K, Sandler K, Cho H, et al. (2012) Mapping Tumor Hypoxia In Vivo Using Pattern Recognition of Dynamic Contrast-enhanced MRI Data. Transl Oncol 5: 437–447.