# A three-dimensional computational analysis of magnetic resonance images characterizes the biological aggressiveness in malignant brain tumours

## Abstract

Glioblastoma (GBM) is the most frequent and aggressive type of primary brain tumour. The development of image-based biomarkers from magnetic resonance images (MRIs) has been a topic of recent interest. GBMs on pre-treatment post-contrast T1-weighted (w) MRIs often appear as rim-shaped regions. In this research, we wanted to define rim-shape complexity (RSC) descriptors and study their value as indicators of the tumour’s biological aggressiveness. We constructed a set of widths characterizing the rim-shaped contrast-enhancing areas in T1w MRIs, defined measures of the RSC and computed them for 311 GBM patients. Survival analysis, correlations and sensitivity studies were performed to assess the prognostic value of the measurements. All measures obtained from the histograms were found to depend on the class width to some extent. Several measures (FWHM and *β*_{R}) had high prognostic value. Some histogram-independent measures were predictors of survival: maximum rim width, mean rim width and spherically averaged rim width. The later quantity allowed patients to be classified into subgroups with different rates of survival (mean difference 6.28 months, *p* = 0.006). In conclusion, some of the morphological quantifiers obtained from pre-treatment T1w MRIs provided information on the biological aggressiveness of GBMs. The results can be used to define prognostic measurements of clinical applicability.

### 1. Introduction

Glioblastoma (GBM) is the most lethal and frequent primary malignant brain tumour. Patients diagnosed with the disease and treated with state-of-the-art therapies (maximal safe surgery, radiotherapy and chemotherapy) have a median survival of only around 14 months [1]. In the last 20 years, improvements in GBM treatments have been only minor.

Magnetic resonance images (MRIs) are routinely used in GBM patients for diagnosis, treatment planning and follow-up. Moreover, technological improvements currently allow images of millimetre resolution to be obtained, allowing for three-dimensional tumour reconstruction with excellent definition of macroscopic lesions [2].

GBM is one of the most proangiogenic tumour types, which leads to a complex unstable vessel network. Some MRI modalities, specifically the post-contrast T1-weighted (T1w) images, rely on this to render an image of the area where an injected contrast agent is released due to the leakiness of the vessels. Thus, tumour areas appear as hyperintense areas over the normal tissue background. Those images provide an estimate of the areas where there is a macroscopically large number of tumour cells, which is also the target for surgery. Some examples of the macroscopic part of a GBM as seen in post-contrast T1w images are shown in figure 1. The inner part of the tumour is typically occupied by necrotic (dead) tissue that has been damaged by the action of the tumour itself. Thus, the patterns observed in figure 1 are characteristic of GBMs and display irregular contrast-enhancing (CE) rim-like structures surrounding a necrotic core. Typically, T1w images are analysed together with other MRI sequences that provide complementary information. For example, T2/FLAIR sequences provide an estimate of the areas of oedema and infiltrative tumour while diffusion images may provide estimates of cellularity. However, T1w images provide the least ambiguous information and are routinely obtained in higher resolution allowing for the best tumour delineation and reconstruction.

Recent studies have proposed to correlate the tumour’s biological aggressiveness with quantitative parameters obtained from different magnetic resonance sequences [3,4] and more specifically in glioblastoma [5]. To cite some examples, Swanson *et al.* [6] used both the tumour volumes in T1w and T2/FLAIR MRI sequences at two time points to estimate the growth speed of GBMs. Several authors have studied the prognostic value of the volume of CE areas [7,8], ratios between CE tumour volume and total volume [9] or tumour diameter [9,10]. However, many of these studies were based on low spatial resolution images. Finally, Wangaryattawanich *et al.* [11] identified the CE volume and the maximal diameter as prognostic parameters using volumetric images, which has not been confirmed in other studies [12].

Simple shape parameters of the T1w areas characterizing tumour shape irregularity have been found to have a clear prognostic value [13]. Also, a mathematical study using a modified Fisher–Kolmogorov equation [14] predicted the potential value of the width of CE areas as a biomarker of biological aggressiveness. An indirect confirmation of the prediction was reported in [12]. Several studies have used machine-learning approaches to construct prognosis models from sets of descriptors of tumour shape, texture, etc., obtained from different MRI sequences [15–20]. Finally, other authors have explored different types of spatial imaging biomarkers based on the quantification of tumour subregions [21–23].

However, the complexity of the CE areas in T1w MRIs has not been studied in detail and correlated with a tumour biological aggressiveness.

In this study, we aimed to construct geometric measurements allowing the complexity of CE regions in volumetric (i.e. high-resolution) post-contrast T1w MRIs to be quantified computationally and the values of those measurements correlated with surrogates of biological aggressiveness, i.e. patient survival. We did this in a very large patient sample of 311 patients, indeed, larger than any dataset used in previously published clinical studies.

### 2. Material and methods

#### 2.1. Tumour segmentation and three-dimensional reconstruction

Active tumour areas are contained between two surfaces. The outer, *Σ*_{E} is the interface between the ‘normal’ tissue and CE areas. The inner, *Σ*_{I}, is the interface between the tumour tissue and the necrotic areas. In practice, these surfaces are not known explicitly as images are discretized into voxels and the algorithms have to work on the coordinates obtained from the tumour segmentation from the MRIs. In our case, The Digital Imaging and Communications in Medicine files were imported into the scientific software package Matlab and were preprocessed by using an in-house semi-automatic segmentation procedure. Tumours were automatically delineated by using a grey-level threshold chosen to identify the largest CE tumoural volume, described [13]. Then, the image expert manually corrected each segmentation section by section. The segmentations were supervised by experienced radiologists and the methodology found to be reproducible.

Let ${\mathcal{P}}_{I,E}$ be the discrete set of points *x*, *y*, *z* ∈ *Σ*_{I,E} defining the discrete approximations to the surfaces *Σ*_{I,E}. The segmentation of the tumours as shown in figure 2 provides the coordinates of points in ${\mathcal{P}}_{I,E}$ and allows for the reconstruction of the surfaces.

#### 2.2. Quantification of the contrast-enhancing ‘rim’ morphology

