2.1.

Participants

Institutional review board approval was obtained for this study. We first analyzed a discovery cohort of 31 patients with de novo glioblastoma, who were diagnosed at the Hospital of the University of Pennsylvania (UPenn) between 2006 and 2013. Inclusion of subjects was based on the following criteria: (i) age over 18 years, (ii) histopathological tissue diagnosis of glioblastoma (World Health Organization grade IV), (iii) clinical diagnosis of tumor recurrence proven with histopathologic analysis after repeat resection, (iv) no previous history of tumor and resection before the first surgery at University of Pennsylvania, (v) availability of 3 T scanner data, (vi) available preoperative MRI consisting of precontrast axial T1-weighted (T1), T1CE, T2, T2-FLAIR, DTI, and DSC-MRI, and (vii) available postoperative and recurrence time-point MRIs comprising at-least T1, T1CE, T2, and T2-FLAIR. More detail on image acquisition and data demographics can be found in Appendix A. Subjects who had residual tumor after surgical resection were excluded. Residual tumors were defined as any contrast-enhancing areas identified by neuroradiologist (M.B., 16 years of experience) on the immediate postoperative MRI scans captured no later than 48 h after the surgery. The replication cohort comprised 59 patients with de novo glioblastoma who met the aforementioned criteria. The same treatment, i.e., gross total resection of ET followed by temozolomide and radiotherapy was given to all the patients involved in the study.

2.2.

Preprocessing, Segmentation, and Calculation of Perfusion and Diffusion Derivatives

The fractional anisotropy (FA), radial diffusivity (RAD), axial diffusivity (AX), and apparent diffusion coefficient (ADC) were derived from DTI, and relative cerebral blood volume (rCBV) was derived from DSC-MRI. All preoperative MRIs of each patient were coregistered, smoothed, corrected for magnetic field in-homogeneities, and skull stripped.^23–25 The computer-based glioma image segmentation and registration (GLISTR) algorithm²⁶ was used to segment peritumoral ED and two distinct regions of tumor core (TC), i.e., ET and nonenhancing tumor core (NET) (Fig. 1).

Fig. 1

Brain tumor segmentation via GLISTR by using input (a) T1CE, (b) T2-FALIR, (c) T1, and (d) T2 sequences. The segmented pertitumoral edema/tumor infiltration, enhancing tumor, and nonenhancing tumor core are seen overlaid on (e) axial, (f) sagittal, and (g) coronal view of T1CE sequence.

The GLISTR algorithm uses T1, T1CE, T2, and T2-FLAIR sequences. It simultaneously registers a probabilistic atlas of a healthy population to the MRIs with glioblastoma and segments the sequences into ET, NET, ED, and healthy tissue labels. This method uses the expectation maximization (EM) algorithm and incorporates a tumor growth model, a process which modifies the original atlas into a sequence with ET, NET, and ED adapted to best match a given set of patient sequences. The modified atlas is registered into the patient sequences, and then GLISTR utilizes this process for estimating the posterior probabilities of various tissue labels. EM iteratively refines the estimates of the posterior probabilities of tissue labels, the deformation field, and the tumor growth model parameters. The deformation field takes into account the mass effect produced by the tumor. Therefore, GLISTR has a framework for joint segmentation and nonlinear registration that is guided by the atlas. The inverted deformation field warps the patient sequences into the common atlas space. The mass effect gets relaxed in the atlas-normalized sequences; hence, the ventricles get decompressed compared to the original patients’ sequences. This method thus incorporates a tumor growth model to adapt the normal atlas into the anatomy of the patient brain.

2.3.

Radiomic Features

Ideally, one would want to have ground truth of peritumoral tissue that is not infiltrated and tissue that is most-heavily infiltrated, in order to train a classifier. In the absence of such a ground truth, we took an alternative approach, which his discussed next. The heterogeneity of edematous tissue was assessed by defining two regions of interest (ROIs) within the peritumoral edema: near-tumor ROI (N-ROI), adjacent to TC, and far-from-tumor ROI (F-ROI), farthermost from the TC but still within the peritumoral bright-FLAIR ED (Fig. 2). The hypothesis behind using N-ROI and F-ROI was that they are expected to have relatively higher and lower infiltration and, therefore, can serve as references to model relatively more- and relatively less-infiltrated regions, respectively.^7,27 There was not any fixed length of the ROIs; however, F-ROI was generally three times larger in size compared to N-ROI because more space was available for drawing farther from the tumor. The F-ROIs were subsampled at equal intervals to have a balanced training dataset with equal number of voxels of N-ROI and F-ROI.

