Cluster analysis of dynamic contrast enhanced MRI reveals tumor subregions related to locoregional relapse for cervical cancer patients

Abstract

Background: Solid tumors are known to be spatially heterogeneous. Detection of treatment-resistant tumor regions can improve clinical outcome, by enabling implementation of strategies targeting such regions. In this study, K-means clustering was used to group voxels in dynamic contrast enhanced magnetic resonance images (DCE-MRI) of cervical cancers. The aim was to identify clusters reflecting treatment resistance that could be used for targeted radiotherapy with a dose-painting approach.

Material and methods: Eighty-one patients with locally advanced cervical cancer underwent DCE-MRI prior to chemoradiotherapy. The resulting image time series were fitted to two pharmacokinetic models, the Tofts model (yielding parameters K^trans and ν_e) and the Brix model (A_Brix, k_ep and k_el). K-means clustering was used to group similar voxels based on either the pharmacokinetic parameter maps or the relative signal increase (RSI) time series. The associations between voxel clusters and treatment outcome (measured as locoregional control) were evaluated using the volume fraction or the spatial distribution of each cluster.

Results: One voxel cluster based on the RSI time series was significantly related to locoregional control (adjusted p-value 0.048). This cluster consisted of low-enhancing voxels. We found that tumors with poor prognosis had this RSI-based cluster gathered into few patches, making this cluster a potential candidate for targeted radiotherapy. None of the voxels clusters based on Tofts or Brix parameter maps were significantly related to treatment outcome.

Conclusion: We identified one group of tumor voxels significantly associated with locoregional relapse that could potentially be used for dose painting. This tumor voxel cluster was identified using the raw MRI time series rather than the pharmacokinetic maps.

A wide range of studies suggests that tumor heterogeneity could become an important biomarker for tumor grading and treatment outcome [1,2]. There is spatial heterogeneity within each tumor, as well as heterogeneity between tumors. Intratumor heterogeneity suggests regions with different physiological properties within the tumor, which may respond differently to treatment. Early identification of such regions can be used to visualize treatment-resistant areas and could therefore be useful in planning targeted therapy.

Functional imaging, such as dynamic contrast enhanced magnetic resonance imaging (DCE-MRI), gives insight into the tumor microenvironment and allows us to assess the heterogeneity within the tumor [1,3]. Pharmacokinetic models can be fitted to DCE-MRI data to ease interpretation of the dynamic image series and to link the information captured in the images to physiological parameters. Two models often used in studies of cervical cancer patients are the Tofts and the Brix models [4,5]. These models are based on different underlying assumptions. The Tofts model requires a measurement of the arterial input function, whereas the Brix model assumes a monoexponential venous clearance of contrast agent. Both Tofts and Brix parameters have been shown to correlate with treatment outcome for cervical cancer [6,7].

Unsupervised clustering techniques such as the K-means clustering algorithm can be used to group voxels according to similarity, without requiring any prior assumptions or knowledge about the voxels [8]. The K-means algorithm is fast and robust, and is appropriate for large data sets, consisting of thousands to millions of voxels [8]. K-means clustering and similar methods may be used in segmenting MR images of the brain [9], but have also been utilized for segmentation of cervical cancers [10]. In a previous study, we used K-means analysis of Tofts parameter maps to identify a voxel cluster associated with treatment outcome [11]. As this cluster was associated with local control, not relapse, it was not a good candidate for targeted therapy. Similar cluster analysis based on Brix parameter maps is particularly appealing, as the Brix amplitude parameter A_Brix is known to reflect hypoxia in cervical cancer [12,13]. Unprocessed DCE-MRI data may also be attractive for clustering. Thus, a more comprehensive analysis of cluster features based on DCE-MRI for cervical cancer is warranted.

In this study, we investigated the association between tumor voxel clusters and treatment outcome for patients with locally advanced cervical cancer. Tumor voxels with similar properties were grouped into clusters using the K-means algorithm. Clustering was based on the unprocessed DCE-MRI time series or on physiological parameters derived from the Tofts or Brix models. The cluster volume fractions as well as the spatial distribution of the clusters were then tested for associations with treatment outcome. The significant voxel clusters were evaluated for their potential to be used in treatment planning.