Quantifying the space contained between surfaces is a topic of interest in many scientific fields [24–28]. Typically, the goal is to measure the deviation of one surface with respect to the other to estimate morphological changes. Our problem is different since we are interested in quantifying the rim complex structure. In what follows, we will denote by $\mathcal{R}\subset {\mathbb{R}}^{3}$ the space contained between the surfaces *Σ*_{E} and *Σ*_{I}.

To describe $\mathcal{R}$, we first constructed the set of distances between the *N*_{E} points of ${\mathcal{P}}_{E}$ and those *N*_{I} points in ${\mathcal{P}}_{I}$. To do so, we first computed the distances from every point in ${\mathcal{P}}_{E}$ to those in ${\mathcal{P}}_{I}$. For every point in ${P}_{E}^{k}\in {\mathcal{P}}_{E}$, we looked for the point ${P}_{I}^{q}$ providing the minimal distance of all points in ${\mathcal{P}}_{I}$ and defined ${\delta}_{E}^{k}={P}_{E}^{k}-{P}_{I}^{q}$. Next, we constructed the set

Figure 3 shows, through a colour map plotted on the surface *Σ*_{E}, the set of distances Δ_{E} associated with each point.

Next, the process was repeated for all the points in ${\mathcal{P}}_{I}$. For each of the *N*_{I} points in ${\mathcal{P}}_{I}={\{{P}_{I}^{q}\}}_{q=1}^{{N}_{I}}$, we obtained the minimal distance vector to ${\mathcal{P}}_{E}$ $({\delta}_{I}^{q})$. In addition, we required that this vector did not cross the surface *Σ*_{I}. Since we did not have an explicit parametrized form for *Σ*_{I}, this condition was implemented by requiring that the vector ${\delta}_{I}^{q}$ be sufficiently far from all points in ${\mathcal{P}}_{I}$. To make this more precise, the angles *θ*^{q,j} spanned by the vector ${\delta}_{I}^{q}$ and the one generated by the point ${P}_{I}^{q}$ and the other points ${P}_{I}^{j}$ were constrained to be larger than a small number *ε* (taken here to be 0.035). In other words, we define

Figure 4 shows a schematic case (in two-dimension) in which the crossing condition is not satisfied.

Finally, we define the final set of distances Δ characterizing $\mathcal{R}$:

*Σ*

_{E,I}and one of the minimal distances between them included in the set Δ.

Our technique has some similarities to the one described in [24], but with substantial modifications to make it more suited to our problem. The improvements imply that: all points on both surfaces are analysed, segments from *Σ*_{I} to *Σ*_{E} do not cross *Σ*_{I}, segments that are obtained twice are counted only once and new measurements are constructed characterizing $\mathcal{R}$.

#### 2.3. RSC distance distribution

From the set Δ, we constructed histograms for the distribution of distances for different class widths *t*. Those histograms Δ_{t} will be used to compute some of the measurements of the study to be defined later.

The choice of the distribution class width was made taking into account the minimal discrete spatial scales characterizing the image resolution: i.e. the MRI pixel resolution, the slice thickness and the spacing between slices (named according to the DICOM standard).

We chose *t* = 0.5 mm and *t* = 0.75 mm, as class widths for our study, because: (i) 0.5 is the best MRI pixel spacing resolution and (ii) 0.75 is the intermediate value between the best and the worst resolution (0.5 mm and 1 mm, respectively). We explored also 1 mm as class widths but there were very few classes grouping together many different data. The choice of two values for *t* allowed us to perform a sensitivity analysis to check whether the parameters derived from the distributions Δ_{0.5} and Δ_{0.75} were robust under changes in the class width.

#### 2.4. Centralization measures of RSC distances

*Mean distance*(*δ*_{mean}). Defined as the mean of the set of distances in Δ, i.e. ${\delta}_{\mathrm{mean}}=\sum _{n=1}^{N}{\delta}_{n}/N$.*Maximum distance*(*δ*_{max}). The largest measurement in the set*δ*_{max}= max_{j=1, … ,N}*δ*_{j}. It corresponds to the largest rim width.*Modal distance*(*δ*_{mode,t}). The mode of the distribution Δ_{t}. This measurement is dependent on the class width*t*.*Spherical rim width*(*δ*_{s}). Let us first define the total tumour volume (*V*_{T}), and the necrotic (or inner) volume (*V*_{I}) as the volumes enclosed by the surfaces*Σ*_{E}and*Σ*_{I}, respectively. The CE volume is then defined as*V*_{CE}=*V*_{T}−*V*_{I}. These volumes can be easily obtained from the tumour segmentations and the voxel dimensions. Then, the spherical rim width is a rough approximation to the averaged rim width computed by the equation$${\delta}_{\mathrm{s}}={\left(\frac{3}{4\pi}\right)}^{1/3}[{({V}_{\mathrm{CE}}+{V}_{I})}^{1/3}-{{V}_{I}}^{1/3}].$$2.4This quantity measures the average width of the contrast-enhancing areas by assuming sphericity of the necrotic volume and also that CE areas are placed peripherally with an annular shape. This measurement can be computed directly from the tumour volumes, thus its computation does not require the construction of the set Δ.

#### 2.5. Geometric irregularity (dispersion) measurements