Fig. 2

Regions of interest: N-ROI, adjacent to TC; and F-ROI at the most distant edge of edema.

A set of radiomic features was extracted at each voxel of the N-ROI and F-ROI by using all the imaging sequences. The features were categorized into five different groups; (i) intensity features: signal intensity of preoperative MRIs; (ii) distance features: distance of a voxel from TC; (iii) statistical features: local mean and median of a voxel in preoperative MRIs; (iv) texture features: local contrast and entropy of a voxel in preoperative MRIs; (v) perfusion temporal dynamics: perfusion time series of each voxel was summarized by five principal components accounting for more than 95% of the signal’s variance. More detail on the features can be found in Appendix B. Prior to feature extraction, all the images were normalized in 0 to 255 range.

2.4.

Discovery and Replication Evaluation

A multidimensional pattern classifier was trained on the above-mentioned features of the voxels of N-ROI and F-ROI using support vector machines (SVM)²⁸ with a Gaussian kernel function. The trained model was then applied on test voxels to calculate the infiltration scores. The infiltration score of each voxel was converted to pseudoprobability of that voxel belonging to N-ROI or F-ROI class by using the Platt’s calibration method.²⁹ In our modeling, this voxelwise map signifies spatial pseudoprobability of infiltration that was named the infiltration index. It may be noted that the infiltration index of a voxel has a value between 0 and 1 representing estimated least-infiltrated and most-infiltrated tissues, respectively.

In the discovery cohort, leave-one-out crossvalidation strategy was adopted wherein an SVM classifier was trained on $n-1$ ($n=31$ for discovery cohort) subjects and was tested on the left out subject. This process was reiterated $n$ times, each time leaving a different subject out. In order to provide realistic estimates of how well the predictive models were likely to generalize to new populations, a model was prepared on all the subjects of the discovery cohort, i.e., $n=31$ and was tested in a totally independent replication cohort of $n=59$ de novo glioblastoma subjects. In particular, the principal component analysis transform and $z$-score parameters (mean, standard deviation) were learned from the training set and then applied to test set (using crossvalidation on discovery cohort and using direct evaluation on replication cohort).

2.5.

Evaluation Method

The evaluation method following the guidelines of Ref. 7 was used here. In brief, evaluation was performed by comparing the estimated infiltration map to two ROIs: recurring ROI (R-ROI) and nonrecurring ROI (NR-ROI). Figure 3 provides a cartoon representation of a tumor with the annotation of recurrence sites. It shows the regions (R-ROI and NR-ROI, both colored) used in the evaluation process.

Fig. 3

Regions used in the evaluation process of the proposed method. R-ROI: region demarcated by the radiologists and additional 3 mm margin around it. White area in the figure shows regions excluded from the evaluation either due to close proximity to TC or to recurrence. NR-ROI: region other than TC, R-ROI, and the regions not included in evaluation due to close proximity (color figure online).

Fig. 4

Top-row shows preoperative MRI sequences and the R-ROI (red) demarcated by the radiologists. The extended blue region shows the 3-mm margin that was also added to R-ROI to account for registration errors. The bottom-row shows recurrence MRI sequences, which were used as a reference to demarcate R-ROI (color figure online).

The following regions were excluded from the evaluation process due to several reasons.

To assess the performance of the proposed infiltration prediction method, the R-ROI and NR-ROI served as positive and negative regions. The binary label maps were obtained by applying threshold of 0.5 on the continuous infiltration probability map (Fig. 5), and true positive (TP), false positive (FP), true negative (TN), and false negative (FN) were calculated by comparing the label map with positive and negative regions. Sensitivity, specificity, accuracy [Eqs. (1)–(3)], and area-under-the-curve (AUC) were calculated from the TP, FP, TN, and FN. In order to account for the different sizes of R-ROIs and NR-ROIs, we calculated balanced accuracy, which is the average value of sensitivity and specificity [Eq. (4)].

Fig. 5

