Abstract
A novel post-processing methodology able to assess whole-body tumor heterogeneity in patients with metastatic disease is proposed. The method is demonstrated on paired pre- and post-treatment data sets obtained from an initial cohort of six patients with metastatic disease from primary prostate or ovarian cancers. Whole-body diffusion-weighted imaging and T1-weighted contrast-enhanced imaging data were acquired covering the chest, abdomen, and pelvis. Joint histograms of Apparent Diffusion Coefficient and Fractional Enhancement values were calculated within volumes of interest and were modeled as a Gaussian mixture of two classes. Probability maps and volumetric estimates of the magnetic resonance data-derived classes providing visualization of pre- and post-treatment data are shown in three patient examples. This technique provided spatially heterogeneous characterization of regions following treatment as defined by the combined analysis of apparent diffusion coefficient and fractional enhancement. A new whole-body magnetic resonance data analysis has been demonstrated enabling visualization of intra-patient response heterogeneity in patients with metastatic cancer. Changes in the parameters of each subpopulation derived from this technique (apparent diffusion coefficient and fractional enhancement) reflect changes in the tissue properties of each subpopulation following treatment. Furthermore, the volume change of each population can be quantified. Such techniques may be essential for personalized anti-cancer therapy where there is a need to detect early drug-resistance and monitor heterogeneous response.
Keywords
Introduction
Tumor heterogeneity relates to biological differences that may exist within and between tumors (intra- and inter-patient heterogeneity).1–3 These differences arise as a result of clonal evolution in the genetic and micro-environmental characteristics of tumor-cell subpopulations, as the cancers grow and metastasize.2–4 As oncologic practice evolves toward more precise and individualized therapy, consideration of tumor heterogeneity is gaining importance when assessing tumor response to treatment.5,6 More precise treatment targeting will rely on identification of the specific disease phenotype as undetected drug-resistant subclones can lead to disease that is unresponsive to treatment. Current practice utilizes surgical resection, sequential and/or multi-region biopsies to study and demonstrate tumor heterogeneity. However, these techniques are invasive, and biopsies are also prone to sampling bias, which may underestimate heterogeneity within and between tumors.1,7
Imaging techniques are non-invasive and can be performed in vivo at high-spatial resolution, making them ideal tools to study tumor heterogeneity, both prior to and during therapy. 8 Magnetic resonance imaging (MRI) is attractive as it does not incur ionizing radiation and can be safely repeated for sequential examinations. Furthermore, a range of functional imaging techniques (e.g. dynamic contrast-enhanced MRI [DCE-MRI], diffusion-weighted MRI [DWI] and intrinsic susceptibility MRI [IS-MRI]) can be applied to inform on biologically relevant tumor characteristics such as vascularity, cellularity and tissue oxygenation. In addition, recent MR system hardware developments have enabled extended imaging coverage over a wide area of the body within a clinically acceptable time frame, making it feasible to perform anatomical and functional whole-body MRI to evaluate patients with cancers involving both bones and soft tissues.9–12
The vast majority of patients with metastatic castration-resistant prostate cancer (mCRPC) develop metastasis in bone. Bone biopsies are painful, not always accessible and are subject to sampling variation in regions of heterogeneity. Patients with metastatic ovarian carcinoma exhibit peritoneal and liver metastasis beyond the locally advanced disease in the pelvis.
In this manuscript, we propose a novel image analysis strategy for quantification and visualization of whole-body response heterogeneity in an exploratory cohort of six patients with metastatic ovarian and prostate cancers.
We focus our attention on the analysis of whole-body MRI examinations that combine DWI11,12 with a semi-quantitative T1-weighted contrast-enhanced imaging protocol. Both techniques provide voxel-wise quantification of the tumor microenvironment via co-registered calculation of the apparent diffusion coefficient (ADC) from DWI and the fractional enhancement (FE) from contrast-enhanced studies. Previous pre-clinical studies have demonstrated the correlation between increase in tumor ADC and cell-kill following successful treatment. 13 Furthermore, decrease in FE may reflect damage to the tumor vasculature. 14 By observing contemporaneous changes of these parameters by two-dimensional histogram analysis, we can visualize tumor response that conforms to the expected treatment behavior versus those areas that behave differently, thereby providing insights into response heterogeneity. We evaluate the ability of our new approach to identify clusters of voxels with defined characteristics, visualizing the spatial distributions in the body, and demonstrate the approach using paired pre- and post-treatment whole-body MRI data in six patients (a mix of responding and non-responding patients) following the same treatment with a novel-targeted agent. The analysis reported here is directed at exploring this newly proposed methodology for visualizing treatment effect.
Materials and methods
Patients
Six patients with metastatic ovarian and castrate resistant prostate cancers enrolled in a phase II clinical trial testing, a novel targeted agent were reviewed before and after treatment. The whole cohort consisted of two women and four men, with age ranging from 48 to 70 years (mean age 64.2 years). All patients were scanned prior to and at six weeks after the start of treatment. In addition to their primary cancer location, patients had bone or liver metastases and nodal involvement.
All patients gave written informed consent, and the local ethical committee approved the research study.
Two patients (one ovarian and one prostate cancers) showed benefit from therapy as evidenced by clinical and biochemical parameters, while the remaining four patients demonstrated no significant benefit from therapy.
MR acquisition
All data were acquired on a 1.5T Avanto scanner (Siemens, Erlangen, Germany) using a neck coil plus two body-array coils to fully cover the chest, abdomen and pelvis of the patient. The MR protocol consisted of a sequential acquisition of three axial stations from chest to pelvis, repeated for each of the following sequences: (a) DWI, (b) T1-weighted pre-Gadolinium contrast administration, and (c) T1-weighted post-Gadolinium contrast (see Figure 1).
Illustration of the order of MR image acquisition used in studies. Patients were scanned over three contiguous stations covering the chest (1), abdomen (2), and pelvis (3) for each sequence: DWI, T1-w pre-contrast, and T1-w post-contrast. The effect of the delay after Gadolinium administration on FE values estimated at each of the three stations was assessed visually by comparing FE distributions at each station individually (data not shown). The variation between stations was observed to be smaller than that between tissue populations, suggesting our results were not affected by these acquisition delays.
Diffusion sequence
DWI was performed during free-breathing using a twice refocused spin echo Inversion recovery 2D echo planar imaging with the following parameters: 46 contiguous slices per station, voxel resolution = 2.5 × 2.5 × 5 mm3, number of excitations (NEX) = 4, field of view (FoV) = 380 × 380 mm2, TR = 12,700 ms, Ti = 160 ms, TE = 69 ms, 6/8 partial Fourier encoding in the phase-encode direction, Grappa acceleration factor 2 (32 reference lines), bandwidth = 1960 Hz/Px, b values = 50/900 s/mm2 acquired as a three scan trace, acquisition time ∼6 min/station. The parameter of interest derived from these data was the ADC, calculated assuming a mono-exponential model:
T1-weighted sequence (pre- and post-Gadolinium)
T1-w images were acquired using a 3D gradient echo with the following parameters: 52 contiguous slices per imaging station, voxel resolution = 1.6 × 1.6 × 5 mm3, NEX = 1, FoV = 400 × 300 mm2, TR/TE = 5.68/2.65, 24° flip angle, SPAIR fat suppression, Grappa acceleration factor 2 (32 reference lines), bandwidth = 250 Hz/Px, acquisition time = 21 s/station. The patients were instructed to hold their breath during the acquisition of each station. A dose of contrast agent (Dotarem) of 0.2 ml/kg at 2 ml/s flow rate, followed by 20 ml of saline flush at the same flow rate, was delivered by power injector. The post contrast images were acquired at specific delayed time-points for each station following the first pass of contrast, i.e. 30 s (chest), 65 s (abdomen) and 120 s (pelvis) after injection (see Figure 1). Following acquisition, the T1-weighted data were re-sampled to match the resolution of the DWI data using the OsiriX image-processing suite,
15
and a FE parameter map was calculated:
Image segmentation
Regions of interest were defined for disease by an experienced radiologist with five years of experience in body DW-MRI using computed high b-value images with the aid of a previously described segmentation technique.16,17 Briefly, the radiologist visually selected a computed b-value that maximized the contrast between disease and background tissues on Maximum Intensity Projection displays of the Whole-Body Diffusion-Weighted Imaging (WBDWI) data. An initial threshold was then manually selected to provide maximal contrast between disease and background classes by the expert reader. The classification was then smoothed using a Markov random field prior, after which the resulting tumor Volumes Of Interest (VOIs) ICR2 were manually corrected if necessary. Segmentation of patient data was performed before and after treatment by the same individual to avoid inter-operator variation. All VOIs were stored and transferred onto calculated ADC and FE maps to enable voxel-wise statistics of ADC and FE to be recorded (Figure 2(a)).
Outline of image analysis steps (prostate patient, axial images). (a) Volumes of interest (VOIs – shown in green) are transferred from WBDWI data to the resolution-matched FE maps. (b) Co-registered FE and ADC voxels within the VOIs are visualized as a two-dimensional scatter plot. GMM is used to infer two subpopulations from the data (isocontours represent the logarithm of the fitted marginalized probability distribution, given by equation 6). The results from GMM fitting (c) may be used to calculate the classification probability of each voxel, visualized as a color scale overlaid on the DWI images for each visit (d).
Image analysis
Our aim was to segregate the derived joint distributions of ADC and FE into distinct subpopulations that may be biologically meaningful to reflect treatment response heterogeneity within and between tumors (i.e. cellularity as characterized by ADC and vascularity characterized by the FE values). For example, we expected viable tumors to return lower ADC values and higher FE values. The converse should be observed in responding treated tumors. We modeled joint histograms of ADC and FE of the total patient data set (i.e. pre- and post-treatment combined) using a Gaussian Mixture Model (GMM) of j = 1 … M subpopulations, each of which may be characterized by a normal distribution N with mean vector µj and covariance matrix Σj. Given that a voxel i comes from subpopulation j, the probability of its value yi = (ADCi, FE
i
) is given by
The marginalized probability of the values at each voxel is then given by
This classification can be visualized as a color scale to identify regions of heterogeneous tumor environment (see Figure 2).
The image processing steps taken in our image analysis approach for quantification and visualization of response heterogeneity are demonstrated in Figure 2:
Pre- and post-contrast-enhanced images are interpolated to match the resolution and FoV of the DW-imaging studies using the OsiriX image-processing suite.
15
The interpolated data are used to calculate FE maps according to equation (2), which are assumed to be inherently co-registered with the calculated ADC maps (this could be tested via visual inspection of fused post-contrast and low b-value images). Regions of interest defined on the DW-images are translated onto the FE maps (Figure 2(a)). After visual inspection of the two dimension histogram, these data were modeled by a two-Gaussian Mixture (M = 2) using the EM algorithm
18
(Figure 2(b) and (c)). The posterior probability for classification (equation (6)) is calculated at each voxel location and visualized on a color scale that may be overlaid on reference anatomical images (Figure 2(d)). Volumetric estimates of each subpopulation can be obtained and visualized by classifying each voxel i with the labeling that maximizes the posterior probability as per equation (6).
These methods provide three quantitative indices that may be used to monitor heterogeneous change throughout treatment: (i) changes in mean ADC provide an indication of cellular density therapeutic effects on each subpopulation, (ii) changes in calculated FE values indicate the effects on the vascular properties of the tumor and (iii) changes in the proportional volume of each subpopulation indicate tumor shrinkage and/or growth following treatment.
Results
The proposed methodology was successfully applied in all patient data sets; Table 1 presents the changes in ADC and FE for both of the sub-components in all patients, while Table 2 displays the changes in the proportional volume. Typical examples are demonstrated in Figures 3–5.
Image analysis for an ovarian patient (peritoneal metastases) who showed clinical improvement following treatment (Patient 2, axial acquisitions, coronal reconstruction). GMM successfully resolves two subpopulations before and after treatment. One class (red) is represented by areas of low ADC/relatively high FE, characterizing viable disease in both pelvic and peritoneal nodes. Another class (green) shows areas of low FE and high ADC indicating regions of lysis in the pelvic tumor (see rightmost images). The post-treatment analysis clearly demonstrated treatment response in regions of viable disease (FE decrease and volume shrinkage of red areas). Image analysis for a prostate patient (bone metastases) with no clinical improvement following treatment (Patient 6, axial acquisitions, coronal reconstruction). GMM successfully resolves two subpopulations before and after treatment. Little change was observed in the classification map and in the relative volume of each class following treatment. Furthermore, the mean values of each subpopulation show little change indicating non-responding/stable disease following therapy. Changes in mean apparent diffusion coefficient (ADC) and fractional enhancement (FE) following treatment for both estimated classes in all patient cases. Note that small changes in ADC are observed in all cases except for Patient 2. mCRPC: metastatic castrate resistant prostate cancer; mOC: metastatic ovarian cancer. Changes in proportional tumor volume in all patient cases. Note the dramatic changes in proportional volume for patients 2 and 4, both of whom demonstrated clinical response following treatment.

