2.1.

Dataset

A retrospective cohort of 101 patients diagnosed with primary (de novo) glioblastoma at the Hospital of the University of Pennsylvania (HUP) between 2006 and 2013 was used in this analysis.²⁹ The molecular characterization of these tumors was not available, since they were diagnosed prior to the most recent World Health Organization (WHO) classification scheme of 2016, hence they can now be categorized as diffuse astrocytoma not-otherwise-specified, according to the WHO 2016 and the C-IMPACT classification.^39,40 These patients were scanned preoperatively using a standardized Adv-mpMRI acquisition protocol of six different modalities, comprising native (T1) and contrast-enhanced (T1-Gd), T1-weighted, T2-weighted (T2), T2-FLAIR, DTI, and DSC MRI volumes. Mean and median patient age was 62.5 and 61.4 years, respectively (range: 22 to 88.6 years). The protocol was approved by the Institutional Review Board at HUP, and informed consent was obtained from all subjects. No randomization method was used for allocating samples to experimental groups.

Patients were treated with standard maximal safe resection followed by adjuvant chemoradiation using temozolomide and six weeks of daily fractionated radiation, five days per week. The details of this cohort have been described previously,²⁹ with a median overall survival (OS) of 13 months, consistent with historical series. The described cohort was divided in patients that survived $<12\text{months}$ (short-survivors, $n=45$) against others, and to patients that survived more than 14 months (long-survivors, $n=45$) against others. We chose these thresholds based on equal quantiles from the median OS ($\sim 13\text{months}$) to avoid potential bias toward one of the survival groups (short-survivor versus long-survivor) and while considering that discrimination of groups should be clinically meaningful. The median OS of the described cohorts is not significantly different from the median survival of glioblastoma patients in several randomized phase III trials, noting that our cohort consists of unselected patients rather than those eligible for such trials.^41,42

2.2.

Statistical Consideration for Survival Analysis

OS was defined as the duration of time between the establishment of diagnosis and the date of death. Mean and median OS was 14.17 and 12.63 months, respectively (range: 0.1 to 68.27, standard deviation: 11.4), and no patient was censored, as all patients have reached the endpoint of death. Kaplan–Meier (KM) curves were generated for the depiction of OS based on the result of each predictive model, as well as for the true classification. The Cox proportional hazards model was used to estimate the hazard ratio of death between groups.⁴³ Based on our sample size, we are powered to detect predictive performance for a classifier of short-survivors (or, equivalently, long-survivors) with area under the curve (AUC) values of 0.65 or larger with 80% power assuming a 5% type I error rate.

2.3.

Equipment and Acquisition Protocol Details

All the MRI volumes used in this study were acquired in the axial plane using a 3-Tesla Siemens Magnetom Trio A Tim clinical MRI system (Erlangen, Germany), according to the standardized advanced acquisition protocol followed at HUP. The acquired T1 volumes had dimensions of $192\times 256\times 192\text{pixels}$, with spatial resolution of $0.976\times 0.976\times 1{\mathrm{mm}}^3$, slice spacing of 1 mm, and inversion time (TI), repetition time (TR), and echo time (TE) equal to 950, 1760 and 3.1 ms, respectively. The dimensions of the axial 3-D T1-Gd volumes were $192\times 256\times 192\text{pixels}$, with spatial resolution of $0.976\times 0.976\times 1{\mathrm{mm}}^3$, slice spacing of 1 mm, and TI, TR, and TE equal to 950, 1760, and 3.1 ms, respectively. The contrast material used in the scans included in this study was either gadodiamide (Omniscan, GE Healthcare, Mickleton, NJ) or gadobenate dimeglumine (MultiHance, Bracco SpA, Milan) and was administered intravenously (IV). The total dose of contrast material was divided into two IV injections to help minimize errors due to potential contrast leakage out of intravascular space. The initial loading dose described the 25% of the total injected contrast material and the second bolus injection the remaining 75%, with a delay of 5 min. The exact contrast material dosage was dependent on patient weight and given on a relative proportion of $0.3\mathrm{mL}/\mathrm{kg}$. The T2 volumes were acquired prior to any contrast administration, and their dimensions were $208\times 256\times 64\text{pixels}$, with spatial resolution of $0.938\times 0.938\times 3{\mathrm{mm}}^3$, slice spacing of 3 mm, and TR and TE equal to 4680 and 85 ms, respectively. The T2-FLAIR volumes were acquired between the initial IV injection and the DSC acquisition. The dimensions of the T2-FLAIR images were $192\times 256\times 60\text{pixels}$, with spatial resolution of $0.938\times 0.938\times 3{\mathrm{mm}}^3$, slice spacing of 3 mm, and TI, TR, and TE equal to 2500, 9420, and 141 ms, respectively. The DSC acquisition was performed during the second IV injection of the contrast material. The dimensions of the DSC-MRI images were $128\times 128\times 20\text{pixels}$, with spatial resolution of $1.72\times 1.72\times 3{\mathrm{mm}}^3$, slice spacing of 3 mm, and TR and TE equal to 2000 and 45 ms, respectively. DTI scans [axial two-dimensional (2-D)] were acquired using a single-shot spin echo planar imaging sequence (variant: segmented $k$-space\spoiled, options: partial Fourier-phase\fat saturation), with 95 phase encoding steps. Following acquisition at $b=0{\mathrm{s}/\mathrm{mm}}^2$ (repeated 3 times), diffusion weighted images were acquired ($b=1000{\mathrm{s}/\mathrm{mm}}^2$) with diffusion gradients applied in 30 directions. The dimensions of the DTI volumes were $128\times 128\times 40\text{pixels}$ ($\text{matrix size}=128\times 128$, $\text{field of view}=220\times 220{\mathrm{mm}}^2$), with spatial resolution of $1.72\times 1.72\times 3{\mathrm{mm}}^3$, slice spacing of 3 mm, flip angle of 90 deg, imaging frequency of 123, and TR and TE equal to 5000 and 86 ms, respectively.