*Full width at half maximum*(*FWHM*). This is a measurement frequently used in different scientific fields to characterize the ‘width’ of a distribution [29–31]. It is given by the distance between points on the distribution at which the function reaches half its maximum value. To compute it for the distribution Δ_{t}, we take the frequency polynomial*f*and the distances*δ*_{l}and*δ*_{r}closest to*δ*_{mode}, on the left- and right-hand sides, respectively, such that $f({\delta}_{\mathrm{l}})=f({\delta}_{\mathrm{r}})=\frac{1}{2}f({\delta}_{\mathrm{mode}})$, then$$\text{FWHM}={\delta}_{\mathrm{r}}-{\delta}_{\mathrm{l}}.$$2.5*a**FWHM*_{N}*and FWHM*_{M}are two measurements of FWHM normalized to the maximal distance*δ*_{max}and to the full range of distances within the tumour*δ*_{max}−*δ*_{min}, respectively, i.e.$${\text{FWHM}}_{N}=\frac{\text{FWHM}}{{\delta}_{max}}$$2.5and*b*$${\text{FWHM}}_{M}=\frac{\text{FWHM}}{{\delta}_{max}-{\delta}_{min}}.$$2.5c*Standard modal deviation*(*D*_{mode}). This measurement quantifies the deviation of all the distances with respect to*δ*_{mode}. It is given by the equation$${D}_{\mathrm{mode}}=\sqrt{\frac{1}{n}\sum _{i}^{N}{({\delta}_{i}-{\delta}_{\mathrm{mode}})}^{2}}.$$2.6*Mode variation coefficient*(CV_{mode}). This relates the value of*δ*_{mode}and its dispersion*D*_{mode}. Larger values of this quantity correspond to a larger dispersion in the distribution. This is a dimensionless variable and it is typically expressed as a percentage as given by the equation$${\text{CV}}_{\mathrm{mode}}(\text{\%})=\frac{{D}_{\mathrm{mode}}}{{\delta}_{\mathrm{mode}}}\times 100.$$2.7*Standard mean deviation*(*D*_{mean}). This quantifies the amount of variation of all distances with respect to*δ*_{mean}:$${D}_{\mathrm{mean}}=\sqrt{\frac{1}{n}\sum _{i}^{N}{({\delta}_{i}-{\delta}_{\mathrm{mean}})}^{2}}.$$2.8This quantity is computed directly on Δ and thus does not depend on the class size*t*.*Mean variation coefficient*(CV_{mean}). This quantity relates the value of*δ*_{mean}and its dispersion*D*_{mean}as given by the equation$${\text{CV}}_{\mathrm{mean}}(\text{\%})=\frac{{D}_{\mathrm{mean}}}{{\delta}_{\mathrm{mean}}}\times 100.$$2.9*Geometrical heterogeneity*(*G*_{h}y*G*_{H}). These two measurements quantify the ratio of the dispersion of the largest measurements over the whole set of those in Δ. Their mathematical expressions are$${G}_{H}=\frac{{\delta}_{max}-{\delta}_{Q3}}{{\delta}_{max}}$$2.10and*a*$${G}_{h}=\frac{{\delta}_{max}-{\delta}_{Q3}}{{\delta}_{max}-{\delta}_{min}}.$$2.10*b**Width between mode and mean*(WMM). Is the distance between*δ*_{mode}and*δ*_{mean},$$\text{WMM}={\delta}_{\mathrm{mean}}-{\delta}_{\mathrm{mode}}.$$2.11*a*From this measurement, related normalized quantities can be defined as

$${\text{WMM}}_{N}=\frac{\text{WMM}}{{\delta}_{max}}$$2.11and*b*$${\text{WMM}}_{M}=\frac{\text{WMM}}{{\delta}_{max}-{\delta}_{min}},$$2.11*c*

#### 2.6. Decay measures

We were also interested in the decay rate of the distributions for large distances as estimates of the dispersion of distances in the tumours. We defined several measurements of decay computed from the mode *δ*_{mode}. The interest of those measurements lies in the fact that high frequencies in classes above *δ*_{mode} may identify tumours with anisotropic growth and/or the presence of regions with different biological aggressiveness and infiltration.

*Exponential decay*(*β*_{E}). We fitted the parameters*K*and*β*_{E}to the model$${\beta}_{E,t}\hspace{0.17em}:\hspace{0.17em}{\mathrm{\Delta}}_{t}(\delta )\approx K{\text{e}}^{{\beta}_{E}\delta}\phantom{\rule{1em}{0ex}}\mathrm{\forall}\delta \in {\mathrm{\Delta}}_{t},\delta \ge {\delta}_{\mathrm{mode}},$$2.12using the values of frequencies larger than the modal class for Δ_{0.5}and Δ_{0.75}. The strength of decay is given by the exponent*β*_{E}.*Rational decay*(*β*_{R}). This is given by the model$${\beta}_{R,t}\hspace{0.17em}:\hspace{0.17em}{\mathrm{\Delta}}_{t}(\delta )\approx K{\delta}^{-{\beta}_{\mathrm{R}}}\phantom{\rule{1em}{0ex}}\mathrm{\forall}\delta \in {\mathrm{\Delta}}_{t},\delta \ge {\delta}_{\mathrm{mode}},$$2.13where the parameters*K*and*β*_{R}in equation (2.13) are fitting parameters. The decay rate is given by the exponent*β*_{R}.*Half width at half maximum*(HWHM). This is defined as the distance between*δ*_{r}and*δ*_{mode}$$\text{HWHM}={\delta}_{r}-{\delta}_{\mathrm{mode}}.$$2.14Normalized versions of this quantity by*a**δ*_{max}and*δ*_{min}are given by the following equations$${\text{HWHM}}_{N}=\frac{\text{HWHM}}{{\delta}_{max}},$$2.14and*b*$${\text{HWHM}}_{M}=\frac{\text{HWHM}}{({\delta}_{max}-{\delta}_{min})}.$$2.14*c*

Figure 6 shows examples of the distribution Δ_{0.5} for a particular patient and the location of some of the measurements described above.

#### 2.7. Patients

This study was based on data from the observational studies GLIOMAT and TOG. Both studies were first approved by the Institutional Review Board (IRB) of the Hospital General Universitario de Ciudad Real (approval numbers 03/2014 and 10/2016). All patients data were analysed anonymously and according to the principles expressed in the Declaration of Helsinki. Informed consent was obtained from all patients participating in the TOG study. For the GLIOMAT study, informed consent was waived as patients had either previously provided authorization for use of their medical records for research or were deceased.

A total of 311 patients from 10 local medical institutions diagnosed with GBM in the period 2006–2017 were included in the study. Inclusion criteria were: pathologically proven GBM according to the 2007 World Health Organization (WHO) Classification of Tumours of the Central Nervous System [32], availability of a pre-treatment volumetric contrast-enhanced (CE) T1w MRI sequence (slice thickness ≤ 2.00 mm, pixel spacing 1.20 mm, no gap), no substantial imaging artefacts, presence of contrast-enhancing areas and unifocal lesions. These sequences were gradient echo using three-dimensional spoiled gradient recalled echo or three-dimensional fast field echo. Availability of the relevant clinical variables was also required: age, treatment scheme and survival information. Table 1 summarizes patient data.

