Radiomics-based prediction of local control in patients with brain metastases following postoperative stereotactic radiotherapy

AURORA study

MR imaging and clinical data was collected as part of the “A Multicenter Analysis of Stereotactic Radiotherapy to the Resection Cavity of Brain Metastases” (AURORA) retrospective trial. The trial was supported by the Radiosurgery and Stereotactic Radiotherapy Working Group of the German Society for Radiation Oncology (DEGRO). The inclusion criteria were: known primary tumor with resected BM and SRT with a radiation dose of > 5 Gy per fraction. Exclusion criteria were: interval between surgery and RT > 100 days, premature discontinuation of RT, and any previous cranial radiation therapy (RT). Synchronous non-resected BMs had to be treated simultaneously with SRT. Ethical approval was obtained at each institution (main approval at the Technical University of Munich: 119/19 S-SR).

LF was determined by individual radiologic review or by histologic results after recurrence surgery. FFLF was calculated as the time difference between the end of SRT and LF. If no LF occurred, patients were right-censored after the last available imaging follow-up.

Dataset

In total, we collected data from 481 patients from seven centers. We analyzed four preoperative imaging sequences of each patient: a T1-weighted sequence with and without contrast enhancement (T1-CE and T1), a T2-weighted sequence (T2) as well as a T2 fluid-attenuated inversion recovery sequence (T2-FLAIR). Except for T1-CE, a missing sequence was allowed.

The required data were available for 352 patients. We split the patients into a training cohort with 253 patients from two centers and an external, multicenter, international test cohort with 99 patients from five centers.

Preprocessing

The DICOM (Digital Imaging and Communications in Medicine file format) images were converted to NIfTI (Neuroimaging Informatics Technology Initiative file format) using dcm2niix¹⁹. The MRI sequences were then further preprocessed using BraTS-Toolkit²⁰. First, the sequences were co-registered using niftyreg²¹ and these were then transformed into the T1-CE space. A brain mask was created using HD-BET²² and applied to all sequences to extract only the brain without the surrounding skull. The skull-stripped sequences were transformed into the BraTS space using the SRI-24 atlas²³. Overall, the preprocessing provided co-registered, skull-stripped sequences in a 1 millimeter isotropic resolution in BraTS space.

The missing sequences were then synthesized using a generative adversarial network (GAN). The GAN takes the three available sequences as input and generates the matching missing fourth sequence. We used a GAN which was originally developed for missing sequences in glioma imaging²⁴, but has been proven to work for metastasis imaging^4,5.

Segmentation

