Diffusion-weighted imaging-based probabilistic segmentation of high- and low-proliferative areas in high-grade gliomas.

The apparent diffusion coefficient (ADC) derived from diffusion-weighted imaging (DWI) correlates inversely with tumor proliferation rates. High-grade gliomas are typically heterogeneous and the delineation of areas of high and low proliferation is impeded by partial volume effects and blurred borders. Commonly used manual delineation is further impeded by potential overlap with cerebrospinal fluid and necrosis. Here we present an algorithm to reproducibly delineate and probabilistically quantify the ADC in areas of high and low proliferation in heterogeneous gliomas, resulting in a reproducible quantification in regions of tissue inhomogeneity. We used an expectation maximization (EM) clustering algorithm, applied on a Gaussian mixture model, consisting of pure superpositions of Gaussian distributions. Soundness and reproducibility of this approach were evaluated in 10 patients with glioma. High- and low-proliferating areas found using the clustering correspond well with conservative regions of interest drawn using all available imaging data. Systematic placement of model initialization seeds shows good reproducibility of the method. Moreover, we illustrate an automatic initialization approach that completely removes user-induced variability. In conclusion, we present a rapid, reproducible and automatic method to separate and quantify heterogeneous regions in gliomas.

Introduction

Diffusion-weighted imaging (DWI) has been extensively used as an imaging method that yields surrogate parameters for tumor cellularity and proliferation and may be used to monitor treatment response in the brain and the abdomen[]. The apparent diffusion coefficient (ADC) value depicts the average diffusion or the Brownian molecular motion of water within a voxel[,]. A low ADC correlates, among many other factors, with a high cellular density and increased proliferation, whereas a high ADC is related to low cellular density and low proliferation[,,].

Since quantitative measures of the ADC in tumor regions may be of great clinical importance both in initial tumor grading as well as for therapy follow-up, such measures must be extracted from the data in a reproducible fashion. Commonly, DWI data are evaluated using manual region of interest (ROI) analysis. Using such analysis, a negative correlation between the ADC and the micro vessel density, as a marker of tumor proliferation was found in the tumor bulk[]. In another study, the minimum ADC in astrocytic tumors (all WHO Grades I–IV) was found to be in general negatively correlated to the Ki-67 labeling index (LI), which is an index for cell proliferation[]. However, since the borders within heterogeneous tumor regions on the one hand, and between cerebrospinal fluid (CSF) and necrotic areas on the other hand, are ill defined, the ROI placement is inherently user dependent. Thus, an objective automated method for tissue segmentation beyond manual ROI placement is needed to enable a reliable ADC quantification in heterogeneous gliomas.

Whereas previously proposed methods such as region growing, live wire or shape-based segmentation require and produce sharp boundaries, statistical segmentation approaches are more suitable to incorporate soft boundaries by yielding class probabilities. Gaussian mixture models with the help of expectation maximization (EM) optimization[] or similar optimization strategies are one approach to achieve a probabilistic classification based on histogram analysis[].

An important first requirement for statistical segmentation approaches is that the distributions of the measurement value (e.g. the ADC) within the expected different segments (e.g. high vs low cellularity) can be separated into a superposition of 2 Gaussians. If this is given for the tumor volume, the probability can be calculated that this voxel belongs to one of the expected clusters. In Fig. 1, we illustrate a theoretical distribution of the ADC in a tumor region in the brain.

Figure 1

Theoretical histogram: red depicts the high-proliferating area, green the low-proliferating region, brown the necrosis area and blue the CSF. The volume fractions of the different compartments were set to 0.4:0.3:0.15:0.15 for high-proliferating, low-proliferating, necrosis, and CSF, respectively.

Based on these individual probabilities, using the EM approach[], the assumed tissue groups can be found. The overall probability of the found cluster, which can be seen as a quality measurement of the complete approach, is then given by the so-called posterior probability. Furthermore, the given probabilities within the tumor area can be used to take partial volume effects into account[]. These likelihoods can be used to come to a probabilistic quantification that allows for a reproducible quantification in regions of tissue inhomogeneity[,]. The probabilistic quantification does not necessarily ensure a more valid delineation. However, the objectivity of the quantification increases by applying probabilistic assumptions.

In this article, we describe a method to objectively separate high-proliferating from low-proliferating parts in heterogeneous gliomas represented in diffusion-weighted images. We evaluate the plausibility and the reproducibility of the implemented approach that includes multimodal image overlay, signal intensity thresholding and probabilistic ADC-based clustering.

Materials and methods