patients | |

number of patients (censored) | 311 (27) |

age (years) median [range] | 63 [19 − 86] |

sex (M, F) | $(44\text{\%},56\text{\%})$ |

survival (months), median [range] | 12.76, [0.13, 82.97] |

type of resection | |

total | $149\hspace{0.17em}(47.91\text{\%})$ |

subtotal | $113\hspace{0.17em}(36.33\text{\%})$ |

biopsies | $49\hspace{0.17em}(15.76\text{\%})$ |

type of treatment | |

chemotherapy (CT) | $5\hspace{0.17em}(1.61\text{\%})$ |

radiotherapy (RT) | $27\hspace{0.17em}(8.68\text{\%})$ |

CT and RT | $241\hspace{0.17em}(77.49\text{\%})$ |

no treatment | $38\hspace{0.17em}(12.22\text{\%})$ |

MRI characteristics | |

pixel spacing (mm), mean [range] | 0.91 [0.46–1.09] |

slice thickness (mm), mean [range] | 1.48 [1.00–2.00] |

spacing between slices (mm), mean [range] | 1.10 [0.50–2.00] |

number of slices mean [range] | 174 [80–360] |

field strength (n) |
1.5T (278), 3T (33) |

TR ms mean [range] | 13.78 [4.88–30] |

TE ms mean [range] | 3.65 [1.48–9.64] |

flip angle degrees mean [range] | 20.83 [8–45] |

manufacturer (n) |
GE (136), Philips (108), Siemens (67) |

#### 2.8. Statistical methods

Kaplan–Meier analysis was used to identify individual variables associated with prognosis, using the log-rank test to assess the significance of the results. Kaplan–Meier curves measure the fraction of patients living for a certain amount of time [33,34]. A two-tailed significance level of *p*-value lower than 0.05 was applied. For each variable, we searched for every threshold value splitting the sample into two different subgroups. Then, we chose as best the non-isolated significant value giving the lowest log-rank *p*. Univariate Cox proportional hazards regression analysis [35] was used to obtain the hazard ratio (HR) and its adjusted 95% confidence interval (CI) for each threshold.

The Spearman correlation coefficient was used to assess the dependencies between each pair of variables, in order to identify measurements with similar information [36,37]. Normality of the variables was assessed by the Kolmogorov–Smirnov test. Significant correlation coefficients over 0.75 or below −0.75 were regarded as strong correlations between variables.

To quantify the measurement’s sensitivity to changes in the distribution class width, we used Wilcoxon and Lin tests. The Wilcoxon signed-rank test is a non-parametric statistical hypothesis test used to compare two related samples, or repeated measurements on a single sample to assess whether their population mean ranks differ. The concordance Lin test quantifies data similarity [38,39]. It is considered that Lin coefficient values over 0.99 correspond to perfect correlations. Values above 0.95 (but below 0.99) are considered to denote substantial correlations. Values below 0.95 are considered to be indicators of poor concordance.

SPSS software (v. 22.0.00) was used for the statistical analysis.

### 3. Results

Firstly, the 18 measurements defined in the Material and methods section were computed for the set of 311 patients. Next, we analysed the prognostic value of those measurements by using the Kaplan–Meier analysis with threshold search methods as described in the ‘statistical methods’ section. Results are summarized in table 2 and figures 7 and 8.

variable | units | t |
threshold | p |
HR (CI95%) | dif. surv. |
---|---|---|---|---|---|---|

tumour size measurements | ||||||

V_{T} |
cm^{3} |
n.a. | 28.55 | 0.107 | n.s. | n.s. |

V_{CE} |
cm^{3} |
n.a. | 17.12 | 0.059 | n.s. | n.s. |

V_{I} |
cm^{3} |
n.a. | 5.71 | 0.223 | n.s. | n.s. |

d_{max}3D |
cm | n.a. | 5.15 | 0.068 | n.s. | n.s. |

RSC measurements | ||||||

δ_{max} |
cm | n.a. | 1.02280 | 0.040 |
1.280 (1.010–1.621) | 4.08 |

δ_{mean} |
cm | n.a. | 0.30063 | 0.033 |
1.326 (1.023–1.718) | 2.83 |

δ_{s} |
cm | n.a. | 0.37081 | 0.006 |
1.503 (1.122–2.016) | 6.28 |

δ_{mode} |
cm | 0.75 | 0.33750 | 0.414 | n.s. | n.s. |

δ_{mode} |
cm | 0.50 | 0.3750 | 0.271 | n.s. | n.s. |

dispersion measures | ||||||

G_{H} |
n.a. | 0.61839 | 0.099 | n.s. | n.s. | |

G_{h} |
n.a. | 0.61667 | 0.116 | n.s. | n.s. | |

D_{mean} |
cm | n.a. | 0.12502 | 0.154 | n.s. | n.s. |

CV_{mean} |
% | n.a. | 41.5340 | 0.099 | n.s. | n.s. |

D_{mode} |
cm | 0.75 | 0.21986 | 0.180 | n.s. | n.s. |

cm | 0.50 | 0.12930 | 0.137 | n.s. | n.s. | |

CV_{mode} |
% | 0.75 | 45.1020 | 0.251 | n.s. | n.s. |

% | 0.50 | 57.3280 | 0.270 | n.s. | n.s. | |

FWHM_{N} |
0.75 | 0.14422 | 0.072 | n.s. | n.s. | |

0.50 | 0.24548 | 0.068 | n.s. | n.s. | ||

FWHM_{M} |
0.75 | 0.34129 | 0.030 |
1.371 (1.029–1.826) | 2.07 | |

0.50 | 0.24696 | 0.032 |
1.305 (1.023–1.663) | 2.07 | ||

WMM_{N} |
0.75 | 0.11711 | 0.044 |
0.730 (0.536–0.992) | 5.40 | |

0.50 | 0.10573 | 0.045 |
1.375 (1.005–1.882) | 3.72 | ||

WMM_{M} |
0.75 | 0.11946 | 0.097 | n.s. | n.s. | |

0.50 | 0.12936 | 0.043 |
0.715 (0.516–0.991) | 5.40 | ||

