Assessing the scale of tumor heterogeneity by complete hierarchical segmentation of MRI

The SHARP algorithm takes as input a 2D or 3D image of a tumor, and outputs:

• (a)  
  a binary tree that divides the tumor into regions;
• (b)  
  for each node of the tree, the difference in SI histograms between the two children of that node, as measured by the Earth mover's distance (EMD);
• (c)  
  a list of distance scales (e.g. {32 mm, 16 mm, 8 mm, 4 mm}), and the mean SI heterogeneity at each scale as measured by the mean EMD as calculated in (2)

In practice, the SHARP algorithm has two steps (visual schematic is shown in figure 1, pseudocode is listed in figure 2). First, the tumor is recursively divided in a binary fashion into regions until each voxel has been isolated. The dividing lines are chosen with the goal of maximizing the difference in SI between each pair of regions, using the Normalized Cuts algorithm (Shi and Malik 2000). This algorithm treats the input region as a graph, with one node representing each voxel. The weight between nodes i and j is:

$$\begin{eqnarray}&&{w_{ij}} = {{\rm{e}}^ - }^{{{\left( {\frac{{|{\rm{SI}}\left( i \right) - {\rm{SI}}\left( j \right)|}}{{{\sigma _I}}}} \right)}^2} - {{\left( {\frac{{{{\left\| {X\left( i \right) - X\left( j \right)} \right\|}_2}}}{{{\sigma _X}}}} \right)}^2}},\end{eqnarray} \tag{ 1 }$$

where SI(i) is the signal intensity at node i, X(i) is the spatial location of node i, and σ_I and σ_X are constants defining the sensitivity of the partitioning to SI and distance difference between voxels. The input region is divided into two partitions A and B with the goal of maximizing the weight of the edges within each partition, and minimizing the weight of the edges between partitions. An eigenvalue problem is solved so as to maximize N_cut, defined as:

$$\begin{eqnarray}&&{{N}_{\text{cut}}}(A,B)=2-{{{\sum}^{}}_{u\in A,t\in B}}w(u,t).\end{eqnarray} \tag{ 2 }$$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To reduce execution time and memory requirements, we use the Nyström approximation to the Normalized Cut, which first solves the grouping problem for a random subset of the input voxels and then extends the two groups to include all the input voxels (Fowlkes et al 2004). Because SI on MRI does not have any absolute meaning akin to Hounsfield units on CT, the SI term in equation (1) is normalized by dividing SI by the interquartile range of the whole tumor's SI (pseudocode line 1). Similarly, the distance term in equation (1) is normalized to the approximate diameter of the input region: σ_X is scaled by the diameter of the sphere with equal volume to the input region (pseudocode line 8). This step is performed because when σ_X was kept constant despite varying volume of the input region, we found that for very small input regions the distance term in equation (1) was dominated by the SI difference term. The distance normalization ensures that the same arrangement of signal intensities will produce the same partitions regardless of the scale of the input region.

The Normalized Cut algorithm does not guarantee that at each partition step, each of its two partitions is connected (i.e. a path exists from each voxel to every other voxel without leaving the region). We therefore added an evaluation for each partition step: if there are more than two connected components, the components are reassigned so that there are only two connected partitions. Each connected component other than the largest is reassigned to a component to which it is in contact (by a 4-connected neighborhood for 2D images and a 6-connected neighborhood for 3D images).

The output of the recursive partitioning is a binary tree that hierarchically divides the tumor; at the greatest tree depth each node represents a single voxel. We quantify SI heterogeneity at each partition by measuring how dissimilar the two partition regions are from each other. Dissimilarity is measured by the EMD (Rubner et al 2000) between their normalized SI histograms. If the frequency bars on a histogram plot are thought of as piles of dirt, the EMD is the quantity of dirt that must be moved, multiplied by the horizontal distance it is moved, to make the two histograms match. Mathematically it is defined as:

$$\begin{eqnarray}&&{\rm{EMD}} = \sum\nolimits_{i = 1}^n {|{\rm{CD}}{{\rm{F}}_a}\left( i \right) - {\rm{CD}}{{\rm{F}}_b}\left( i \right)|} ,\end{eqnarray} \tag{ 3 }$$

where histograms a and b each have n bins and CDF_a represents the cumulative distribution function of histogram a.

