The proposed pipeline (Fig. 1) for GBM overall survival prediction integrates imaging and clinical data through a multi-step process. Multimodal MRI volumes are preprocessed and transferred to a 3D SegResNet for tumor segmentation [24]. High-dimensional features are extracted from the encoder’s latent space, capturing relevant spatial patterns from the MRI data. These features are then reduced via principal component analysis (PCA) to retain the most informative components.

The reduced imaging features, combined with age as clinical variable, are input into three machine learning classification models: (i) random forest, (ii) XGBoost, and (iii) a multi-layer perceptron (MLP). This framework leverages advanced imaging analysis and machine learning to provide precise prognostic predictions for glioblastoma patients.

Fig. 1

Proposed pipeline for glioma feature extraction and survival prediction. Multimodal 3D MRI volumes are pre-processed according to dataset requirements, processed by a 3D SegResNet for feature extraction, and reduced via PCA. The resulting features, combined with clinical data, are used for survival classification

Input MRI 3D volume

Two distinct datasets are used to evaluate the performance of the proposed framework for OS prediction: the BraTS2020 dataset and a proprietary dataset from the Fondazione IRCCS Istituto Neurologico Carlo Besta (Milan, Italy). Both datasets provided imaging and clinical data essential for glioma segmentation and survival analysis. A detailed summary of the datasets will be provided in the subsequent sections. Table 1 outlines the technical specifications of the MRI data, including imaging modalities, resolution, and available segmentation masks. Table 2 provides an overview of the clinical data, summarizing key attributes such as patient demographics and overall survival outcomes, offering a comprehensive basis for the analysis.

Table 1 Technical specifications of MRI data used in this study

Clinical data are presented as median ± standard deviation (Table 2) and their distributions are illustrated in Fig. 2.

Table 2 Technical specifications of clinical data used in this studyFig. 2

Distribution of age at diagnosis and overall survival (OS) for the two cohorts. Upper row: age histograms with median and percentiles for (left) BraTS data and (right) IRCCS Besta data. Lower row: OS (in months) distributions with median and percentiles for (left) BraTS data and (right) IRCCS Besta data

BraTS2020 dataset

The BraTS 2020 dataset is part of the Multimodal Brain Tumor Segmentation Challenge. It comprises pre-operative MRI scans of 369 patients acquired from multiple institutions using a variety of MRI scanners and protocols. Each case includes four MRI sequences—T1, T1Gd, T2, and T2-FLAIR—and tumor annotations identify enhancing tumor (ET), peritumoral edema (ED), necrotic and non-enhancing tumor core (NEC) regions, with expert-validated ground truth segmentations (Fig. 3). All images have been co-registered on the same anatomical template, interpolated to uniform isotropic resolution of 1 mm\(^3\) and dimension of \(240\times 240\times 144\) as reported in Table 1, and finally skull-stripped.

Fig. 3

Example slices from the BraTS 2020 MRI dataset. a T1-weighted post-gadolinium (contrast-enhanced) (T1Gd); b T1-weighted imaging (T1); c T2-weighted imaging (T2); d T2 fluid-attenuated inversion recovery (T2-FLAIR). The segmentation mask is shown on the far right, with the following labels: green (NEC, necrotic tumor core), yellow (ED, peritumoral edema), and red (CE, contrast-enhancing tumor)

In addition to imaging data, the dataset provides clinical metadata, including overall survival (OS) and patient age. For survival prediction, a subset of 236 cases with complete OS labels and corresponding imaging data was selected. This ensured compatibility with the study’s aim of robust training and evaluation of predictive models.

IRCCS Besta dataset

The IRCCS Besta dataset includes MRI scans from 420 glioma patients, collected retrospectively at the IRCCS Besta Neurological Institute of Milan. Imaging data were acquired using either Siemens MAGNETOM Avanto or Philips Ingenia 1.5 T MRI scanners. For each patient, volumetric T1Gd and FLAIR sequences were provided, with a section thickness of 1 mm (Fig. 3). Ground-truth tumor annotations include necrosis (NEC), contrast-enhanced regions (CE), and hyperintense regions (ED), labeled by expert neuroradiologists (Fig. 4; Table 3).

Fig. 4