Ten patients, 8 males and 2 females, with histologically proven grade III–IV glioma were included in the study. The group consisted of 2 grade III and 8 grade IV patients; the mean age was 69 ± 14 years. Patient names were rendered anonymous and the resulting code used to uniquely link data and patient. Diffusion tensor imaging (DTI), axial T2-weighted fluid attenuated inversion recovery (FLAIR), and contrast-enhanced T1-weighted magnetic resonance imaging (MRI) data were acquired and processed in order to obtain a quantification of ADC in high and low malignant tissue. The extent of the gross tumor volume (GTV) was delineated manually by a radiologist on the T2-weighted FLAIR data. T1-weighted contrast-enhanced imaging was used for delineation of the high-proliferating area within the GTV. Fig. 2 gives an overview of the steps performed during the analysis.

Figure 2

Overview of the processing steps from image acquisition to quantification. CSF, cerebrospinal fluid; CA, contrast agent; EM, expectation maximization.

Imaging parameters

Measurements were performed on a 1.5-T whole-body clinical MRI scanner equipped with a quadrature head coil (Magnetom Symphony, Siemens Healthcare, Erlangen, Germany, gradient strength = 40 mT/m).

DTI

Fifty axial slices were acquired; using a thickness of 2.5 mm, 6 gradient directions, 2 b-values (0 and 1000 s/mm2) and 10 averages. Repetition time (TR)/echo time (TE) 8400/108; voxel size 2.5 mm3; resolution 96 × 96; field of view 240 × 240 mm2.

T2-weighted FLAIR

Twenty-four axial slices; TR/TE 9000/114; voxel size 0.5 × 0.5 × 5.0 mm3; resolution 448 × 512; field of view 210.11 × 210.11 mm2; inversion time 2500; slice thickness of 5 mm.

T1-weighted contrast-enhanced imaging

Twenty-four axial slices; TR/TE 700/18; voxel size 0.5 × 0.5 × 5.0 mm3; resolution 512 × 512; field of view 240.13 × 240.13 mm2.

ADC calculation

Post-processing was performed using the MeVisLab development environment (Fraunhofer MEVIS, Bremen, Germany). The tensor components were estimated by linear least squares calculation. The ADC is determined by the following equation: where λ is the ith eigenvalue of the tensor.

Preprocessing

To improve the signal to noise ratio of the diffusion data, the 10 averages were matched spatially using a linear registration model, normalized cross-correlation as a similarity measure, and trilinear interpolation. To further improve the matching, we re-sampled the data to a finer grid of an isotropic voxel size of 1.25 mm using a cubic B-spline for the interpolation.

The axial T1-weighted contrast-enhanced images as well as the T2-weighted FLAIR images were registered onto the ADC map using normalized mutual information and trilinear interpolation.

Exclusion of CSF and necrosis

CSF is clearly hypointense in T2-weighted FLAIR images and a mask was created in the corresponding images to exclude the voxels belonging to the lower quantile in the T2-weighted FLAIR histogram. Similarly, an ADC threshold >0.001 mm2/s was used to exclude necrotic areas.

Manual placement of conservative ROIs

To obtain reference ADC values for high- and low-proliferative areas, a radiologist manually placed conservative ROIs in regions that clearly represented one of the two classes. All 3 acquired MRI datasets were taken into account.

Test for gaussianity

To justify the later use of a Gaussian mixture model for automated clustering of the 2 regions, a χ2 test was performed on the ADC distributions using ROOT (CERN, Switzerland).

Probabilistic segmentation

The probabilistic segmentation was performed by fitting a mixture model of 2 Gaussians to the normalized histogram within the masked GTV using the EM algorithm. The EM clustering with Gaussian mixture models is described in further detail in Refs.[,]. The derived posterior probabilities of the segmented voxels were used to quantify the ADC of the low- and high-proliferative regions by calculating the weighted sum by the posterior probability of the corresponding ADC measures.

Manual initialization of the model

To manually initialize the EM algorithm, the starting parameters were derived from the manually drawn ROIs. The stability of the approach with regard to the initialization was assessed by systematically choosing ADC values from the ROIs as class means. This procedure was repeated 10 times. The ADC from the conservative ROIs was averaged and the mean of the maximum values of the high-proliferating ROIs +σ and the mean of the minimum values of the low-proliferating ROIs −σ were used as a criterion to exclude runaway values in this test. To evaluate the reproducibility, dot plots were generated and intra-class coefficients (ICCs) were calculated from 10 initializations from all 10 patients using the software Stata (College Station Release 11; StataCorp LP, TX, USA).

Automated initialization of the model

In order to automatically initialize the EM algorithm, the fully automated model initialization as given in Ref.[] was implemented. The GTV describes the gross tumor volume, K is the number of clusters, μ is the mean, σ is the variance and π is the fraction of cluster k belonging to the whole Gaussian mixture model (Eq. (5)).

To evaluate the plausibility of the automated clustering, the ADC in both regions found using the automated approach was compared with the value from the seed-based initialization. The results were included in the dot plots mentioned above. Further details on the EM clustering approach are given in the Appendix.