The EMD between identical histograms is 0; EMD increases as the shape of the two histograms diverges, up to a maximum of the number of histogram bins minus one. Compared to simpler metrics such as difference in mean SI, EMD has the advantage of being more sensitive to differences between regions. For example, consider two tumors A and B, each divided into regions #1 and #2. In region A₁, half the voxels have SI of 1 and half have SI of 3. In region A₂, all voxels have SI of 2. In regions B₁ and B₂, all voxels have SI of 2. For both tumors, the difference in mean SI between regions is 0. However, the EMD is higher for tumor A than tumor B: for histogram bin size of 1, EMD_A = 1 and EMD_B = 0. The EMD corresponds better to the visual appearance of the tumor regions—the two regions of tumor A would appear quite different but the two regions of tumor B would look the same.

For each parent node of the partition tree, we plot EMD against region scale in millimeters (diameter of the sphere with equal volume to the input region) in a scatterplot (figure 1(b)). Data are summarized by binning region scale on a logarithmic scale. From this summary plot, we can extract features such as scale of maximum heterogeneity and slope of heterogeneity versus scale. We can calculate other features from the partition tree, including the area of the tumor with the most/least heterogeneity at each scale, and for a given point in the tumor, the amount of heterogeneity at different scales (figure 3).

Zoom In Zoom Out Reset image size