frequency distribution decay measures | ||||||

β_{E} |
cm^{−1} |
0.75 | 2.72520 | 0.186 | n.s. | n.s. |

cm^{−1} |
0.50 | 3.03910 | 0.109 | n.s. | n.s. | |

β_{R} |
0.75 | 0.02764 | 0.024 |
1.460 (1.052–2.027) | 6.38 | |

0.50 | 0.01450 | 0.025 |
1.389 (1.042–1.852) | 3.94 | ||

HWHM_{N} |
0.75 | 0.21315 | 0.053 | n.s. | n.s. | |

0.50 | 0.17761 | 0.028 |
1.444 (1.039–2.007) | 3.50 | ||

HWHM_{M} |
0.75 | 0.23678 | 0.022 |
1.470 (1.054–2.051) | 6.80 | |

0.50 | 0.18333 | 0.053 | n.s. | n.s. |

First for the dispersion measurements, we obtained significant results for the variable FWHM_{M}. The threshold found FWHM_{M,th} = 0.34129 for a class width *t* = 0.75 mm split patients into two groups with a median survival difference of 2.07 months (*p* = 0.030 and HR = 1.371). The same variable was also statistically significant for class width *t* = 0.5 mm (FWHM_{M,th} = 0.24696, *p* = 0.032, HR = 1.305 and same survival difference). The corresponding Kaplan–Meier plots are shown in figure 7*a*,*b*. The variable WMM_{N} was also statistically significant but its prognostic value was inverted when changing the class width (HR = 0.730 for *t* = 0.75 mm and HR = 1.375 for *t* = 0.5 mm). Finally, WMM_{M} was significant only for *t* = 0.5 mm.

Statistically significant results were found for the variable *β*_{R} both for *t* = 0.75 mm (threshold *β*_{R} = 0.02764, *p* = 0.024, HR = 1.460 and median difference in survival 6.38 months) and *t* = 0.5 mm (threshold *β*_{R} = 0.014502, *p* = 0.025, HR = 1.389 median difference in survival 3.94 months). The respective Kaplan–Meier curves are shown in figure 7*c*,*d*.

Both HWHM_{N} and HWHM_{M} were statistically significant for one of the class widths but became only marginally significant when changing to the other class used for the sensitivity analysis (table 2).

Finally, we obtained significant results for several non-histogram-based measurements. Firstly, *δ*_{max} (with threshold 1.0228 cm) was found to have *p* = 0.040, HR = 1.280 and a difference in median survival of 4.08 months (figure 8*a*). *δ*_{mean} obtained a *p* = 0.033, HR = 1.326 and a difference in survival of 2.83 months (variable threshold was 0.30063 cm). Kaplan–Meier curves for *δ*_{mean} are shown in figure 8*b*. Finally, we would like to highlight the results obtained for *δ*_{s}. Taking *δ*_{s} = 0.37081 cm as the threshold for this variable gave *p* = 0.006, HR = 1.503 and a difference in median survival of 6.28 months, this being the smallest *p*-value achieved throughout the study. Kaplan–Meier plots are shown in figure 8*c*.

Other measurements did not reach statistical significance (table 2).

Next, we performed a sensitivity analysis to detect those variables more dependent on the class width choice. Firstly, we compared the median of the values obtained for each measurement. We obtained *p* = 0.049 for *β*_{E} when changing the class width from *t* = 0.5 cm to 0.75 cm. For the other variables, we found *p* < 0.001. Thus, all class-dependent variables were sensitive to the choice of the class width. Next, we computed the concordance for all class-dependent variables. Only the variable *D*_{mode} showed substantial concordance, *ρ*_{L} = 0.9736.

We computed the Spearman correlations between all of the study variables to identify potential dependences between them. The analysis was also performed for different class widths for class width-dependent variables. It is remarkable that *β*_{R} was independent of all the other study variables. The spherical rim width *δ*_{s} was correlated only with *δ*_{mean} with a correlation coefficient *ρ*_{L} = 0.84. Figure 9 summarizes the results of the correlation analysis.

### 4. Discussion

The goal of this study was to develop a computational procedure to fully characterize CE ‘rim-like’ tumour areas in glioblastoma. We wanted to construct a set of non-trivial measurements characterizing the rim width and correlate them with patient outcome, i.e. to find their potential utility as surrogate markers of a tumour biological aggressiveness.

Several studies have identified simple geometric variables (either volumetric or longitudinal) obtained from pre-treatment MRIs having prognostic value. Among them we would like to highlight: the ratio between volumes obtained in T1 and T2/FLAIR MRI sequences [6], the CE volume [7,8], the ratio between total volume and CE volume [9] or the tumour diameter [9,10]. However, these studies had some limitations: some of them have reported contradictory results [12], and most were based on limited-resolution images, used databases with limited numbers of patients or segmented tumours with methodologies of limited accuracy.

CE rim region correlates with the most active/proliferative areas and where the blood vessels are not-functional. Angiogenesis takes place specially here and the contrast injected (gadolinium) spreads out of the vessels. A previous theoretical study suggested that CE shape contained tumour biology information and its rim morphology mainly depended on the tumour cell’s infiltration speed, proliferation rate and the tumour cell loss factor rate [14]. The numerical simulations based on that mathematical model showed that tumours with largest rim width grew faster. Following this idea, those tumours showing large CE rim width could be those with fast proliferation rate and low internal death rate. Then, the morphological heterogeneity of CE regions may reflect the biological heterogeneity within the tumour. However, no previous study has tried to characterize the complex spatio-temporal geometry of CE areas of the tumour using three-dimensional measurements.

This kind of study requires high-resolution images as recommended for clinical trials [2], to obtain reliable values of the measurements. Other authors have performed volumetric studies on high-resolution images [11,12], but none of them have studied the RSC in full.

From the methodological point of view, the first contribution of our study was the definition of the set of distances Δ as a way of characterizing the RSC and the construction of new measures based on it. Our method had some similarities with the method of Kim *et al.* [24]. However, we analysed all points in both surfaces, excluded distances corresponding to vectors intersecting the inner surface *Σ*_{I} and removed one instance of vectors counted twice.