In this study, we refer to basic-MRI acquisition protocol as the one that comprises T1, T1-Gd, T2, and T2-FLAIR, and to advanced protocol as the one that in addition to the basic protocol includes DSC and DTI volumes.

2.4.

Image Preprocessing

All acquired mpMRI volumes were initially converted to NIfTI,⁴⁴ reoriented to the left–posterior–superior coordinate system, and resampled to spatial resolution of $1\times 1\times 1{\mathrm{mm}}^3$ by coregistration to the SRI atlas⁴⁵ using an affine registration, through the Oxford Center for Functional MRI of the Brain (FMRIB) Linear Image Registration Tool (FLIRT)^46,47 of the FMRIB Software Library (FSL).^48–50 All volumes were then skull stripped using an automated method based on a multiatlas registration and label fusion using a template library of 216 MRI scans and their brain masks.⁵¹ This library was used for target specific template selection and subsequent registrations using a framework of multiatlas segmentation utilizing ensembles,⁵² followed by a region-growing step guided by T2, to obtain the final brain mask including the intracranial cerebrospinal fluid (CSF). All mpMRI volumes were then smoothed using a low-level image processing method, namely Smallest Univalue Segment Assimilating Nucleus, in order to reduce high-frequency intensity variations in regions of uniform intensity profile, while preserving the underlying tissue structure.⁵³ The intensity histograms of all modalities of all patients were matched to the corresponding modality of a single reference patient, who was randomly selected and has an average size tumor when compared to the remaining cohort population. This histogram matching approach was repeated in this study for four different randomly selected patients with an average size tumor (compared to the cohort population) and reproducible results were noted, showing that the final results were not dependent on the selection of a specific reference patient. Furthermore, this histogram matching approach has been notably used routinely in other similar studies.^29,54

A set of commonly used DTI measures, namely the tensor’s fractional anisotropy DTI(FA), radial diffusivity DTI(RAD), axial diffusivity DTI(AX), and apparent diffusion coefficient DTI(ADC), were extracted from the DTI volumes.⁵⁵ Furthermore, the DSC-MRI volumes were used to extract maps of the relative cerebral blood volume, peak height, and percentage signal recovery (PSR).^56–58 Hereafter, we treat these DTI measurements and DSC maps as individual imaging volumes.

2.5.

Segmentation of Tumor Subregions