Example slices from the IRCCS Besta MRI dataset. a T1-weighted post-gadolinium (contrast-enhanced) MRI (T1Gd); b 2 fluid-attenuated inversion recovery (T2-FLAIR). The segmentation mask is shown on the far right, with the following labels: green (NEC, necrosis), yellow (ED, peritumoral edema), and red (CE, contrast-enhancing tumor)

Clinical metadata such as OS and patient demographics are also available, ensuring compatibility for survival prediction tasks (Table 2). To ensure alignment with available OS data and focusing exclusively on high-grade gliomas, a subset of 228 patients was ultimately selected for downstream analyses.

Class distribution for OS prediction

For survival prediction, overall survival was encoded as a binary variable based on a threshold of 400 days (13 months). Cases with survival times \(\le \) 400 days were labeled as “short survival” (Class 0), while those > 400 days were labeled as "long survival" (Class 1). This threshold represents a trade-off between clinical relevance and class balance: it aligns with established survival benchmarks while yielding well-separated groups of comparable size.

To train and evaluate the predictive framework, the dataset was split into training and testing sets using a 70–30 ratio. The training set was used to develop the machine learning models, while the independent test set provided a reliable assessment of performance. Table 3 summarizes the class distribution for each dataset and provides details of the training and testing splits.

Table 3 Class distribution and dataset split for OS prediction among the two datasetsPre-processing and data preparation

Pre-processing was essential to minimize variability introduced by differing imaging protocols, scanner characteristics, and acquisition settings across datasets. Both datasets underwent pre-processing to ensure consistency and compatibility with the 3D SegResNet model, including a resizing of all images to a standardized spatial dimension of \(224\times 224\times 144\).

The BraTS2020 dataset provided pre-processed images, except for intensity normalization, ensuring standardized inputs across scans. In contrast, the IRCCS Besta dataset required additional pre-processing to align with the BraTS2020 specifications. Steps included multi-modal imaging co-registration, skull-stripping, resampling and intensity normalization. All these steps were performed via known softwares for brain MRI manipulation; in particular, multi-modal image registration was performed using the Elastix [25] package, skull-stripping was carried out with the HD-BET [26] tool, and image resampling was executed via the SimpleITK [27] library’s ResampleImageFilter.

A flowchart of the preprocessing workflow applied to the IRCCS Besta dataset is shown in Fig. 5.

Fig. 5

Preprocessing workflow applied on 3D MRI volumes of IRCCS Besta dataset

The 3D SegResNet model is designed to process input volumes of fixed dimensions (224 \(\times \) 224 \(\times \) 144). To ensure compatibility, all MRI volume from both datasets were uniformly cropped to these dimensions during pre-processing, aligning the spatial resolution of the data with the network’s architectural requirements while maintaining critical anatomical information.

As a data augmentation technique for both datasets, we applied random flipping along the axial, coronal, and sagittal axes with a 0.5 probability, and by performing these flips at runtime during each training epoch—presenting each sample in multiple augmented forms—we enriched dataset variability, simulated diverse imaging conditions, and ultimately improved the model’s ability to generalize to unseen data.

Segmentation and feature extraction: 3D SegResNet with latent space analysis

The segmentation task and feature extraction are performed using the 3D SegResNet architecture, available within the MONAI framework [28]. Pretrained on the BraTS2018 dataset, this network is specifically designed for brain tumor segmentation from MRI, producing voxel-wise segmentation maps of glioma subregions, such as necrosis, peritumoral edema, and contrast-enhancing tumor regions.

The network processes multimodal MRI inputs (T1, T1Gd, T2, T2-FLAIR) as separate channels, using complementary information to accurately segment the glioma subregion. The 3D SegResNet adopts a modified 3D U-Net architecture with residual connections, enhancing feature propagation and training stability. Its encoder-decoder structure captures hierarchical spatial features through downsampling and reconstructs segmentation masks via upsampling with skip connections. Residual blocks mitigate vanishing gradients, ensuring effective learning in deep networks.

3D SegResNet fine-tuning