All contrast-enhancing metastases and their surrounding edema were individually segmented using the open-source software 3D-Slicer (version 4.13.0, stable release, https://www.slicer.org/)²⁵ by a medical doctoral student (JAB) after undergoing extensive training by a board-certified radiation oncologist (JCP) (7 years of experience). To ensure accuracy, all segmentations for the test cohort were reviewed and manually adjusted by JCP.

To test the feasibility of a fully automated workflow, segmentations generated by a previously trained neural network^4,5 were used as alternative segmentations and compared to the manual segmentations.

As around 25% of patients had multiple BMs, but usually only the largest is resected²⁶, we also determined the largest metastasis with a connected component analysis²⁷ in all patients with multiple BMs and used only that metastasis and its surrounding edema as segmentations for an additional analysis.

Radiomic feature extraction

Radiomic features were extracted with pyradiomics (version 3.0.1, https://github.com/AIM-Harvard/pyradiomics)²⁸ using the Python implementation. The metastasis segmentation was used to extract the T1-CE features, while the edema segmentation was used for the T2-FLAIR features. In total, we extracted 104 original features per segmentation (see Supplemental Table 1 for a list of features and extraction parameters).

Further analysis and modeling were performed in the programming language R 4.2.3²⁹. To adhere to the Image Biomarker Standardisation Initiative (IBSI) standard³⁰, the kurtosis was adjusted by −3. We created nine feature sets in total. Three of these included only radiomic features. The metastasis and edema feature sets were created by extracting the features from the T1-CE sequence and T2-FLAIR sequence, respectively. Both feature sets were merged into a combined feature set. We also created three clinical feature sets with the following clinical features:

• pre-OP feature set: patient age at RT start, Karnofsky performance status (KPS), histology of the primary tumor, location of BM

• post-OP feature set: pre-OP + resection status

• RT feature set: post-OP + concurrent chemotherapy, concurrent immunotherapy and equivalent dose in two Gray fractions (EQD2)

As a seventh feature set, we combined all radiomic features (combined) with the pre-OP feature set to comb+pre-OP.

Multiple publications suggest the predictive value of the brain metastasis volume (BMV) for predicting LF^31–33. Therefore, we created two additional feature sets by adding the cumulative BMV of each patient as an additional feature to the pre-OP set (pre-OP+BMV) and the comb+pre-OP set (comb+pre-OP+BMV).

Intraclass correlation

To identify radiomic features that were susceptible to small changes in segmentation, we generated additional segmentations of all patients in the training cohort using the previously mentioned neural network⁴. Intraclass correlation (ICC (3,1)) was calculated using the R package “irr”³⁴. According to Koo et al., an ICC above 0.75 is considered “good”³⁵. Consequently, this value was employed as a cut-off threshold. Of the 208 features, 173 (83%) had an ICC of > 0.75 and were selected for all further steps. Of the 35 excluded features, the majority (27) were extracted from the edema mask, while only eight excluded features were extracted from the metastasis mask.

All selected radiomic features were normalized by z-score standardization and by applying the Yeo-Johnson transformation³⁶ to transform the distribution of a variable into a Gaussian distribution.

Feature reduction

We applied a minimum redundancy - maximum relevance (MRMR) ensemble feature selection framework implemented in R³⁷ initially proposed by Ding et al.³⁸ as an efficient method for the selection of relevant and non-redundant features.

We created multiple smaller feature sets of the metastasis, edema, and combined feature sets with three, five, seven, nine, eleven, thirteen, and fifteen features each.

We used bootstrapping³⁹ to obtain more reliable results: Feature reduction was repeatedly applied to 1000 bootstrap samples for each set and each number of features. For our final set of features, we ranked the features based on the number of times they were selected.

The best number of features was later determined by nested cross-validation in the training set.

Batch harmonization

To account for differences created by 29 different MRI scanners in our multicenter dataset, we used batch harmonization implemented by neuroCombat⁴⁰. In total, 10 batches were created according to the MRI model names by combining related models. According to Leithner et al.⁴¹, ComBat harmonization without Empirical Bayes estimation provided slightly higher performance in similar machine learning tasks. Therefore, Empirical Bayes was deactivated. Besides the non-harmonized dataset, we created two harmonized datasets: one by only adjusting the means and the other by adjusting means and variances.

Model creation, testing, and patient stratification

For model creation and evaluation, the R package MLR3⁴² was used as a basis. Our prediction target was right-censored time-to-event data, where we used LF as the event and the FFLF or time-to-last imaging follow-up as the time variable for patients with and without event, respectively. We compared three different learners: random forest (RF), extreme gradient boosting (xgboost), and generalized linear models with elastic net regularization learner (ENR).

We implemented nested cross-validation to select the best mode of batch harmonization and the best number of features: For batch harmonization selection, all three datasets were compared while always using the combined feature set with nine features. Five iterations of five-fold nested cross-validation for dataset selection showed no significant difference between the sets with and without batch harmonization (p = 0.3, Kruskal-Wallis rank sum test). Therefore, all further analyses were performed on the base dataset without batch harmonization to avoid unnecessary and potentially distorting preprocessing steps. To select the ideal number of features in each feature set, the nested cross-validation was conducted without batch harmonization. The best average performance was achieved with seven, three, and seven features in the metastasis, edema, and combined sets, respectively. The comb+pre-OP set, which included the seven combined and four clinical features, therefore, had 11 features. The features are listed in Supplementary Table 2.

The parameter tuning was performed using repeated cross-validation. All tuning and selection steps were performed on the training set. To account for the class imbalance (around 1:5 event:no-event), synthetic minority over-sampling was implemented using SMOTE⁴³. We used an implementation in R which is capable of handling numeric and categorical data. The number of samples in the minority class was increased by creating synthetic samples to reach a ratio of 1:2. We only used SMOTE on the training folds in each step of our (nested) cross-validation. This way we ensured that our models were only validated on real patients.

The final models were trained with the best parameters determined by the cross-validation on the whole training set while also using SMOTE to balance the classes. The models were then tested on our multicentric external test cohort.

The 33rd and 66th percentiles of the continuous risk ranks in the training cohort were used as cutoffs for patient stratification. These cutoffs were used to divide the test cohort into three groups according to their predicted continuous risk rank and compare their survival with Kaplan Meier analysis.

Metrics

To account for both timing and outcome, the learners’ performance was quantified using the concordance index (CI)⁴⁴. The 95% confidence intervals are based on 10,000 bootstrap samples. A decision curve analysis was performed to consider clinical consequences with a time endpoint of 24 months⁴⁵. The threshold range was chosen as suggested by Vickers et al.⁴⁶ based on these considerations: Since LF is a severe event and its detection is critical, a lower threshold of 5% seems appropriate. Especially in elderly and multimorbid patients, where additional imaging may be burdensome, an upper threshold of 30% is reasonable.

Baseline clinical models

To create a baseline for comparison with our radiomic models, we first tested the predictive value of two established clinical indices with univariate Cox analysis: the Recursive Partitioning Analysis (RPA)⁴⁷ and the Graded Prognostic Assessment (GPA)⁴⁸ index. They reached a CI of 0.47 and 0.52 in the internal validation, respectively. In external testing, RPA again performed worse with a CI of 0.39 compared to GPA with a CI of 0.44.

Model performance

The performances in the internal validation, as well as in the multicentric external test cohort, are shown in Table 2. To determine the best overall learner, we ranked the performance across all feature sets and found that ENR ranked best, followed by RF and xgboost with mean ranks of 1.4, 1.6, and 2.9, respectively. Therefore, all further experiments were conducted with ENR. For completeness, the results obtained by RF and xgboost are shown in Supplementary Tables 3 and 4. The highest mean CI across all five folds and ten iterations of the cross-validation was achieved with the comb+pre-OP feature set (CI = 0.67).

The comb+pre-OP feature set also led to the highest performance in the external test cohort and achieved a CI of 0.77. While the T1-CE feature set achieved a CI of 0.76, FLAIR was only able to reach 0.50. The three clinical feature sets performed slightly worse than our radiomic feature sets or the combined feature sets: the pre-OP, post-OP, and RT feature sets reached a CI of 0.64, 0.63, and 0.63 in the internal validation, respectively. In external testing, they achieved a CI of 0.70, 0.65, and 0.70, respectively. While adding the BMV to the pre-OP feature set did not change the predictive performance, adding it to comb+pre-OP led to worse results with a CI of 0.72.

For reproducibility, we listed the beta values used by our best model (comb+pre-OP ENR) in Supplementary Table 5. The corresponding calibration curve to this model is shown in Figure 2 (right panel). Furthermore, we calculated the time-dependent area under the receiver operating characteristic curve (AUC) by transforming the crank to an event probability distribution. The proposed model reached a mean of 0.80. Supplementary Figure 1 shows the plotted time-dependent AUC.

• Download figure
• Open in new tab

Figure 2: Kaplan Meier analysis.

We created dichotomous predictions of the comb+pre-OP ENR model by using the 66th percentiles of the continuous risk ranks in the training cohort as cutoffs for patient stratification. We found a significant difference in freedom from local failure (FFLF) between the predicted low- and high-risk groups (p < 0.001) in the multicenter external test cohort. After 24 months, we found a FFLF of 91% and 26% in the groups, respectively.

Patient stratification

Using the cutoffs determined by the training cohort as described above, our comb+pre-OP ENR model was able to significantly stratify the patients into three risk groups with a low, medium, and high risk of local failure (p = 0.0001, Chi-squared Test). A Kaplan-Meier analysis with all three groups is shown in Supplementary Figure 2.

By combining the low- and medium-risk groups into one, we created dichotomous predictions. Kaplan-Meier analysis (Figure 2) illustrates the survival in each risk group. Decision curve analysis using these predictions showed a net benefit of our predictive model compared to treating all patients in the relevant threshold range (Figure 3).

• Download figure
• Open in new tab

Figure 3: Decision curve analysis (left) and calibration curve (right).

Using the same groups as in Figure 2, we found a net benefit of our predictive model compared to treating all patients in the relevant threshold range from five to 30% through decision curve analysis (left). A decision model shows a clinical benefit if the respective curve shows larger net benefit values than reference strategies. The combination of radiomic features derived from the metastasis, edema, and clinical features (comb+pre-OP) resulted in a higher net benefit compared to using only the clinical Pre-OP features and treating all patients or none. The calibration curve on the right was created by transforming the continuous risk rank predicted by the best comb+pre-OP ENR model (in orange) and by the clinical pre-OP ENR model (in blue) to event probabilities at 24 months. Although both models seem to overestimate the actual risk of our patients, the comb+pre-OP model predicts the risk closer to the actual risk.

The relevance of brain metastasis volume

The predictions of our comb+pre-OP ENR model did weakly correlate with the cumulative BMV or BMV of the largest BM (Spearman’s rank correlation: r = 0.246 (p = 0.014) and 0.254 (p = 0.011), respectively).

While cumulative BMV alone was highly predictive in the test cohort, with a CI of 0.76 in a univariate Cox analysis, it only achieved a CI of 0.53 in internal validation. Using the BMV of only the largest BM increased the internal validation and external testing performance to 0.55 and 0.77, respectively. There was no significant difference in the BMV between the training and test cohort (p = 0.64, Wilcoxon rank sum test).

Stratifying our test set into small and large BMs by dividing the set at the median volume resulted in groups with three and 13 events, respectively. Our best model scored a CI of 0.58 and 0.78 in the respective groups. Interestingly, the model significantly risk-stratified the patients in the small BMV group, but not in the high BMV group (corresponding Kaplan-Meier analysis are depicted in Supplementary Figures 3 and 4).

When repeating the training and testing with the radiomic features extracted only from the biggest BM, the ENR learner was able to reach a CI of 0.75 (comb+clinical, Table 3). The results obtained by the RF model even surpassed our previously best model by 0.01 (CI = 0.78, Supplementary Table 4).

End-to-end model using neural network-based automatic segmentations

To test the predictive value of neuronal network-based segmentations and therefore test the feasibility of a fully automated workflow, we conducted an additional parameter tuning and training run with radiomic features extracted from the automatically created segmentations. The results for our ENR learner are shown in Table 3. The best test results with this data were again obtained with the comb+pre-OP feature set (CI = 0.72). Overall, we observed an average decrease in performance by 0.06.