Results

Analysis of the manually placed ROIs

The ADC distributions within the conservative ROIs with absolute frequency were averaged and showed 2 distinct peaks in all patients (see Fig. 3). The χ2 test yields the following parameters:

Figure 3

Original averaged full frequency histograms consisting of conservative areas. The high-proliferating area is highlighted in red and the low-proliferating area is shown in green. The Gaussian curve fits the area displayed (dashed curves).

High-proliferating tissue χ2 (0.95; 100−3 = 97) = 121.0 > 89.2

Low-proliferating tissue χ2 (0.95; 100−3 = 97) = 121.0 > 85.9

The following absolute ADC values were found by averaging over all patients:

High-proliferating tissue 4.13 × 10–4 ± 9.1 × 10−5 mm2/s

Low-proliferating tissue 7.26 × 10–4 ± 9.9 × 10−5 mm2/s

Histogram clustering results

The Gaussian mixture model appropriately represented the histogram of the ADC values in the GTV in all 10 patients. This is exemplarily shown on 3 patients in Figs. 4i, 5i and 6i. The class ADC values found using the clustering method corresponded well with the mean ADC values in the manually drawn ROIs (compare (f) and (g) in Figs. 4–6). Qualitatively, the found ADC distributions (purple lines) correspond well to the theoretically introduced model shown in Fig. 1.

Figure 4

Imaging data, histogram analysis, and resulting segmentation of patient 1 (WH) (a) Axial ADC map. The arrow points to a low-proliferative area of the tumor. (b) Axial contrast-enhanced T1-weighted image. The arrow indicates a high-proliferating part of the tumor. (c) Axial T2-weighted FLAIR image. The arrow points to a low-proliferative area of the tumor. (d) ADC map with GTV outlined in yellow. The arrow indicates CSF. (e) Segmentation result shows the posterior probability of high- (red) and low-proliferative (green) areas. (f) Enlarged outline of the posterior probability as shown in (e). (g) ADC histogram within the manually drawn ROIs. High- and low-proliferating distributions are represented in red and green, respectively. (h, i) ADC distribution within the GTV. The dotted line shows the ADC distribution after necrosis and CSF exclusion. The fitted Gaussian distributions for the high- and low-proliferating tissue classes are depicted in red and green, respectively. The sum of these distributions represents the data appropriately (see purple line).

Figure 5

The results for patient 2 (SH) organized as in Fig. 4. An interesting aspect in this patient is the almost complete lack of clearly definable low malignant areas. This can be appreciated by comparing the large size of the contrast-enhancing area in (b, arrow) with the hyperintense T2-weighted FLAIR area of similar size (c, arrow). The central necrotic areas of the tumor are outlined in (d, arrow and outlines). In the resulting segmentation (e), these areas are excluded due to the necrosis cut-off. The overall clustering result is plausible since mainly high malignant areas are found. Only some of the partial volume between necrosis and high-proliferative areas (e, arrow) is misclassified as low proliferating. The sum of the Gaussian distributions represents the data appropriately (i, purple line).

Figure 6

The results for patient 3 (SHE) organized as in Fig. 4. The GTV outline in this patient includes normal brain tissue (d, arrows). This leads to a misclassification of normal brain tissue as a high-proliferative area (e, arrows). The sum of the Gaussian distributions represents the data appropriately (h, purple line).

Fig. 4 demonstrates the successful use of the CSF suppression, yielding a good separation of CSF and low-proliferative areas using the T2-weighted FLAIR cut-off. Fig. 5 similarly demonstrates the usefulness of the ADC-based necrosis cut-off. Note, however, that in the border region between necrosis and high-proliferative areas, a small region of misclassification occurred (Fig. 5e, arrows). Fig. 6 again demonstrates that high- and low-proliferative areas could be separated appropriately. However, it also shows a potential pitfall in cases where the GTV includes healthy brain tissue that is erroneously classified as high malignant (Fig. 6e).

Reproducibility and automated segmentation

For all 10 patients, the clustering routine was repeated 10 times using representative initialization seed voxels (Fig. 7). Each segmentation took only 0.1 s and thus, is virtually a real-time approach. For the low- and high-proliferative area, ICCs of 0.920 and 0.996 were found, respectively, indicating good agreement[]. The automated clustering results are within the range of results found using the seed-based initialization and are indicated as gray crosses. The automated clustering was successful in all 10 cases.

Figure 7