Some of the measures included in this study are used in other areas to characterize histograms [24–28]. For these measures, the class width is a key parameter substantially influencing the results. This is why we performed a sensitivity analysis.

Special care was given to obtaining and processing patient data to ensure results of high quality. Firstly, all patients included in the study had volumetric spatial resolutions allowing faithful reconstruction of tumour geometries in three-dimension. All the tumours were segmented following a time-consuming semi-automatic procedure. Software allowing tumours to be corrected on a tablet computer, using a digital pencil, was developed to ensure the maximum accuracy of the tumour outline delineation. The number of patients was the largest used for this type of study by a considerable margin. The combination of all these factors provided both a set of state-of-the-art measurements and data to perform the statistical studies.

Of all the RSC heterogeneity measurements, the variable FWHM_{M} had a substantial prognostic value that was preserved when changing the histogram class width.

The variable *β*_{R} is a measure of how the frequencies of the widths decay as their values become larger. Small values of this variable correspond to CE areas with values located close to the median. Large values of this variable correspond to rim-like structures with a broad distribution of widths, specifically, a substantial contribution of large widths. Thus, it is a measure of rim heterogeneity. This decay rate was the outstanding measurement of the group of variables analysing frequency decay for large rim widths. It was found to be uncorrelated to all of the other study variables and was robust when changing the class width from 0.5 to 0.75 mm: 1.460 vs 1.389 for HR and *p* = 0.024 vs 0.025. This is one of the most significant results of the study, since this non-trivial measurement requires the construction of both the set Δ and the histograms Δ_{t} and thus provides an example of the applicability of the methodology described in the paper.

Other non-histogram-based measurements were found to be significant in our study, all of them in agreement with what would be intuitively expected. First, *δ*_{max} is the maximal infiltration distance found in the tumour. This measurement may provide a hint of the presence of very aggressive tumour areas. Its significance supports the hypothesis that the local width may correlate with the local tumour growth speed [14]. Also *δ*_{mean}, found to be independent of the previous measurement, was a predictor of survival. In the same way, this measurement provides an estimate of mean tumour growth speed. We found it to be correlated with survival, and thus tumour aggressiveness.

Finally, *δ*_{s} is a variable that is based only on tumour volumetric measurements, thus, it does not require the construction of the set Δ. It was surprising that this ‘simple’ quantity, based on the assumption that both the whole tumour and its necrotic part are spherical, was found to be statistically significant in that study. In our study, we confirm the strong predictive value of that measurement. Interestingly, this measurement was correlated with the more meaningful one *δ*_{mean}, which could provide a justification of its prognostic value.

Although our study focused on glioblastoma, our methodology should be applicable to other cancer types. Brain metastases have a similar structure on T1w MRIs, with a necrotic core surrounded by an active area when tumours are sufficiently large. Other types of tumour (e.g. breast) also develop central necrotic areas.

Regarding the limitations of this study, the IDH1 and other molecular markers were not available as these patients were classified according to WHO 2007. On the other hand, to collect the largest patients sample it was necessary to collect images from several hospitals (multicenter study), which may introduce minor heterogeneity in *z* axes.

It would be interesting in future studies to collect post-surgery images and clinical information to quantify the residual tumour volume after surgery, in addition to measures shown in this study, and to analyse its effect on survival patients, and combine the genomic information with these novel biomarkers to development individual prognosis model.

### 5. Conclusion

In this paper, we have described a computational methodology to characterize the macroscopic tumour areas in GBM as observed in T1w MRIs, typically known as ‘rim’-like. Tumour necrotic and CE areas were used to construct a set of distances Δ. A set of 18 measurements characterizing the three-dimensional rim width was defined and computed for the 311 patients included in the study.

Several variables were found to correlate with tumour biological aggressiveness, including: the maximum macroscopic infiltration distance within the tumour (*δ*_{max}), the rational decay exponent *β*_{R}, the FWHM_{M}, the mean value of all distances in the rim *δ*_{mean} and an averaged rim width obtained under the hypothesis that both necrosis and the tumour are spherical (*δ*_{s}). The last significant parameter was a measure of the presence of many large distances within the rim, giving some information on its heterogeneity.

The variables identified are novel imaging biomarkers of direct applicability to predicting patient survival in GBM. Moreover, the methodology can be extended to other tumour types displaying necrotic inner parts surrounded by active tumour areas.

### Data accessibility

Data are available at http://matematicas.uclm.es/imaci/molab/home.

### Authors' contributions

Conceptualization, investigation and methodology by J.P.-B., A.M.-G. and V.M.P.-G.; data curation by J.P.-B. and V.M.P.-G.; formal analysis, software and writing original draft by J.P.-B.; funding acquisition, project administration and resources by V.M.P.-G.; supervision, validation and writing review and editing by A.M.-G. and V.M.P.-G.

### Competing interests

We declare we have no competing interests.

### Funding

This research has been supported by the Ministerio de Economía y Competitividad/FEDER, Spain (grant no. MTM2015-71200-R), and James S. McDonnell Foundation 21st Century Science Initiative in Mathematical and Complex Systems Approaches for Brain Cancer (Collaborative award 220020450).

## Acknowledgements

We would like to thank J. M. Borrás and M. Claramonte (Neurosurgery Department), C. López, M. Calvo (Radiology Department) and E. Arregui (Radiation Oncology Department) from Hospital General Universitario de Ciudad Real; J. A. Barcia and L. Iglesias (Neurosurgery Department), J. Avecillas (Radiology Department) from H. Clínico San Carlos, Madrid; C. Velásquez and J. Martino (Neurosurgery Department) from Hospital Universitario Marqués de Valdecilla and Fundación Instituto de Investigación Marqués de Valdecilla, Santander; P. C. Lara and R. Cabrera (Radiation Oncology Department) from Hospital Universitario de Gran Canaria Doctor Negrín; B. Asenjo and I. Herruzo (Radiology Department), M. Benavides (Medical Oncology Department) from Hospital Regional Universitario de Málaga; D. Albillo (Radiology Department), M. Navarro, J. M. Villanueva and J. C. Paniagua (Medical Oncology) and L. Pérez-Romasanta (Radiation Oncology Department) from Hospital Clínico Universitario de Salamanca; B. Meléndez (Pathology Department), Á. Rodríguez (Neurosurgery Department) and R. Moreno (Radiology) from Hospital Virgen de la Salud (Toledo); T. Revert (Radiology Department) from Hospital de Manises; and E. Arana (Radiation Oncology Department) from Instituto Valenciano de Oncología, for providing the data used in this study.