An estimated map for tumor infiltration from preoperative MRIs. Top row (preoperative MRIs): MRIs, comprising T1CE, T2-FLAIR, T1, and T2, on the right, and 3-D view of an estimated infiltration map overlaid on the T1CE image on the left. Red arrow points to the area of high risk. Bottom row (recurrence MRIs): MRIs, comprising T1CE, T2-FLAIR, T1, and T2, on the right, and 3-D view of real recurrence via T1CE image on the left. Middle row: Detailed view of preoperative T1CE MRI with the overlaid predicted infiltration map on the left, and detailed view of recurrence T1CE MRI on the right (color figure online).

In addition, odds ratio³⁰ was calculated to quantify how strongly the estimated preoperative infiltration maps were associated with subsequent recurrence [Eq. (5)]. We considered the infiltration map as a risk factor (exposure) and calculated its association with the outcome, which is a recurrence in our study. The R-ROI and NR-ROI served as outcome of our model, whereas binary label map provided exposed (EX, $>=0.5$) and not exposed (NEX, $<0.5$) regions (Table 1).

Table 1

Outcome and exposure of the proposed infiltration prediction model.

The exposure and outcome were then used to calculate odds ratio.

2.6.

Software Availability

The proposed method is implemented in Insight Toolkit and is integrated in the Cancer Imaging Phenomics Toolkit (CaPTk), which is available in Ref. 31.

3.1.

Infiltration Scores

Table 2 shows the results separately for the subjects of the discovery and replication evaluations. Results for the discovery cohort were obtained using leave-one-out crossvalidation, whereas the results for replication cohort were based on single model trained on the discovery cohort. Classification results include a mean sensitivity of 80.65%, specificity of 87.63%, accuracy of 87.51%, and recurrence odds ratio estimates of 10.22 for tissue predicted to be relatively more infiltrated in the discovery cohort. Similar results were observed for the replication cohort. The proposed method showed a mean sensitivity of 97.06%, specificity of 76.73%, accuracy of 89.54%, and recurrence odds ratio estimates of 13.66 for the replication cohort.

Table 2

Performance of the proposed infiltration prediction model for discovery and replication cohorts.

Figure 5 shows an infiltration map generated for an example subject. The left part of the middle row shows a preoperative T1CE sequence of the subject, and the predicted infiltration map is overlaid on it. The right part shows the T1CE sequence of the same subject with real recurrence. A side-by-side view of both the sequences shows that the predicted recurrence (red region in left) matches well with the real recurrence (bright region in right).

The predictive power of the developed infiltration model was also accessed via receiver operating characteristic (ROC) curve and corresponding AUC. The AUC was 0.83 and 0.91, for discovery and replication cohort, respectively (Fig. 6).

Fig. 6

ROC curves in (a) discovery and (b) replication cohorts. The ROC curves are compared with chance (the black diagonal line).

3.2.

Imaging Features of Regions Involved in Training and Evaluation

Considering the variation in the imaging profile of R-ROI and NR-ROI, we sought to better understand the primary imaging characteristics (features), which distinguish these regions (Fig. 7). The imaging features were significantly different between R-ROI and NR-ROI [Fig. 7(b)], and between N-ROI and F-ROI [Fig. 7(a)]. The imaging features were consistent between N-ROI (hypothesized to be more infiltrated) and R-ROI (recurrent region), and between F-ROI (hypothesized to be less infiltrated) and NR-ROI (nonrecurrent region). The major finding in this direction is that R-ROI, when compared with NR-ROI, revealed higher vascularity and cellularity. In addition, R-ROI has lower signal intensity on T2 and T2-FLAIR sequences and higher signal intensity on T1, thereby suggesting lower water concentration compared to NR-ROI. Each imaging feature contained significant discriminating information as can be seen by the separation between different regions in Figs. 7(a) and 7(b), which when integrated via a machine learning algorithm yielded the reported infiltration predictions. The discriminating power of various individual intensity features calculated by ROC curve analysis was: T1 ($\mathrm{AUC}=0.65$), T1CE ($\mathrm{AUC}=0.63$), T2 ($\mathrm{AUC}=0.70$), T2-FLAIR ($\mathrm{AUC}=0.63$), ADC ($\mathrm{AUC}=0.72$), FA ($\mathrm{AUC}=0.71$), AX ($\mathrm{AUC}=0.69$), and RAD ($\mathrm{AUC}=0.69$).

Fig. 7

Imaging characteristics of the regions involved in training and evaluation of the proposed model. (a) N-ROI and F-ROI, (b) R-ROI, and NR-ROI. The $x$-axis shows the regions and the $y$-axis shows the values of the voxels (scaled in the range of 0 to 1) in each imaging sequence. The $P$ value of ANOVA test showed a significant difference between the two groups in all sequences ($P<0.001$).