Material and methods

Patients

Eighty-one patients with locally advanced cervical cancer underwent DCE-MRI prior to chemoradiotherapy (Table 1). We obtained written informed consent from all patients, and the Regional Committee for Medical Research Ethics in southern Norway approved the study.

Table 1. Summary of the patient cohort. Some statistics are reported as median value (range).

Number of patients	81
FIGO stage
IB:	2
II:	44
III:	29
IVA:	6
Pathologic lymph nodes
Present	31
Absent	50
Age (years)	57 (32–85)
Tumor volume (cm^3)	32 (1–160)
Progression-free survival	49
Relapses (total)	32
Locoregional	9
Metastatic	28
Follow-up (months)	62 (3–149)
Time to event (months)
Locoregional	11 (5–23)
Metastatic	14 (3–65)

All patients received chemoradiotherapy with curative intent. External beam radiotherapy was given over five weeks (25 fractions; 2 Gy/fr to tumor, 1.8 Gy/fr to pelvis). In addition, each patient was given 5–6 fractions of intracavitary brachytherapy (4.2 Gy/fr) and cisplatin chemotherapy (40 mg/m² per week).

The clinical end-point was locoregional control, defined as control within the irradiated pelvis. As Table 1 shows, 32 patients relapsed, nine with local failure and 28 with metastasis. Five patients had both locoregional and distant failure.

Imaging

MRI was performed using a 1.5 T Signa Horizon LX tomograph (GE Medical Systems, Milwaukee, WI, USA) with a pelvic phased array coil. For tumor delineation, we used an axial T2-weighted fast spin echo sequence (FSE; TR, 4220 ms; TE, 85 ms; matrix size, 512 × 512; field of view, 200 mm; slice thickness, 5 mm; slice gap, 1 mm). The DCE-MRI series was acquired using an axial T1-weighted fast spoiled gradient recalled sequence (FSPGR; flip angle, 90°; echo time, 3.6 ms; repetition time, 80 ms; matrix size, 256 × 256; field of view, 200 mm; slice thickness, 5 mm; slice gap, 1 mm). The contrast agent Gd-DTPA (Magnevist^®; Schering, Berlin, Germany) was administered manually as a fast bolus injection with a total dose of 0.1 mmol/kg body weight followed by a saline flush.

The DCE-MRI series consisted of one pre-contrast image and up to 13 images acquired after contrast agent injection. The sampling interval was 15 seconds for the first 11 images and one minute for the last two images.

Pharmacokinetic analysis