### References

- 1.
Stupp R *et al.*2005 Radiotherapy plus concomitant and adjuvant temozolomide for glioblastoma.**N. Engl. J. Med.**, 987-996. (doi:10.1056/NEJMoa043330) Crossref, PubMed, Web of Science, Google Scholar**352** - 2.
Ellingson BM *et al.*2015 Concensus recommendations for a standardized brain tumor imaging protocol in clinical trials.**Neurooncology**, 1188-1198. (doi:10.1093/neuonc/nov095) PubMed, Web of Science, Google Scholar**17** - 3.
Gatenby RA, Grove O, Gillies RJ . 2013 Quantitative imaging in cancer evolution and ecology.**Radiology**, 8-15. (doi:10.1148/radiol.13122697) Crossref, PubMed, Web of Science, Google Scholar**269** - 4.
Gillies RJ, Kinahan PE, Hricak H . 2016 Radiomics: images are more than pictures, they are data.**Radiology**, 563-577. (doi:10.1148/radiol.2015151169) Crossref, PubMed, Web of Science, Google Scholar**278** - 5.
Zhou M *et al.*2017 Radiomics in brain tumor: image assessment, quantitative feature descriptors, and machine-learning approaches.**ANJR Am. J. Neuro-radiol.**, 208-216. (doi:10.3174/ajnr.a5391) Crossref, PubMed, Web of Science, Google Scholar**39** - 6.
Swanson KR, Rostomily RC, Alvord EC . 2008 Amathematical modelling tool for predicting survival of individual patients following resection of glioblastoma: a proof of principle.**Br. J. Cancer**, 113-119. (doi:10.1038/sj.bjc.6604125) Crossref, PubMed, Web of Science, Google Scholar**98** - 7.
Ramakrishna R *et al.*2010 Imaging features of invasion and preoperative tumor burden in previously untreated glioblastomas: correlation with survival.**Surg. Neurol. Int.**, 40. (doi:10.4103/2152-7806.68337) Crossref, PubMed, Google Scholar**1** - 8.
Zacharaki EI, Morita N, Bhatt P, O’Rourke DM, Melhem ER, Davatzikos C . 2012 Imaging descriptors improve the predictive power of survival models for glioblastoma patients.**AJNR Am. J. Neuroradiol.**, 1065-1071. (doi:10.3174/ajnr.A2939) Crossref, PubMed, Web of Science, Google Scholar**33** - 9.
Mazurowski MA, Zhang J, Peters KB, Hobbs H . 2014 Computer extracted MR imaging features are associated with survival in glioblastoma patients.**J. Neurooncol.**, 483-488. (doi:10.1007/s11060-014-1580-5) Crossref, PubMed, Web of Science, Google Scholar**120** - 10.
Gutman DA *et al.*2013 MR imaging predictors of molecular profile and survival: multi-institutional study of the TCGA glioblastoma data set.**Radiology**, 560-569. (doi:10.1148/radiol.13120118) Crossref, PubMed, Web of Science, Google Scholar**267** - 11.
Wangaryattawanich P *et al.*2015 Multicenter imaging outcomes study of the Cancer Genome Atlas glioblastoma patient cohort: imaging predictors of overal and progresion-free survival.**Neurooncology**, 1525-1537. (doi:10.1093/neuonc/nov117) PubMed, Web of Science, Google Scholar**17** - 12.
Pérez-Beteta J *et al.*2017 Glioblastoma: does the pretreatment geometry matter? A postcontrast T1 MRI-based study.**Eur. Radiol.**, 1096-1104. (doi:10.1007/s00330-016-4453-9) Crossref, PubMed, Web of Science, Google Scholar**27** - 13.
Pérez-Beteta J *et al.*2018 Tumor surface regularity on MR imaging predicts survival and response to surgery in glioblastoma patients.**Radiology**, 171051. (doi:10.1148/radiol.201171051) Crossref, Web of Science, Google Scholar**288** - 14.
Pérez-García VM, Calvo GF, Belmonte-Beitia J, Diego D, Pérez-Romasanta L . 2011 Bright solitary waves in malignant gliomas.**Phys. Rev. E**, 021921. (doi:10.1103/PhysRevE.84.021921) Crossref, Web of Science, Google Scholar**84** - 15.
Cui Y, Tha KK, Terasaka S, Yamaguchi S, Wang J, Kudo K, Xing L, Shirato H, Li R . 2016 Prognostic imaging biomarkers in glioblastoma: development and independent validation on the basis of multiregion and quantitative analysis of MR images.**Radiology**, 546-553. (doi:10.1148/radiol.2015150358) Crossref, PubMed, Web of Science, Google Scholar**278** - 16.
Ingrisch M *et al.*2017 Radiomic analysis reveals prognostic information in T1-weighted baseline magnetic resonance imaging in patients with glioblastoma.**Invest. Radiol.**, 360-366. (doi:10.1097/RLI.0000000000000349) Crossref, PubMed, Web of Science, Google Scholar**52** - 17.
Cui Y, Ren S, Tha KK, Wu J, Shirato H, Li R . 2017 Volume of high-risk intratumoral subregions at multi-parametric MR imaging predicts overall survival and complements molecular analysis of glioblastoma.**Eur. Radiol.**, 3583-3592. (doi:10.1007/s00330-017-4751-x) Crossref, PubMed, Web of Science, Google Scholar**27** - 18.
Lao J, Chen Y, Li ZC, Li Q, Zhang J, Liu J, Zhai G . 2017 A deep learning-based radiomics model for prediction of survival in glioblastoma multiforme.**Sci. Rep.**, 10353. (doi:10.1038/s41598-017-10649-8) Crossref, PubMed, Web of Science, Google Scholar**7** - 19.
Li Q, Bai H, Chen Y, Sun Q, Liu L, Zhou S, Wang G, Liang C, Li Z-C . 2017 A fully-automatic multiparametric radiomics model: towards reproducible and prognostic imaging signature for prediction of overall survival in glioblastoma multiforme.**Sci. Rep.**, 14331. (doi:10.1038/s41598-017-14753-7) Crossref, PubMed, Web of Science, Google Scholar**7** - 20.
Vera L, Pérez-Beteta J, Molina D, Borrás JM, Benavides M, Barcia JA, Velásquez C, Albillo D, Lara P, Pérez-García VM . 2017 Towards individualized survival prediction in glioblastoma patients using machine learning methods.**Neurooncol.**(Suppl. 3), iii84-iii84. (doi:10.1093/neuonc/nox036.317) Google Scholar**19** - 21.
Zhou M, Chaudhury B, Hall LO, Goldgof DB, Gillies RJ, Gatenby RA . 2017 Identifying spatial imaging biomarkers of glioblastoma multiforme for survival group prediction.**J. Magn. Reson. Imaging**, 115-123. (doi:10.1002/jmri.25497) Crossref, PubMed, Web of Science, Google Scholar**46** - 22.
Zhou M, Hall L, Goldgof D, Russo R, Balagurunathan Y, Gillies R, Gatenby R . 2014 Radiologically defined ecological dynamics and clinical outcomes in glioblastoma multiforme: preliminary results.**Transl. Oncol.**, 5-13. (doi:10.1593/tlo.13730) Crossref, PubMed, Web of Science, Google Scholar**7** - 23.
Lee J, Narang S, Martinez J, Rao G, Rao A . 2015 Spatial habitat features derived from multiparametric magnetic resonance imaging data are associated with molecular subtype and 12-month survival status in glioblastoma multiforme.**PLoS ONE**, e0136557. (doi:10.1371/journal.pone.0136557) PubMed, Web of Science, Google Scholar**10** - 24.
Kim HS, Park SB, Lo SS, Monroe JI, Sohn JW . 2012 Bidirectional local distance measure for comparing segmentations.**Med. Phys.**, 6779-6790. (doi:10.1118/1.4754802) Crossref, PubMed, Web of Science, Google Scholar**39** - 25.
Boyer DM, Lipman Y, St. Clair E, Puente J, Patel BA, Funkhouser T, Jernvall J, Daubechies I . 2011 Algorithms to automatically quantify the geometric similarity of anatomical surfaces.**Proc. Natl Acad. Sci. USA**, 18 221-18 226. (doi:10.1073/pnas.1112822108) Crossref, Web of Science, Google Scholar**108** - 26.
Kim H, Monroe JI, Lo S, Yao M, Harari PM, Machtay M, Sohn JW . 2015 Quantitative evaluation of image segmentation incorporating medical consideration functions.**Med. Phys.**, 3013-3023. (doi:10.1118/1.4921067) Crossref, PubMed, Web of Science, Google Scholar**42** - 27.
Boydev C, Taleb-Ahmed A, Derraz F, Peyrodie L, Thiran J-P, Pasquier D . 2015 Development of CBCT-based prostate setup correction strategies and impact of rectal distension.**Radiat. Oncol.**, 83. (doi:10.1186/s13014-015-0386-8) Crossref, PubMed, Web of Science, Google Scholar**10** - 28.
Biasotti S, Cerri A, Bronstein A, Bronstein M . 2014 Quantifying 3D shape similarity using maps: recent trends, applications and perspectives.**Eurographics**135-159. (doi:10.2312/egst.20141039) Google Scholar - 29.
Cecchin D, Poggiali D, Riccardi L, Turco P, Bui F, De Marchi S . 2015 Analytical and experimental FWHM of a gamma camera: theoretical and practical issues.**PeerJ**, e722. (doi:10.7717/peerj.722) Crossref, PubMed, Web of Science, Google Scholar**3** - 30.
Naresk A, Karuna D . 2015 Half-width at half-maximun, full-width at half-maximun analysis for resolution of asymmetrically apodized optical system with slit apertures.**J. Phys.**, 117-126. (doi:10.1007/s12043-014-0828-0) Google Scholar**84** - 31.
Khan JN, Nazir SA, Horsfield MA, Singh A, Kanagala P, Greenwood JP, Gershlick AH, McCann GP . 2015 Comparison of semi-automated methods to quantify infarct size area at risk by cardiovascular magnetic resonance imaging at 1.5T and 3.0T field strengths.**BMC Res. Notes**, 52. (doi:10.1186/s13104-015-1007-1) Crossref, PubMed, Google Scholar**8** - 32.
Louis DN, Ohgaki H, Wiestler OD, Cavenee WK, Burger PC, Jouvet A, Scheithauer BW, Kleihues P . 2007 The 2007 WHO classification of tumours of the central nervous system.**Acta Neuropathol.**, 97-109. (doi:10.1007/s00401-007-0243-4) Crossref, PubMed, Web of Science, Google Scholar**114** - 33.
Bland JM, Altman DG . 1998 Survival probabilities (the Kaplan–Meier method).**BMJ**, 1572-1580. (doi:10.1136/bmj.317.7172.1572) Crossref, PubMed, Google Scholar**317** - 34.
Zacharaki EI, Morita N, Bhatt P, O’Rourke DM, Melhem ER, Davatzikos C . 2012 Survival analysis of patients with high grade gliomas based on data mining of imaging variables.**AJNR Am. J. Neuroradiol.**, 1065-1071. (doi:10.3174/ajnr.A2939) Crossref, PubMed, Web of Science, Google Scholar**33** - 35.
- 36.
Sprent P, Smeeton NC . 2007**Applied nonparametric statistical methods**, 4th edn. London, UK: Chapman & Hall. Google Scholar - 37.
Altman DG . 1991**Practical statistics for medical research**. London, UK: Chapman & Hall. Google Scholar - 38.
Lawrence L . 1998 A concordance correlation coefficient to evaluate reproducibility.**Biometrics**, 225-268. (doi:10.2307/2532051) Google Scholar**45** - 39.
Lin L . 2000 A note on the concordance correlation coefficient.**Biometrics**, 324-325. (doi:10.1111/j.0006-341X.2000.00324.x) Google Scholar**56**