Figure 3 demonstrates the classification for a responding patient diagnosed with metastatic ovarian cancer (Patient 2). Before treatment, two components can clearly be seen; one color-coded red to represent regions of solid tumor (low ADC and high-contrast uptake indicated by high FE) and another color-coded green to represent regions of tumor lysis (indicated by an ADC close to that of free water and minimal contrast uptake). Following treatment, there is a 20.5% reduction in the proportional fluid volume in the tumor (see Table 2) suggesting clearance of this component, which is in agreement with an observed reduction of measured ADC (−0.59 × 10−3 mm2/s, Table 1). We also observed a reduction in the calculated FE for the viable tumor component (red) following treatment. This indicates degradation of the vascular network to the tumor in these regions, which is in line with the observed clinical response for this patient.
Figure 4 presents an example of GMM analysis resolving two subpopulations before and after treatment for a prostate patient with bone disease and nodal involvement (Patient 4). One population demonstrates wide variation in FE and low ADC (red), suggestive of disease in this region. The second subpopulation demonstrated relatively wide variation in ADC and low FE values (green), indicating increased tumor necrosis and weaker vascular supply in these regions. Following treatment, we observed a reduction in the proportion of voxels within the first subpopulation and an increase in the proportion attributed to the second. Through the use of class probability maps (Figure 4, center column), it was demonstrated that this class-shift was mostly observed within the bone, while no significant changes were seen for the nodal site or prostate, clearly suggesting differential response for this patient.
Image analysis for a prostate patient (bone metastases and pelvic nodes) who showed clinical improvement following treatment (Patient 4, axial acquisitions, coronal reconstruction). GMM successfully resolves two subpopulations before and after treatment; one represented by areas of low ADC and higher FE (red), and another by areas of low FE and higher ADC (green). Before treatment, the majority of voxels belong to the red class, both in the bone and nodal sites. After treatment, we observed an increase in the proportion of voxels represented by the second class (green). This change was visualized to be primarily in the bone (center column), while no changes were seen for the nodal site or prostate. This was highly suggestive of differential response in this patient.
In contrast, Figure 5 shows a prostate cancer patient (Patient 6) for which we observed little change in the measured values of ADC/FE (Table 1) or the tissue classification map, indicating non-responding, stable disease following therapy. Similar results were observed for three other patients, all of whom were classified as clinically not responding to treatment.
Discussion
Clonal evolution of cancer remains a major challenge for the successful management of patients using targeted and conventional therapies.3,19 The ability to accurately determine and distinguish regions that are responding versus those that are resistant to therapy would help to design better drug regimens, especially in patients who have had multiple lines of drug treatment since conventional size measurement criteria are often unhelpful assessing disease status in the presence of differential tumor response.
Imaging can provide non-invasive quantitative measurements to describe pathophysiological changes to the tumor in response to treatment. In this paper, we employed whole-body MRI combining DWI and contrast-enhanced imaging to quantify the tissue ADC and fractional tissue enhancement (FE), which are believed to reflect tumor cellularity and tumor vascularity, respectively.13,14 Through GMM of the joint distribution of voxel ADC and FE values in regions of suspect malignancy, we are able to classify biologically distinct sub-regions throughout the body and monitor changes to each sub-region both visually and quantitatively.
This technique was successfully applied to a small cohort of six patients, treated with a novel targeted therapy. In cases of metastatic bone cancer, it is not possible to monitor post-therapeutic tumor response using conventional size measurement criteria as bone metastases are considered ‘unmeasurable’. 20 Our method provides not only quantification of biological changes occurring to bone as a result of treatment but also provides fully volumetric assessment of tumor burden, visualization and characterization of heterogeneous tumor response in vivo. Although previous studies have addressed the issue of tumor response heterogeneity through either single MR parameter methods (e.g. functional diffusion maps 21 ) or by multispectral analysis (MSA) of multi-parametric MRI,22–24 these methods either relied on good registration of post-treatment data sets or were restricted to limited fields of view. Our technique provides whole-body estimation in a clinically feasible time frame and does not rely on accurate registration of pre- and post-treatment data sets.
In this paper, we have chosen to fit two subpopulations in the examples based on visual inspection of the multi-parametric histograms, and it should be noted that this may not be the general case. Other techniques such as non-parametric inference may be advantageous in this respect. 25 Furthermore, the cases presented in this paper are used to illustrate the potential of this technique, and future studies should include larger patient cohorts to validate this approach and explore the repeatability of the derived parameters. A caveat of our proposed technique is the assumption of good registration between the T1-w and DWI acquisitions, which may be problematic in regions that are subject to substantial respiratory motion such as the liver, which was excluded in the current study. However, given the large number of voxels available to whole-body MRI, we believe that motion will have minimal impact on the final results, especially when evaluating skeletal disease, which is static.
In conclusion, we describe a new multi-modal image processing methodology that can be used to provide quantification and visualization of the heterogeneity of tumors across the body and can be applied to both skeletal and soft tissue disease. Changes in the ADC and FE values inferred from MR imaging may reflect biological changes occurring in the tissue subpopulations following therapy. The volume change of each subpopulation can also be quantified providing clinicians with a treatment response biomarker that takes account of spatial heterogeneity. Understanding response heterogeneity is increasingly important in the development of more precise and personalized targeted cancer therapy. The methods described in this article could be easily be extended to incorporate other imaging modalities such as PET/CT (provided adequate spatial alignment of data sets can be achieved) equipping oncologists with new metrics for monitoring the response of patients to targeted therapies.
Footnotes
Acknowledgments
The authors are grateful for the patients’ availability and cooperation during the course of this study. MOL is a NIHR Senior Investigator. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. This paper presents independent research funded by the National Institute for Health Research (NIHR). The views expressed are those of the authors and not necessarily those of the NHS, the NIHR or the Department of Health.
Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by funding from CRUK and EPSRC support to the Cancer Imaging Centre at ICR and RMH in association with MRC & Dept. of Health C1060/A10334, C1060/A16464 and NHS funding to the NIHR Biomedicine Research Centre and the Clinical Research Facility in Imaging; An Experimental Cancer Medicine Centre Network award (joint initiative, CRUK & UK Department of Health) [grants C51/A7401 & C12540/A15573]). Support was also received from EPSRC Platform Grant EP/H046526/1, NIHR postdoctoral fellowship NHR011X (MDB) and the Spanish Medical Oncology Society “Beca SEOM para la Investigación Traslacional en el Extranjero” (DL).