To address the domain change between the two public BraTS releases and assess performance in our private data set, we initialized the 3D SegResNet with weights pre-trained in BraTS2018 and subsequently fine-tuned the entire network on the BraTS2020 training set for 50 epochs (learning rate = 1 \(\times \) e\(^\), weight decay = 1 \(\times \) 10\(^\)) using the same augmentation pipeline described above. Optimization was performed with Adam and a CosineAnnealingLR scheduler.

Inference was then carried out using these fine-tuned weights under three input configurations, in order to perform a fair comparison of our datasets:

BraTS4CH: BraTS2020 data, 4-channel input (T1, T1Gd, T2, T2-FLAIR).

BraTS2CH: BraTS2020 data, 2-channel input (T1 = 0, T1Gd, T2 = 0, T2-FLAIR).

Besta2CH: IRCCS Besta data, 2-channel input (T1 = 0, T1Gd, T2 = 0, T2-FLAIR).

The dice similarity coefficient (DSC) was calculated for the three output channels: tumor core (TC), whole tumor (WT) and enhancement tumor (ET) across the validation splits. Table 4 reports the mean DSC values for each channel.

Table 4 Mean dice similarity coefficient for the output channels under the three input configurationsPrincipal component analysis

Beyond segmentation, this study uses the encoder’s latent space to extract high-dimensional features from the multimodal MRI inputs. These features provide a compact, abstract representation of the imaging data, capturing complex spatial and intensity patterns that are challenging to derive manually. After extracting the encoder latent features for each cohort, the principal component analysis (PCA) was applied independently to each dataset to retain the variance patterns specific to the cohort while reducing dimensionality. This step ensures efficient integration of encoder-derived features into the survival prediction pipeline, improving computational efficiency while preserving critical prognostic information.

By combining segmentation, feature extraction from the encoder’s latent space, and PCA-based dimensionality reduction, this approach provides a robust framework for both tumor characterization and survival prediction. The encoder of the 3D SegResNet generates 256 high-dimensional feature maps, each with spatial dimensions of 14 \(\times \) 14 \(\times \) 10. Each feature map represents a processed and abstract representation of the input data, emphasizing relevant patterns for the network’s task.

To reduce the dimensionality of these high-dimensional feature maps while retaining essential information, a summarization strategy was applied. Specifically, each volumetric feature map was condensed into a single scalar value by computing its median intensity along the spatial dimensions.

This approach preserves key information while minimizing redundancy and ensuring compact feature representation. The median was chosen because of its well-documented robustness to outliers, particularly in medical imaging data, where variability and noise are common challenges. It has been shown to effectively preserve diagnostic features, such as edges, while reducing noise and dimensionality, ensuring that essential information is retained for downstream analysis [29].

OS classification

The dataset was randomly divided into a 70% training set—with an internal cross-validation loop for hyperparameter selection—and a 30% hold-out test set reserved for final evaluation of performance. Hyperparameter optimization was carried out using scikit-learn’s GridSearchCV [30], restricting all searches to the training folds in order to minimize overfitting, by applying a fivefold cross- validation. F1-score was adopted as scoring metric for hyperparameters selection.

Three classification algorithms were implemented for OS prediction:

1.

Random forest (RF): RF works by constructing multiple decision trees during training and combining their outputs to enhance predictive accuracy and reduce overfitting. GridSearchCV was employed to explore variations in the number of trees, maximum depth, splitting criterion, minimum samples per split and leaf, features considered at each split, and class-weight balancing.

2.

XGBoost: it is a gradient boosting algorithm known for its efficiency and accuracy in structured data classification. For XGBoost, a similar set of parameters was employed with the extension of the search to include learning rate and subsample fraction; during fitting, both training and validation losses were tracked and early stopping halted training after ten rounds without improvement.

3.

Multi-layer perceptron (MLP): a deep learning model designed to capture complex interactions between input features. Hyperparameters were tuned over hidden-layer architectures, L2 regularization strengths, initial learning rates, and activation functions (relu, tanh, logistic). Early stopping—with 20% of the training data held out internally and a patience of ten epochs—prevented overtraining over a maximum of 1000 iterations.

Once optimal hyperparameters were identified for each algorithm, predictive performance on the independent test set was assessed by comparing paired F1 scores via the Wilcoxon signed-rank test, thereby providing a statistically robust evaluation of relative model efficacy.