The four modalities of the basic-MRI acquisition protocol (i.e., T1, T1-Gd, T2, and T2-FLAIR) from the preoperative scans of each patient were used to obtain the segmentation labels of the various tumor subregions (Fig. 1). The method used to produce the segmentation labels is called GLISTRboost^59,60 and it defines a hybrid generative-discriminative brain tumor segmentation method. The generative part incorporates a glioma growth model^61–63 and is based on an expectation–maximization framework to segment the brain volumes into tumor (i.e., ET, NET, and ED), as well as healthy tissue labels (i.e., white and gray matter, cerebrospinal fluid, vessels, and cerebellum). During this step, the patient data are coregistered to a standardized atlas coordinate system allowing for estimation of the tumor proportion in various anatomical regions and hence offering location and spatial distribution information. The discriminative part is based on a gradient boosting^64,65 multiclass classification scheme, which was trained on data provided during the brain tumor segmentation (BraTS) challenge^5,66,67 to refine tumor labels based on information from multiple patients. Finally, a Bayesian strategy⁶⁸ is employed to further refine and finalize the tumor segmentation labels based on patient-specific intensity statistics from the four basic modalities available. The derived segmentation labels were considered final after their evaluation and manual revision, when needed, by an expert board-certified neuroradiologist (M. B.).

Fig. 1

Segmentation labels of brain and glioblastoma subregions in multimodal MRI using GLISTRboost.^59,60

2.6.

Radiomic Features

2.7.

Predictive Modeling

The problem of modeling the prediction of survival is approached as a binary classification problem, between long-survivor and short-survivor, using a formulation of two multivariate SVM classifiers (SVC) with a linear kernel, similar to Ref. 29. Moreover, after evaluating the similarity of the results obtained using a linear kernel and a radial basis function kernel, we chose the linear kernel for computational simplicity and for generalization purposes. One SVC model is constructed to distinguish between patients that survived $<12\text{months}$ against others (referred to as “short-SVC”) and another for discriminating patients that survived more than 14 months against others (referred to as “long-SVC”). The decisions of these classifiers are then combined: (a) patients classified by both classifiers as short-survivors are assigned a label of short-survivors, (b) patients classified by both SVC models as long-survivors are assigned a label of long-survivors, (c) patients classified as long-survivors by the short-SVC and as short-survivors by the long-SVC are given the label of the intermediate survivor, and (d) patients for which the two classifiers disagree (i.e., the short-SVC classified them as short-survivors and the long-SVC as long-survivors) are given a label based on which distance from the two SVM hyperplanes is dominant.

The discrimination of subjects into three groups (short-, intermediate-, and long-survivors) is clinically useful for several reasons. Three-class classification schemes are abundant in oncology. The terminology of “low,” “intermediate,” and “high” is commonly used in risk assessments, staging systems, and pathological grading systems, for their ease of translation into clinical decision-making. Recommending therapies for patients requires a clinician’s assessment of patient performance status and life expectancy, which can be notoriously imprecise; a three-class category of “better than expected,” “as expected,” and “worse than expected” is practical for routine clinical use, i.e., when trying to customize therapeutic decisions for the patient sitting in front of the physician.

2.8.

Cross Validation

The generalization performance of our predictive modeling was quantitatively validated using a nested fivefold cross validation scheme to provide unbiased performance estimates. In this fivefold configuration, short- and long-survivors have been proportionally and randomly divided into five nonoverlapping smaller equally sized datasets such that they are statistically representative of the whole, and during each fold, four of these subsets were considered to be the discovery/retrospective cohort and 1 as the replication/prospective cohort, which is unseen for this specific fold. During each fold, the replication/prospective cohort is a different subset. This cross-validation (CV) scheme is similar to analyzing the given data as having independent retrospective and prospective cohorts, but in a more statistically robust manner by randomly permuting across the provided data. Furthermore, nested within each of these folds the feature selection was conducted as previously described, and immediately after the feature selection the SVC parameter (c) was optimized on another fivefold CV scheme.

2.9.

Code Availability

The tool used for coregistration of the brain scans (i.e., FLIRT) is publicly available from the FSL, in: fsl.fmrib.ox.ac.uk. The approach used for skull stripping⁵² and the segmentation method for preoperative mpMRI brain glioma scans, GLISTRboost,^59,60 has been made available for public use through the Online Image Processing Portal (IPP)⁸⁴ of the Center for Biomedical Image Computing and Analytics (CBICA). CBICA’s IPP allows users to perform their data analysis using integrated algorithms, without any software installation, while also freely using CBICA’s high-performance computing resources. It should be noted that we used the publicly available Python package scikit-learn⁸⁵ for the implementation of the gradient boosting algorithm used in GLISTRboost. The Cancer Imaging Phenomics Toolkit (CaPTk),^54,86,87 developed by the CBICA, was used for the (i) seed-point initialization required for GLISTRboost, (ii) denoising approach (SUSAN), as well as (iii) extraction of the radiomic features (i)–(iv) and (vi) described in Sec. 2.6.1. The code source, as well as the executable installer of CaPTk, is available in Ref. 88. The texture features, described in Sec. 2.6.1 as (v), were extracted using the MATLAB^® “radiomics” package,⁷¹ written by one of the two equally contributing first authors of the IBSI. The tumor growth parameters [features (vii)] were extracted through the glioma growth model^61–63 incorporated in the generative part of GLISTRboost. Finally, results were originally obtained by G. S. using in-house code written by G. S., S. B., and G. E., with consultation of H. A. These results were subsequently replicated by S. B. and A. S. at two independent sets of experiments.