SHARP is implemented in MATLAB R2012b (MathWorks, Natick, MA). Source code and an example dataset are available on GitHub [https://github.com/MGensheimer/SHARP]. We use the implementation of Normalized Cuts in Sameer Agarwal's freely available Spectral Clustering Toolbox (Agarwal 2002). Running time for a tumor with 4000 pixels is 48 s on a Linux system with 16 GB of memory and two Intel Xeon E5620 processors with eight cores.

We performed validation of the algorithm using three classes of synthetic images: one with more coarse-scale heterogeneity ('Large blocks'), one with a mix of coarse- and fine-scale heterogeneity ('Large + small blocks'), and one with more fine-scale heterogeneity ('Small blocks') (figure 4(a)). Image size was 64  ×  64 pixels, with pixel size of 1  ×  1 mm. Each image was initialized with a uniform value of 0. Then, many rectangles were added, with intensity in the rectangle shifted by an amount randomly sampled from a uniform distribution centered at 0. In the three classes of images, the rectangles of different size are weighted differently. For the 'Large blocks' image, the highest weighted rectangles are large, of size 32  ×  64 mm. For the 'Small blocks' image, the largest rectangles are given a very low weight and the smallest rectangles, of size 1  ×  1 mm, are given the highest weight. For the 'Large + small blocks' image, all sizes of rectangles have the same weight.

Zoom In Zoom Out Reset image size

SI was normalized to a scale of 0 to 1, with outlier pixels (intensity below 2nd percentile or above 98th percentile) set to 0 or 1. C_I and C_X were set by running the algorithm on synthetic images and one patient's MRI, varying C_I and C_X, and choosing values that led to partitions that tracked along SI gradients but reliably produced two cohesive regions. The same C_I and C_X values were used throughout this work. 300 samples were used for the Nyström approximation of the Normalized Cut. 16 histogram bins were used when calculating EMD between regions. To report mean EMD at each distance scale, parent region diameter was grouped into six bins ranging from 2 to 64 mm, equally spaced on a logarithmic scale. MATLAB code for generating and analyzing the synthetic images can be found in the synthetic.m file in the GitHub repository.

It is important that the algorithm be robust to artifacts such as noise or intensity gradients that can commonly be present in MR images. Therefore we ran simulations where Gaussian white noise or a smoothly varying intensity gradient was added to the three classes of synthetic images. Various strengths of noise or gradient, defined as a percentage of the variation in the initial image, were applied. For example, 10% Gaussian noise meant adding noise with a Gaussian distribution with standard deviation equal to 10% of the standard deviation of the initial image (σ_I). Since standard deviation is not defined for a linear intensity gradient, for this perturbation the strength of the gradient was measured in percentages of σ_I·4 with the justification that 95% of the values of the original normally distributed image lie within a range of σ_I·4. So a gradient with 100% strength would create added variation equal to the variation already present in the image.

We measured the effects of the image perturbations on the SHARP algorithm output using intraclass correlation, as recommended by the Quantitative Imaging Biomarkers Alliance (Raunig et al 2015). For each of the three synthetic image classes, and for each level of image perturbation, 20 random images were generated. Then, for each distance scale the EMD was compared between images of the same class with different amounts of noise/gradient (within class variation), and between images of the three different classes (between class variation). Intraclass correlation was calculated using the ICCest function in the ICC R package. An intraclass correlation value of >0.8 indicates good within-group agreement (Mosher et al 2011).

Using the SHARP algorithm, we examined the primary tumor of 18 patients with rhabdomyosarcoma on initial diagnosis MRI, to explore the influence of scale of SI heterogeneity on relapse-free survival (RFS). The Seattle Children's Hospital Institutional Review Board approved this retrospective study. Eligibility for the study included diagnosis with rhabdomyosarcoma between 2003–2012, pre-treatment contrast-enhanced MRI available in the digital radiology system, and treatment at Seattle Children's Hospital. All patients were treated with chemotherapy and local treatment (radiation, surgery, or both) according to national standard regimens. Only patients with primary tumor greater than or equal to 2 cm in maximal diameter were included, to ensure that several scales of heterogeneity could be analyzed with SHARP. In total, 18 patients met the eligibility requirements and were included in the analysis. We analyzed the T1-weighted gadolinium-enhanced sequence. In-plane resolution ranged from 0.31  ×  0.31 to 1.32  ×  1.32 mm. Slice spacing ranged from 3 to 8 mm. Median TR was 568 ms (range 34–850), and median TE was 12 ms (range 3–29).

The tumor border was manually drawn by M.G. on three equally spaced parallel image slices using 3D Slicer version 3 (Pieper et al 2004). Tumor SI was normalized to a scale of 0 to 1. Then, the SHARP algorithm was run in 2D on the tumor portion of each of the three slices. The lists of parent region size and corresponding EMD for the three slices were concatenated.⁴ To report mean EMD at each distance scale, parent region diameter was grouped into bins equally spaced on a logarithmic scale. We ran a separate analysis for two bin spacings: powers of 2 (bin centers {0.25 mm, 0.5 mm, 1 mm, ..., 128 mm}), and powers of $\sqrt{2}$ (bin centers {0.25 mm, 0.353 mm, 0.5 mm, 0.707 mm, ..., 128 mm}). To enable unbiased comparison of various sized tumors, distance scale was then divided by the maximal diameter of each tumor.

For calculation of RFS, events were disease relapse (growth of treated lesions or appearance of new lesions) or death. Median follow-up was 35 months. For RFS analysis, patients without an event were censored at last follow-up. Imaging covariates were compared to RFS using univariate Cox proportional hazards regression. We extracted measures expressing the scale of heterogeneity of MRI, rather than the absolute quantity of heterogeneity:

• (a)  
  MaxScale: log₂(scale with maximum heterogeneity), with heterogeneity defined as mean EMD at that scale
• (b)  
  HetSlope: Slope of least squares regression line fitting heterogeneity versus log₂(scale)

For both of these measures, scale was normalized to maximal tumor diameter in order to create comparable results between tumors of different diameters. In summary, four heterogeneity measures per patient were produced: scale of maximum heterogeneity and slope of heterogeneity versus scale, each for two scale bin sizes (powers of 2, powers of $\sqrt{2}$). These are abbreviated as MaxScale₂, HetSlope₂, $MaxScal{{e}_{\sqrt{2}}}$, and $HetSlop{{e}_{\sqrt{2}}}$. The coxph function in R version 3.0.2 was used for survival analysis (R Core Team 2013).

We compared the performance of SHARP to three existing methods: spectral analysis using the 2D discrete Fourier transform of the image, fractal analysis and texture analysis using Haralick features.

2.4.1. Spectral analysis.

2.4.2. Fractal analysis.

2.4.3. Texture features.

We created three classes of synthetic images, meant to approximate tumors with varying appearance, to test whether the SHARP algorithm correctly summarizes the scale of heterogeneity. We randomly generated 20 images of each class and ran the algorithm on each image. Results are presented in figure 4. The algorithm successfully differentiated between the different classes of images. It correctly identified the scale of maximum heterogeneity for the 'Large blocks' and 'Small blocks' cases, and captured the multiscale heterogeneity of the 'Large + small blocks' case.

The three comparison methods were also tested on the synthetic images. Findings are shown in table 1. As image features become more small-scale ('Small blocks' image), the power spectrum becomes flatter, indicating greater signal present at higher frequency levels. Fractal dimension scales inversely with the scale of the image features, which matches intuition as the patterns in the 'Small blocks' image are more space-filling than those in the 'Large blocks' image. Some of the texture results also make intuitive sense. Most adjacent pixels in the 'Large blocks' image have similar values, in contrast to the 'Small blocks' image, so the f9 entropy value increases with smaller image feature size.

Robustness of the algorithm to high frequency image noise or coil/magnetic field inhomogeneity was tested by adding Gaussian white noise or a smoothly varying intensity gradient across the synthetic images (figure 5). 20 images of each class were generated and then 5%/10% noise or 10/20% gradient were added. With increasing amounts of white noise, the amount of high-frequency / small scale heterogeneity increased slightly for the 'Large blocks' and 'Large + small blocks' images. Added noise did not affect the plot of heterogeneity versus scale for the 'Small blocks' image, likely because this image already consisted of mainly high-frequency noise. Conversely, adding an intensity gradient affected mainly the output for the Small Blocks images, presumably because this was a low-frequency/large scale perturbation being added to an image that started with mostly high frequency signal.

Zoom In Zoom Out Reset image size

To test the dependence of the algorithm results on the overall shape of the tumor, we created three tumor outlines with various complexity (figure 6). All three shapes were calibrated to have the same area, which allowed the effects of tumor shape to be tested in isolation. The three synthetic images (Large blocks, etc.) were then masked with the three tumor shapes in turn and the SHARP algorithm was run on 20 randomly generated images of each class. Figure 5 shows that the results were quite stable across tumor shapes.

Zoom In Zoom Out Reset image size

We quantitatively assessed the effects of all the image perturbations on the SHARP algorithm output using intraclass correlation. For each distance scale, this measure compared (1) the variation in EMD between each of the three original synthetic images and the corresponding perturbed image (within-class comparison) and (2) the variation in EMD between synthetic images of the three different classes (between class comparison). The intraclass correlation values for each of the six distance scales were averaged to create a summary measure (table 2). For all amounts of added noise or bias field, the intraclass correlation was >0.8, which indicates good robustness of the results to the perturbations. Similarly, the intraclass correlation between the three tumor shapes was good at 0.89.

We applied the SHARP algorithm to examine the influence of scale of rhabdomyosarcoma tumor heterogeneity on recurrence patterns. Patient characteristics are listed in table 3. Median tumor diameter was 39 mm (range 21–120). The primary tumors tended to harbor more coarse-scale than fine-scale heterogeneity on contrast-enhanced MRI: scale of maximum heterogeneity was most commonly the full tumor diameter (n = 11 using powers of 2 bin spacing, n = 7 using powers of $\sqrt{2}$ bin spacing). There were 10 disease recurrences, and no deaths without prior disease progression. Using T1-weighted gadolinium-enhanced MRI, the SHARP algorithm produced four heterogeneity measures per patient, which were compared to RFS using univariate Cox proportional hazards analysis. One of the four measures, $MaxScal{{e}_{\sqrt{2}}}$, had a statistically significant association with RFS: coarser scale of maximum SI heterogeneity was predictive of shorter RFS (Cox model log-rank p = 0.05, concordance index 0.71). As seen in figure 7, patients with $MaxScal{{e}_{\sqrt{2}}}$ > median (n = 7) had median RFS of 15.5 months, versus not reached when $MaxScal{{e}_{\sqrt{2}}}$ ≤ median (n = 11), log-rank p = 0.02. Figure 8 shows examples of tumors with varying scales of heterogeneity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The other three heterogeneity measures were not predictive of RFS, with Cox model log-rank p = 0.24–0.45. Similarly, Children's Oncology Group (COG) risk group was not predictive of RFS (p = 0.35) (Malempati and Hawkins 2012). Tumor diameter was not predictive of RFS (p = 0.19).

In order to compare our SHARP method to existing spatial analysis methods, we also performed spectral, fractal, and texture analyses. Slope of the Fourier power spectrum was not predictive of RFS (Cox model p = 0.60). Neither fractal dimension nor fractal fit was predictive of RFS (p = 0.34, p = 0.36, respectively). Similarly, the three Haralick texture features were not predictive (p = 0.09–0.90).