Journal of The Royal Society Interface
Published:https://doi.org/10.1098/rsif.2018.0503

    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.

    Figure 1.

    Figure 1. Examples of GBMs as seen in post-contrast T1-weighted MRIs. (a,c) Shown in the axial plane. (b,d) Shown in the coronal plane.

    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 [1520]. Finally, other authors have explored different types of spatial imaging biomarkers based on the quantification of tumour subregions [2123].

    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 PI,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 PI,E and allows for the reconstruction of the surfaces.

    Figure 2.

    Figure 2. Volumetric reconstruction of tumours from MRIs. (a) Axial, sagittal and coronal MRIs of the patient’s head showing the GBM location. (b) Details of the tumour along the three orthogonal planes. (c) Segmentations (green) of the contrast-enhancing regions for all the axial slices containing tumour voxels. (d) Three-dimensional reconstruction of the tumour from the individual slices.

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

    Quantifying the space contained between surfaces is a topic of interest in many scientific fields [2428]. 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 RR3 the space contained between the surfaces ΣE and ΣI.

    To describe R, we first constructed the set of distances between the NE points of PE and those NI points in PI. To do so, we first computed the distances from every point in PE to those in PI. For every point in PEkPE, we looked for the point PIq providing the minimal distance of all points in PI and defined δEk=PEkPIq. Next, we constructed the set

    ΔE={δEk,k=1,,NE}.2.1

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

    Figure 3.

    Figure 3. Three-dimensional rim width colour map of the characterization of R as a set of distances plotted on ΣE.

    Next, the process was repeated for all the points in PI. For each of the NI points in PI={PIq}q=1NI, we obtained the minimal distance vector to PE (δIq). 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 δIq be sufficiently far from all points in PI. To make this more precise, the angles θq,j spanned by the vector δIq and the one generated by the point PIq and the other points PIj were constrained to be larger than a small number ε (taken here to be 0.035). In other words, we define

    ΔI={δIq:θq,j=arccos((δIq,(PIqPIj))δIq(PIqPIj))>ϵ,j=1,,NI,jq},2.2
    where (,) stands for the scalar product of two vectors in R3.

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

    Figure 4.

    Figure 4. Process of exclusion of vectors δIq crossing ΣI. (a) Two-dimensional example with a region of interest where crossings occur. (b) Detail showing an example of vectors that are not included in the set ΔI.

    Finally, we define the final set of distances Δ characterizing R:

    Δ=(ΔEΔI)={δn}n=1N.2.3
    Thus, the final outcome of the procedure is a set of distances characterizing the ‘rim width’. Figure 5 shows an example of surfaces ΣE,I and one of the minimal distances between them included in the set Δ.
    Figure 5.

    Figure 5. An example of surfaces ΣE and ΣI together with one of the minimal distance vectors δi, joining them.

    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 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. δmean=n=1Nδn/N.

    • Maximum distance (δmax). The largest measurement in the set δmax = maxj=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 (VT), and the necrotic (or inner) volume (VI) as the volumes enclosed by the surfaces ΣE and ΣI, respectively. The CE volume is then defined as VCE = VTVI. 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

      δs=(34π)1/3[(VCE+VI)1/3VI1/3].2.4
      This 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 [2931]. 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(δl)=f(δr)=12f(δmode), then

      FWHM=δrδl.2.5a
      FWHMN and FWHMM 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.
      FWHMN=FWHMδmax2.5b
      and
      FWHMM=FWHMδmaxδmin.2.5c

    • Standard modal deviation (Dmode). This measurement quantifies the deviation of all the distances with respect to δmode. It is given by the equation

      Dmode=1niN(δiδmode)2.2.6

    • Mode variation coefficient (CVmode). This relates the value of δmode and its dispersion Dmode. 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

      CVmode(%)=Dmodeδmode×100.2.7

    • Standard mean deviation (Dmean). This quantifies the amount of variation of all distances with respect to δmean:

      Dmean=1niN(δiδmean)2.2.8
      This quantity is computed directly on Δ and thus does not depend on the class size t.

    • Mean variation coefficient (CVmean). This quantity relates the value of δmean and its dispersion Dmean as given by the equation

      CVmean(%)=Dmeanδmean×100.2.9

    • Geometrical heterogeneity (Gh y GH). These two measurements quantify the ratio of the dispersion of the largest measurements over the whole set of those in Δ. Their mathematical expressions are

      GH=δmaxδQ3δmax2.10a
      and
      Gh=δmaxδQ3δmaxδmin.2.10b

    • Width between mode and mean (WMM). Is the distance between δmode and δmean,

      WMM=δmeanδmode.2.11a

    • From this measurement, related normalized quantities can be defined as

      WMMN=WMMδmax2.11b
      and
      WMMM=WMMδmaxδmin,2.11c

    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

      βE,t:Δt(δ)KeβEδδΔt,δδmode,2.12
      using 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

      βR,t:Δt(δ)KδβRδΔt,δδmode,2.13
      where 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

      HWHM=δrδmode.2.14a
      Normalized versions of this quantity by δmax and δmin are given by the following equations
      HWHMN=HWHMδmax,2.14b
      and
      HWHMM=HWHM(δmaxδmin).2.14c

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

    Figure 6.

    Figure 6. (a) Quantities obtained from the distribution Δ0.5. (b) Analysis of the distribution decay using both exponential (equation (2.12)), and rational (equation (2.13)) fits.

    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.

    Table 1. Summary of patient data.

    patients
    number of patients (censored) 311 (27)
    age (years) median [range] 63 [19 − 86]
    sex (M, F) (44%,56%)
    survival (months), median [range] 12.76, [0.13, 82.97]
    type of resection
    total 149(47.91%)
    subtotal 113(36.33%)
    biopsies 49(15.76%)
    type of treatment
    chemotherapy (CT) 5(1.61%)
    radiotherapy (RT) 27(8.68%)
    CT and RT 241(77.49%)
    no treatment 38(12.22%)
    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.

    Figure 7.

    Figure 7. Kaplan–Meier plots for the variables FWHMM and βR. Patient groups are split by thresholds: (a) FWHMM = 0.34129 for t = 0.75 mm. (b) FWHMM = 0.24696 for t = 0.5 mm. (c) βR = 0.02764 for class width t = 0.75 mm. (d) βR = 0.01450 for class width t = 0.5 mm.

    Figure 8.

    Figure 8. Kaplan–Meier curves for patients according to their δmax, δmean and δs values. (a) Patients with δmax > 1.02280 cm (blue), patients with δmax < 1.02280 cm (red). (b) Patients with δmean > 0.30063 cm. (blue), patients with δmean < 0.30063 cm (red). (c) Patients with δs > 0.3781 cm (blue), patients with δs < 0.3781 cm (red).

    Table 2. Summary of the results obtained in the survival analysis. n.s., not statistically significant; n.a., not applicable (variables not dependent on the class width t).

    variable units t threshold p HR (CI95%) dif. surv.
    tumour size measurements
    VT cm3 n.a. 28.55 0.107 n.s. n.s.
    VCE cm3 n.a. 17.12 0.059 n.s. n.s.
    VI cm3 n.a. 5.71 0.223 n.s. n.s.
    dmax3D 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
    GH n.a. 0.61839 0.099 n.s. n.s.
    Gh n.a. 0.61667 0.116 n.s. n.s.
    Dmean cm n.a. 0.12502 0.154 n.s. n.s.
    CVmean % n.a. 41.5340 0.099 n.s. n.s.
    Dmode cm 0.75 0.21986 0.180 n.s. n.s.
    cm 0.50 0.12930 0.137 n.s. n.s.
    CVmode % 0.75 45.1020 0.251 n.s. n.s.
    % 0.50 57.3280 0.270 n.s. n.s.
    FWHMN 0.75 0.14422 0.072 n.s. n.s.
    0.50 0.24548 0.068 n.s. n.s.
    FWHMM 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
    WMMN 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
    WMMM 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
    HWHMN 0.75 0.21315 0.053 n.s. n.s.
    0.50 0.17761 0.028 1.444 (1.039–2.007) 3.50
    HWHMM 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 FWHMM. The threshold found FWHMM,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 (FWHMM,th = 0.24696, p = 0.032, HR = 1.305 and same survival difference). The corresponding Kaplan–Meier plots are shown in figure 7a,b. The variable WMMN 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, WMMM 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 7c,d.

    Both HWHMN and HWHMM 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 8a). δ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 8b. 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 8c.

    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 Dmode 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.

    Figure 9.

    Figure 9. Colour map showing the Spearman correlations between the study variables. Significant (p < 0.001) correlations with |ρ| ≥ 0.75 were considered to be strong.

    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 [2428]. 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 FWHMM 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 FWHMM, 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.

    Footnotes

    Published by the Royal Society. All rights reserved.

    References