The proposed model was also evaluated using conventional imaging sequences such as T1, T1CE, T2-FLAIR, and T2 only, and intensity, texture, statistical, and distance features were used to build the model. The experiment was performed under the same conditions, i.e., leave-one-out crossvalidation on the discovery cohort ($n=31$) and then applying the model built on discovery cohort to an independent replication cohort ($n=59$). The performance of the model for both the cohorts was lower (odds ratio: discovery = 6.59, replication = 8.65) compared to the performance obtained using all the imaging sequences (odds ratio: discovery = 10.22, replication = 13.66). This observation underscores the advantage of using a multiparametric MRI model that integrates synergistic imaging features extracted from various imaging sequences. However, this result also indicates that useful predictions can be obtained via conventional MRI, which is available in all clinics, as long as rich radiomic feature sets are extracted.

4.1.

Implications for Surgical and Radiation Planning

The ability to predict the site of tumor recurrence has numerous potential clinical ramifications. Neurosurgeons have proposed the concept of supratotal resection, in which they administer agents (e.g., Ref. 33) that identify areas of microscopic tumor infiltration at the time of surgical resection of a primary glioma, allowing for increased rates of gross total resection and improved progression-free survival.³⁴ Adding a noninvasive tool to the armament of neurosurgeons may further increase these surrogates for overall survival. The proposed method would enable intensive, yet targeted, surgery and radiotherapy, thereby potentially delaying recurrence and prolonging survival.

Radiation dose escalation trials in the twentieth century uniformly failed to show a survival benefit.^35,36 However, there is renewed interest in these types of trials in the modern era of concurrent chemoradiation with temozolomide. It is theorized that radiotherapy delivery techniques may improve the therapeutic ratio in dose escalation. Identification of the region of interest for radiation dose escalation remains challenging, as a balance is needed between normal tissue toxicity and therapeutic doses. Quantifying a region at highest risk of tumor recurrence provides a potential target for dose escalation based on the very spatial heterogeneity of disease in the peritumoral edema. The high risk areas, found in the infiltration maps, which often appear adjacent to ET or cavity but sometimes are a distance away, may serve as potential regions that could be included in a target volume for radiation dose escalation. We believe that using predictive models to better delineate the region of highest risk can influence the pattern of relapse in patient with glioblastoma, with the potential to improve clinical outcomes such as progression-free survival and overall survival, by directing high-dose radiation primarily toward regions likely to present earlier recurrence while offering relative preservation of lower-risk brain tissue. The proposed method will assist personalization of treatment regimens based on patient-specific features rather than a one-size-fits-all delineation of high-risk regions of interest.

4.2.

Biological Interpretation of the Radiomic Signatures

The main findings of the derived radiomic signature (Fig. 7) indicate that R-ROI (when compared with NR-ROI) has lower signal intensity on native-T2 and T2-FLAIR signals, thereby suggesting lower water concentration. Further, the R-ROI shows relatively higher T1 signal intensity, which would also be consistent with lower water concentration. Finally, the R-ROI shows slightly higher T1CE, which suggests relatively more compromised blood–brain barrier in tissue, also consistent with the characteristics of infiltrating tumor. The diffusion measures provide information that relates to cell density of the peritumoral tissue. Regions with high cellularity tend to have lower ADC³⁷ and higher FA.³⁸ Consistent with the existing literature, the obtained radiomic signature suggests that the R-ROI (when compared with NR-ROI) has lower mean diffusivity and higher FA. The features calculated from DSC-MRI signal, which relate to aspects of tissue vascularization, perfusion, and permeability of blood vessels, were also different among R-ROI and NR-ROI.

Individual assessment of these radiomic signatures displayed differences between R-ROI and NR-ROI, as can be seen via visual inspection in Fig. 7. However, appropriate integration of these radiomic signatures via machine learning methods yielded much better sensitivity and specificity, both in crossvalidation and in independent cohort evaluation, underlining the value of multivariate pattern analysis approaches.

4.3.

Limitations

A limitation of this study is that the data were acquired from a single institution, whereas multicenter data would be beneficial to further and externally validate our infiltration prediction model. However, the use of independent discovery and replication cohorts, along with the use of clinically available imaging sequences, provides confidence that this radiomic signature will generalize well to other institutions and patient populations.