Voxel clustering based on the Tofts model [5] and the Brix model [4] was compared. Details about these models, as well as the model fitting, can be found in the Supplementary Materials (available online at http://www.informahealthcare.com). Both of these models describe the relative signal increase (RSI) as a function of the time t after contrast injection. The Tofts model yields two parameters. K^trans is the transfer constant of the contrast agent from the blood to the extravascular extracellular space (EES) and ν_e is the volume fraction of the EES in the voxel. The Brix model uses three parameters. A_Brix is the amplitude, k_ep is the transfer rate of contrast agent from tissue to blood and k_el is the wash-out rate of contrast agent from the blood [4]. Voxels with very high unphysiological values of the pharmacokinetic parameters, A_Brix >10, k_ep >12 min⁻¹ and/or k_el >0.5 min⁻¹, were excluded from further analysis after the Brix model fitting. About 0.3% of the voxels had unphysiologically high RSI values (RSI >8). These values were considered as outliers, and were removed before the cluster analysis of the RSI maps.

K-means cluster analysis

K-means cluster analysis is an unsupervised method for dividing data into groups, or clusters, with similar properties [8,14]. In this study, we clustered the image voxels based on the pharmacokinetic parameters (Tofts model: K^trans, ν_e; Brix model: A_Brix, k_ep, k_el) or on the RSI time series. The clustering algorithm used the squared Euclidian distance, and 30 replicates with different initializations. Initializations had only very minor effects on the resulting clusters (data not shown). The data was scaled to mean zero and standard deviation one before the cluster analysis. The optimal number of clusters was found, as previously described in [11], by comparing the mean intra-cluster distance with the minimum inter-cluster distance [14]. To check the robustness of the clustering, we divided the voxels into clusters based only on the central slice of each tumor, and compared the resulting cluster centroids with the centroid of the clusters based on all tumor voxels in all slices.

Spatial distribution of voxel clusters

The spatial distribution of the voxel clusters within each tumor was analyzed by comparing the number and size of connected areas or patches formed by each cluster. We examined the spatial distribution of clusters in the central slice of each tumor. For each cluster, we counted the number of connected areas (NumArea) formed by the voxels belonging to the cluster. We also calculated the median size (number of voxels) of the connected areas (SizeArea). These two measures were used to assess the heterogeneity of clusters within the tumors. To enable comparison between tumors, both SizeArea and NumArea were normalized by dividing by the total number of voxels in the tumor slice.

Statistical analysis

The volume fraction of each voxel cluster was calculated for each tumor. The patients were then divided in two groups based on the median cluster volume fraction (CVF). In addition, we grouped the patients according to the median value of the spatial parameters SizeArea and NumArea. Differences in treatment outcome between the groups were estimated by log-rank tests and visualized using the Kaplan-Meier estimator. All p-values reported are two-sided and we used a significance level of 0.05 after Bonferroni correction. The association between tumor volume and CVF or the spatial parameters was examined with Pearson’s correlation analysis.

Calculations

The log-rank tests were conducted in R (v3.2.0) using the survival library. All other analyses were conducted in MATLAB (v8.3, R2014a, The Mathworks Inc., Natick, MA, USA). For the K-means cluster analysis and the correlation tests, we used the MATLAB Statistics Toolbox (v9.0). The spatial distribution features were calculated using the MATLAB Image Processing Toolbox (v9.0).

Results

The optimal K-means clustering based on the Tofts parameters divided the voxels into one large and two smaller clusters, as shown in Figure 1 (top row) for two tumors with different treatment outcomes. With the exception of one patient, who had no voxels belonging to the Tofts3 cluster, all tumors contained voxels from all three clusters. Table 2 shows the centroid of each K-means cluster, as well as the median CVF. The largest voxel cluster, Tofts2, had a median volume fraction of 62%. Although initial significance screening suggested significant relationship between Tofts1 and locoregional relapse (nominal p-value 0.04), this was not significant after correction for multiple testing.

Figure 1. Tofts (top), Brix (middle) and RSI (bottom) voxel clusters in the central slice of tumors from two different patients. These patients had no relapse (left) or locoregional relapse (right), respectively. The white areas in some of the slices are voxels excluded due to unphysiological parameter values.

Table 2. K^trans and ν_e values for the centroids of the three Tofts-based voxel clusters. The last column gives the median cluster volume fraction (CVF) of the cluster in each tumor.

Cluster	K^trans (min^−1)	νe	Median CVF
Tofts1	0.09	0.91	0.07
Tofts2	0.13	0.29	0.62
Tofts3	0.34	0.57	0.18

The optimal clustering based on the Brix parameters resulted in one relatively large cluster, Brix3, with a median CVF of 66% and three smaller clusters (Figure 1, Table 3). Most tumors contained voxels from all four clusters. The exception was the Brix1 cluster, which was not represented in two tumors, and the Brix4 cluster, which was not seen in three tumors. After Bonferroni correction, none of the Brix clusters were significantly related to treatment outcome.

Table 3. A_Brix, k_ep and k_el values for the centroids of the four Brix-based voxel clusters. The last column gives the median cluster volume fraction (CVF) of the cluster in each tumor.

Cluster	ABrix	kep (min^−1)	kel (min^−1)	Median CVF
Brix1	1.61	5.94	0.06	0.09
Brix2	2.10	1.00	0.25	0.08
Brix3	1.84	2.05	0.04	0.66
Brix4	5.74	0.91	0.18	0.02

Clustering based on the RSI data resulted in two clusters shown in Figure 1. The time curves representing the centroids of these clusters, shown in Figure 2, indicated that RSI1 was a cluster of low-enhancing voxels, while RSI2 contained high-enhancing voxels. The RSI1 cluster was the largest, with a median CVF of 58%, while the RSI2 cluster had a median CVF of 42%. With the exception of one patient, all tumors contained voxels from both clusters. This tumor had only RSI1, i.e. low-enhancing voxels. The number of connected areas (NumArea) of the low-enhancing RSI1 cluster was significantly related to locoregional control (adjusted p-value 0.048, Figure 3). Patients with few connected RSI1 areas (NumArea), that is a less patchy distribution of this cluster, had poorer treatment outcomes compared to patients with many low-enhancing patches scattered throughout the tumor (Figure 1). The NumArea parameter for RSI1 had a significant negative correlation to total tumor volume (p-value <0.001), indicating that larger tumors had fewer low-enhancing patches relative to smaller tumors.

Figure 2. RSI time curves representing the centroids of the two RSI voxel clusters. RSI1 (solid line) was a cluster of low-enhancing voxels, while RSI2 (dashed line) consisted of high-enhancing voxels. The median volume fraction of these voxel clusters was 58% for RSI1 and 42% for RSI2.

Figure 3. Plot of the Kaplan-Meier estimators for the tumor voxel cluster showing significant differences in treatment outcome between patient groups. Clustering was based on the RSI values. NumArea was the number of connected areas formed by this cluster.

We also performed log-rank test after dividing the patients into two groups according to tumor volume. The tumor volume was not significantly associated with locoregional control (p-value 0.51).

Discussion

Previous work has suggested a potential role of Tofts, Brix and RSI parameters in the treatment planning of cervical cancer [6,15]. Our work may facilitate this strategy by proposing a method for visualizing treatment-resistant regions within the tumors. In this study, we identified a tumor voxel cluster in DCE-MR images that was significantly related to locoregional relapse for patients with cervical cancer. This cluster was based on the RSI, not on the pharmacokinetic maps.

The RSI-based cluster significantly related to locoregional relapse consisted of voxels with low contrast enhancement. Tumor regions of low enhancement most likely reflect areas of low perfusion, low oxygen concentration and thus high radioresistance. These characteristics have previously been found to be negatively associated with outcome [15,16], in line with the findings of this study. The low-enhancement cluster was also found to be significantly correlated to tumor volume. However, tumor volume was not significantly related to locoregional control. Therefore, the significant association between the low-enhancing cluster and locoregional control was not due to correlation with tumor volume, but rather the physiological characteristics of the cluster.

Despite the significance of the RSI-based cluster, none of the Tofts or the Brix clusters were significantly associated with outcome after correction for multiple testing. This suggests that analyzing the raw RSI data retains information related to outcome that may have been lost by pharmacokinetic model fitting. The choice of pharmacokinetic model is much debated [17,18], but our results indicate that pharmacokinetic modeling can be bypassed and the RSI values can be used directly.

It should be noted that Andersen et al. [11] also found three Tofts clusters when analyzing this dataset, but that these clusters are not identical to the ones identified in the current study, due to different scaling of the data. Andersen et al. found that the Tofts-based clusters were mainly separated by ν_e, and not K^trans, while the Tofts-based clusters in the current study were separated by both pharmacokinetic parameters, as shown in Table 2. The data scaling of the current study thus ensured that both Tofts parameters were reflected in the voxel clusters.

The spatial distribution of the voxel clusters is important if they are to be used for targeted radiotherapy in the form of dose painting [19]. Such tailored treatment will be easier to implement for a continuous tumor region than for more scattered clusters [20]. Voxel clusters related to locoregional relapse would be most relevant to use for dose painting, as such clusters would indicate radiotherapy-resistant areas in the tumor [19]. Of the clusters identified in the current study, the low-enhancing RSI1 cluster appears to be the best candidate for dose painting. In tumors with poor prognosis the low-enhancing RSI1 voxels were gathered in few patches, which potentially could receive a radiation boost. Regions showing low contrast enhancement have in fact been proposed as targets for dose painting [21]. Our study confirms that these regions are indeed linked to locoregional relapse.

In conclusion, voxels from pretreatment DCE-MRIs of cervical cancers were grouped into clusters, one of which was related to locoregional recurrence. This cluster, based on the RSI time series, could potentially be used for targeted radiotherapy.

Conflicts of interest

The authors report no conflicts of interest. The authors alone are responsible for the content and writing of the paper.