3.1.

Accuracy Results

The obtained cross-validated accuracies for classifying glioblastoma patients in short-, intermediate- and long-survivors, when using (a) the Adv-mpMRI and IFP (“IFP advanced” model), (b) the Bas-mpMRI and IFP (“IFP-basic” model), and (c) the Bas-mpMRI and augmented feature panel (AFP) (“AFP-basic” model) were equal to 60,89%, 72.77%, and 74.26%, respectively (Table 1).

Table 1

Accuracy for predicting survival in glioblastoma patients based on various configurations of MRI acquisition protocols and features. IFP and AFP stand for initial and augmented feature panels, respectively.

We calculated individual performance metrics (sensitivity and specificity) for the individual classifiers (short-SVC and long-SVC) within each model, allowing for a receiver operating characteristic analysis, the interpretation of which should be conducted in pairwise manner (short-SVC and long-SVC) for the three-class classification evaluated in this study (Fig. 2). The classification results from these two classifiers were then combined to build a three-class classifier (short-, intermediate-, and long-survivors), as described in Sec. 2.7. We first evaluated the performance of this three-class classification using the IFP advanced model, similar to our previous work.²⁹ This model returned an accuracy of 72.77%, which was improved relative to the performance of the individual short-SVC and long-SVC (Tables 1 and 2). We next evaluated the performance of this three-class classification using the IFP-basic model, which resulted in a notable decrement in performance, returning an accuracy of 60.89%. Finally, we utilized the AFP-basic model, which resulted in a classification accuracy of 74.26%.

Fig. 2

ROC analysis for the short-/long-SVCs of each of the predictive models.

Table 2

Cross-validated performance evaluation metrics for the short-SVC and long-SVC for the three assessed predictive models. IFP and AFP stand for initial and augmented feature panels, respectively.

3.2.

Feature Selection

We analyzed the features selected by the nested fivefold cross validation to determine if the various predictive models actually select features from the varying configurations, as well as to identify which features are selected. Indeed, it is noted that the IFP advanced model selected 69.57% of its features from the basic-MRI protocol and 30.43% from the advanced MRI protocol (Fig. 3). Furthermore, the AFP-basic model revealed that 73% of its selected features were obtained from the AFP and the remaining 27% from the IFP (Fig. 3). The list of selected features as well as their type for the IFP advanced and the AFP-basic model are given in Tables 3 and 4 (Appendix), respectively.

Fig. 3

Proportion of the selected features for each of the predictive models.

Table 3

The exact features selected from the IFP-advanced model.

Table 4

The exact features selected from the AFP-basic model.

3.3.

Kaplan–Meier Analysis

Cox proportional hazard modeling for OS was performed for the classification between short-, intermediate-, and long-survivors based on their real survival classes, as well as on the predictive models. Using the IFP advanced, the $\text{median}\pm \text{standard deviation}$ OS of the short- and long-survivor groups was $9.1\pm 7.7\text{months}$ and $14.9\pm 11.3\text{months}$, respectively, and the hazard ratio for death, for short versus long, was 2.02 (95% CI:1.67 to 2.44). Using IFP-basic, the $\text{median}\pm \text{standard deviation}$ OS of the short- and long-survivor groups was $9.2\pm 8.8\text{months}$ and $14.2\pm 11.8\text{months}$, respectively; hazard ratio for death was 1.74 (95% CI:1.45 to 2.10). Using the AFP-basic, the $\text{median}\pm \text{standard deviation}$ OS of the short- and long-survivor groups was $5.1\pm 3.9\text{months}$ and $16.6\pm 14.8\text{months}$, respectively; hazard ratio for death was 2.84 (95% CI:2.42 to 3.34) (Fig. 4).

Fig. 4

KM curves on the provided patient data for their classification on short- ($<12\text{months}$), intermediate- (between 12 and 14 months), and long-survivors ($>14\text{months}$), based on (a) the real labels, (b) the IFP-advanced model, (c) the IFP-basic model, and (d) the AFP-basic model.