Dot plots of the clustering results using conservative ROIs and automated clustering. (a) Results of the seed point segmentation within the high-proliferating region (black dots) and automated segmentation (black crosses). Fractional outliers were found in BM 4/10 and in JA 1/10. The ICC for the seed-based initialization was 0.92. Note that the ADC from the automated clustering is well within the range of the seed point-based measurement. (b) Results of the seed point segmentation within the low-proliferating region. No complete outliers were found. Fractional outliers were found in BM 4/10, WF 1/10. The ICC for the seed-based initialization was 0.996. Again note that the ADC from the automated clustering is well within the range of the seed point-based measurement.

Discussion

In this work, we show that low- and high-proliferating areas in high-grade gliomas can be objectively delineated and probabilistically quantified by an EM algorithm. We propose an automated and rapid approach that anticipates the heterogeneous nature of high-grade gliomas. The application of a Gaussian mixture model is validated by a χ2 test of the averaged high-proliferating and low-proliferating ROIs, yielding Gaussian distributions in the different areas.

Quantitatively, we demonstrate that the derived mean values of this theoretical model are plausible. Through the application of a statistical algorithm with a Gaussian mixture model, we are able to probabilistically quantify the ADC in areas of low- and high-proliferation in high-grade gliomas. This unsupervised algorithm is plausible because the overall distribution fits well to the raw data and the found clusters correspond well to the conservative ROIs. The high ICCs indicate a good reproducibility of the clustering approach. Furthermore, automated clustering was successful in all 10 patients.

By application of a binary T2-weighted FLAIR mask, CSF could be excluded robustly (Fig. 4). Necrotic areas could be suppressed by application of an ADC-based histogram cut-off. However, the necrosis threshold is not completely satisfactory because it cannot suppress necrosis totally without additionally cutting off apparent low-proliferating tumor tissue (Fig. 5).

Limitations

A limitation of our current approach is that the signal overlap between normal tissue and high-proliferating areas[] is not solved. Potentially, the addition of a partial volume class may mitigate this error. A further, related limitation is that the method only clusters on a predefined GTV drawn on a two-dimensional slice. Future research should extend the method onto the whole brain, resulting in a fully automated three-dimensional whole-brain clustering algorithm that includes and defines contextual three-dimensional constraints. Another limitation of our study was that no biopsies from different tumor regions could be obtained. Thus, we can only rely on the available image data and the unique ADC distribution to obtain a plausible separation of high- and low-proliferative areas.

The resulting problems of intensity overlap as discussed above and the definition of three-dimensional constraints may be addressed by incorporation of algorithms that include prior knowledge of spatial and anatomical constraints of normal anatomy and tumor tissue[,]. In these publications, a discriminative variant of the usually generative Markov random fields in combination with a support vector machines (SVM) algorithm was used, which may be combined with the approach proposed here.

A potential drawback of the EM algorithm may be data overfitting[]. Although deterministic, the algorithm may converge towards a local maximum in the log-likelihood step. This is due to singularities that can occur if one mean of the Gaussian mixture model μ equals x which denotes an observed data point. In this case, the log-likelihood term goes to infinity and thus has to be excluded from the maximization procedure and only a local maximum can be achieved. However, in our case, we did not encounter these potential difficulties regarding the objectivity of the retrieved ADC values, considering the applied test. In Ref.[] several EM techniques were evaluated for brain segmentation to overcome singularities and compared with underlying ground truth. A variational EM algorithm and 2 hybrid approaches, using genetic algorithms in combination with EM and the variational approach plus genetic algorithms, respectively, were evaluated. In Ref.[], it can be clearly seen that the simple hybrid algorithm (GA plus EM) outperforms the classic EM algorithm in terms of a higher log-likelihood result and moreover a faster convergence at equal computational costs. In Ref.[], the hybrid variational approach in conjunction with EM outperforms all other proposed algorithms. Regarding complexity pertaining to the variational approach, Ref.[] describes a little complexional overhead.

Thus, the proposed variational approach in combination with genetic algorithms may improve objective parameter estimation, also in our setting. Future research should explore this combined approach.

Future research and practical applications

The probabilistic ADC measurement in gliomas as proposed here can be used for the planning of resections and biopsies, as high-proliferative areas can be identified. Moreover, it can be applied in modern radiotherapy techniques where intensity modulation is used to adapt the local dose to the biological activity in heterogeneous tumor parts[,]. The identified low-proliferating and high-proliferating parts could be also used in a tumor growth model as described in Ref.[] to more precisely estimate the proliferation rate in different tumor areas. Moreover, this approach could be established in other heterogeneous tumor types where the ADC is used for tumor therapy monitoring.

Conclusion

In conclusion, we have presented a method for the objective delineation of high- and low-proliferating areas in high-grade gliomas. The probabilistic quantification allows for an objective and reproducible extraction of the ADC values in the different clusters. These objective values, reflecting high-proliferating and low-proliferating regions in gliomas may enable the detection of slight changes in ADC under therapy, and can, for example, aid the differential diagnosis between pseudoresponse and pseudoprogression.