Association of graph-based spatial features with overall survival status of glioblastoma patients

Glioblastoma is the most common malignant brain tumor with less than 15 months median survival. To aid prognosis, there is a need for decision tools that leverage diagnostic modalities such as MRI to inform survival. In this study, we examine higher-order spatial proximity characteristics from habitats and propose two graph-based methods (minimum spanning tree and graph run-length matrix) to characterize spatial heterogeneity over tumor MRI-derived intensity habitats and assess their relationships with overall survival as well as the immune signature status of patients with glioblastoma. A data set of 74 patients was studied based on the availability of post-contrast T1-weighted and T2-weighted fluid attenuated inversion recovery (FLAIR) image data in The Cancer Image Archive (TCIA). We assessed the predictive value of MST- and GRLM-derived features from 2D images for prediction of 12-month survival status and immune signature status of patients with glioblastoma via a receiver operating characteristic curve analysis. For 12-month survival prediction using MST-based method, sensitivity and specificity were 0.82 and 0.79 respectively. For GRLM-based method, sensitivity and specificity were 0.73 and 0.77 respectively. For immune status, sensitivity and specificity were 0.91 and 0.69, respectively, for the GRLM-based method with an immune effector. Our results show that the proposed MST- and GRLM-derived features are predictive of 12-month survival status as well as the immune signature status of patients with glioblastoma. To our knowledge, this is the first application of MST- and GRLM-based proximity analyses for the study of radiologically-defined tumor habitats in glioblastoma.

methods such as a minimum spanning tree (MST) construction and graph run-length matrices (GRLM).In previous breast cancer studies, graph (MST-derived) features have been used to understand the proximity relationships between distinct immunohistochemical entities (cell types) within hematoxylin and eosin (H&E) pathology slides 18 .These MST-based features successfully distinguished samples of high and low lymphocytic infiltration extent with a classification accuracy greater than 90% in breast cancer 18 .Based on the success of this approach, we hypothesized that characterizing the spatial relationship between radiologically-distinct habitats of a tumor using MST or GRLM approaches might have predictive value for underlying clinical outcome in glioblastoma.
Another graph-based characterization called GRLM 19 was also used to characterize the spatial heterogeneity of a tumor to predict clinical outcome in glioblastoma, based on the idea of gray-level run-length matrices.The gray-level run-length method is one of popular methods extracting high order statistical features in texture analysis 20 .In this study, we used GRLM to compute the runs of radiologically defined tumor habitats on a region of interest (ROI) image to obtain a run-length matrix instead of counting runs of voxel intensities.In addition, Researchers have been increasingly studying the role of the immune system in the development, progression, and potential treatment of glioblastoma.The immune microenvironment within the tumor, which includes various immune cells and signaling molecules, plays a critical role in shaping the tumor's behavior and response to treatment.In this study, the immune gene signature status was determined for each patient using single-sample gene set enrichment analysis (ssGSEA), based on the established immune effector (IE) and immune suppressor (IS) responses [21][22][23][24][25] .The ssGSEA process assigned an enrichment score to measure how genes within a specific gene set are collectively upregulated or downregulated within a sample.
The purpose of this study is to evaluate the prognostic significance of higher-order spatial proximity characteristics from habitats using quantitative metrics derived from graph-based methods.We aim to investigate the association of MST and GRLM-derived spatial proximity features of these tumor habitats with overall survival of glioblastoma patients as well as predicting immune signature status.In the broader context, this study aims to investigate relationships between imaging and phenotypic characterizations of the tumor, thereby augmenting the foundation for population-based correlation studies in glioblastoma.

Data
A data set of 74 patients was studied based on the availability of post-contrast T1-weighted and T2-weighted fluid attenuated inversion recovery (FLAIR) image data in The Cancer Genome Atlas (TCGA), The Cancer Imaging Archive (TCIA-http:// www.cance rimag ingar chive.net/).The dataset consisted of 25 female and 49 male patients with de-novo (primary) glioblastoma, all of whom were classified as IDH-wild type astrocytoma.The patient demographics are summarized in Table 1.The mRNA expression data and clinical data such as survival information for these cases 26 were obtained from the cBioPortal for Cancer Genomics (http:// www.cbiop ortal.org).In this study, the previously defined immune effector and immune suppressor response 21 were used to derive the immune gene signature status using single-sample gene set enrichment analysis: ssGSEA for each patient [22][23][24][25] .Basically, the IE response activates immune cells to eliminate threats like pathogens and abnormal cells.Immune cells like T cells, B cells, and natural killer cells are recruited to target and remove harmful entities.Conversely, the IS response moderates immune reactions to prevent excessive active activity, avoiding tissue damage and chronic inflammation.Regulatory T cells and immune checkpoint molecules contribute to maintaining immune balance.The study employs ingle-sample gene set enrichment analysis (ssGSEA) to evaluate these responses in patients.This technique assesses how genes linked to these responses are regulated within patient samples, providing insights into their immune gene signature status.
For predicting immune signature status, 34 patients were used based on availability among 74 patients.The MR images were preprocessed with registration, non-uniformity correction using N3 27 , voxel isotropic resampling, and intensity normalization 28 before subsequent analysis.Registration of the T1 post-contrast image and T2 FLAIR image along with non-uniformity correction for MRI-artifacts were performed using the Medical Image Processing, Analysis, and Visualization (MIPAV) software 29 .The FLAIR MR image is registered to the T1 post-contrast image using affine transformation with 12 degrees of freedom along with trilinear interpolation.Segmentation was a semi-automated process for which MITK3M3 Image analysis toolkit was used.The clinicians used this tool to contour the tumor region on multiple slices with interpolation performed to obtain a 3D volumetric tumor mask.This step was performed independently on T1-post contrast as well as T2-FLAIR.In all the processes, readers were blinded to clinical/molecular characteristics.Voxel isotropic resampling was performed using the NIFTI toolbox in MATLAB to make voxel sizes isotropic (1 mm).The resulting T1 postcontrast image, T2 FLAIR image after preprocessing, and the segmentation results on T2 FLAIR image using MITK is shown in Fig. 1.Overall survival (months) 5.9 ± 3.0 24.8 ± 12.7 16.9 ± 13.6

Radiologically defined tumor habitats
In this study, we used two different MR sequence images (T1 post-contrast image and T2 FLAIR image).For each MRI image-sequence, the voxels within the tumor ROIs were separated into low and high intensity group (habitat) using the two-component GMM.Four binary masks were prepared from these groups.A two-dimensional grid line was overlaid on each binary mask.The grid lines were equally spaced with a distance of 8 voxels (8 mm × 8 mm), chosen empirically.The voxel intensity values within the ROI for each patient were scaled to lie between zero and one by linear transformation, and then fitted using a two-component Gaussian mixture model (GMM) 30,31 .The threshold between two Gaussian groups was determined by calculating the average of the means of the two Gaussian populations underlying the low intensity voxel group and high intensity voxel group within the tumor, following a process similar to prior work 31 .

Minimum spanning tree (MST) and graph run-length matrices (GRLM)
A spanning tree of a graph represents a tree that connects all the vertices.Each edge has an associated weight (or length), the weight of a tree represents the sum of all weights of its edges.Then, the MST can be defined as a spanning tree with the minimum weight of a tree among any other spanning trees.
The gray-level run-length matrix M(i, j) is a popular method in texture analysis that measures the variation of the voxel intensities to quantify intuitive qualities such as smoothness, coarseness, and roughness.Gray-level run-length can be defined as the number of runs with voxels of gray level i and run length j in a given direction.Basically, run-length matrix provides the coarseness of a texture in a specified direction.Runs of data represent sequences in which the same data value, gray level intensity, occurs in many consecutive voxel elements.In general, fine texture or high frequency tends to have more short runs with similar gray level intensities and coarse texture or low frequency tends to have longer runs of similar gray level intensities.

Feature extraction with minimum spanning tree (MST)
We computed the coordinate of the centroid that specifies the center of mass of the region from each habitat inside of the small bounding grid box 32 .Coordinates from all grid boxes in each habitat were combined into one map.An MST was constructed across all these co-ordinates as vertices.Thus, we have four MSTs (one for each habitat: T1 high intensity group, T1 low, T2 high, or T2 low intensity group, respectively).Figure 2 shows MST of each of the four groups on an overlapped T1 and T2 ROI map.The mean, median, standard deviation, skewness, kurtosis, min/max ratio, and disorder of the branch lengths in MST were computed to obtain a set of seven features 18 (for each habitat.Expressions for the features are listed below: Mean edge weight, f µ : where w i is an individual weight or branch length in MST.Median edge weight, f median , is the value separating the higher half of branch lengths from the lower half.
Standard deviation of the distribution of edge weights, f σ : (1) Skewness of distribution of edge weights, f skewness : Kurtosis of distribution of edge weights, f kurtosis ,: (2) (3) www.nature.com/scientificreports/Finally, the min/max ratio, f r , is the ratio between maximum of W divided by the minimum value of W where W is a set of weights (branch lengths), W = {w 1 , w 2 , . . ., w n }.'Disorder' is the standard deviation f σ divided by the mean value of W, f µ .

Feature extraction with graph run-length matrices (GRLM)
Similar to the construction of the MST, we computed the coordinate of the centroid from each habitat inside of the rectangular grid box.For the extraction of GRLM features, we combined all coordinates from all four habitats into one set (map).Again, these four habitat groups represent T1 high-and low-intensity voxel groups and T2 high-and low-intensity voxel groups.Figure 3 shows an ROI map with all coordinates (across all four habitats) as vertices.Then, we constructed a Delaunay triangulation 33 to connect these vertices.There are several ways to triangulate any given set of points and a Delaunay triangulation is one of the most widely used in scientific computing of various applications.In the Delaunay triangulation, all triangles for a set of points will have empty circumscribed circles 33 .
After constructing the ROI map with the Delaunay triangulation, a run-length matrix was computed.In this study, we used the GRLM method in the manner proposed by Tosun et al. for histopathological image segmentation 19 because this aims to represent the spatial separation between point set entities.GRLM G(t, l) can be defined as the number of graph-edge runs with an edge type t and a path length l for a single node.The algorithm starts from the initial node to the furthermost node in the path within a circular window.Figure 4 shows the calculation of a graph run-length matrix for a single node.For the entire region, the algorithm accumulates the run-length matrices of the nodes in an ROI.
For extracting GRLM-based features, we used six measurements: short path emphasis (SPE), long path emphasis (LPE), edge type nonuniformity (ETN), and path length nonuniformity (PLN).The expressions for these features are listed below: Short path emphasis (SPE) and SPE(t) for each edge type t: where n r is the total number of runs in the GRLM and n r (t) is the total number of runs corresponding to edge type t.

Classification of immune signature status
The immune signature scores from the TCGA cohort were dichotomized at the median value; either an up regulated as signature score > median value or down regulated as signature score ≤ median value.This binary designation is used as the class label in the classification task.These gene signatures associated with immune status in glioblastoma, immune effector response and immune suppression response, have been previously validated within glioblastoma and were evaluated in every patient in our dataset 21 .

Statistical analysis
A total of 28 MST based features (seven MST-features from each of the four habitats) and 24 GLCM based (SPE, LPE, ETN, PTN, and 10 SPE(t) and 10 LPE(t)) features for 10 edge types were extracted and analyzed in this study.Survival was dichotomized at the 12 month time point (i.e.≤ or > 12 months) based on the median overall survival duration, which is ranging from 9 to 15 months 1-3 and imbalance in sample size of the two survival groups was controlled using class-proportional sampling 34 in Waikato Environment for Knowledge Analysis (WEKA v3.7.12) 35 .Classification of the survival labels using all of the MST features was performed using random forest (RF) classification 36 (10,000 trees using random forest classifier within WEKA).The important features were identified by using the Gini index, which calculates the amount of probability of a specific feature that is classified incorrectly when selected randomly.
In this study, we split the data into five subsets (5-folds) and rotating the training and validation among them using the fivefold cross-validation.Basically, cross-validation is one of the most powerful statistical techniques used to assess the performance of a predictive model by dividing the available dataset into subsets and using these subsets for both training and testing.The primary goal of cross-validation is to estimate how well a model will generalize to new, unseen data.The model performance was evaluated using the receiver operating characteristic (ROC) curve and the area under ROC curve (AUC).The ranking of the features was also obtained using a classifier-based attribute evaluator (where RF was set as the classifier) within WEKA.Sample size imbalance between classes is handled using class-proportional sampling 34 .The GRLM and MST features were computed in MATLAB (R2020a, The MathWorks, Inc.) using the formulas listed in this paper.The ROC, AUC, and classification using the RF classifier were computed using R (R Foundation for Statistical Computing, Vienna, Austria).

Results
Classifier models were obtained using the random forest classification approach 36 with fivefold cross-validation for the prediction of 12-month overall survival status based on all of the features derived from MST and GRLMbased methods, respectively.Figure 5 shows the ROC curves for the classifiers for MST and GRLM.The optimal cutoff point was determined by maximizing the sum of sensitivity and specificity.The results for the area under the ROC curves, the true positive rate, and true negative rate were summarized in Table 2.The area under the ROC curves was 0.832 for MST and 0.773 for GRLM.The accuracy was 80.7% for MST and 74.7% for GRLM computed using Eq. ( 11) where TP, FP, TN, and FN represent true positive, false positive, true negative, and false negative, respectively.
The most important features based on rank (top five) are ratio of T1-low MST, ratio of T2-high MST, ratio of T1-high MST, disorder of T1-low MST, and standard deviation of T1-high for MST and LPE1, SPE4, LPE4, LPE10, and SPE10 for GRLM.
For classification of immune signature status, the accuracies from the random forest were 72.0% and 77.3% for IE and IS in MST, respectively, and 66.3% and 80.0% for IE and IS in GRLM, respectively.The AUC was 0.747 and 0.743 for IE and IS in MST, respectively, and 0.681 and 0.789 for IE and IS in GRLM, respectively.95% confidence intervals (CI) of ROC are 93.7-96.4% for IS in MST and 88.5-92.3%for IS in GRLM, respectively.Table 3 summarized the true positive rate, true negative rate, classification accuracy, and 95% CI for MST and GRLM, and Fig. 6 shows the ROC curves for the classifiers for MST and GRLM, respectively.The optimal cutoff points were determined by maximizing the sum of sensitivity and specificity.

Discussion
In this study, we identified four distinct groups of voxel intensities (habitats) within the tumor ROI obtained from different MR sequences and separated them into high and low intensities using Gaussian mixture model.All studies were performed in a 2D slice with the maximum tumor area in both T1 post-contrast image and the T2 FLAIR image.These four habitats are considered as distinct entities within an ROI and their spatial relationships are characterized using MST and GRLM approaches.The MST is one of the most common characterization of the spatial proximity of points distributed topologically in space 37 .Gray-level run-length matrix is also  www.nature.com/scientificreports/resulting in differences in voxel resolution (256 × 256 or 512 × 512), slice thickness (1.4-5.0 mm), repetition time (4.9-3285.6 ms for T1 and 400-1100 ms for T2), and echo time (2.1-20 ms for T1 and 14-155 ms for T2).In this study, we performed image preprocessing steps such as voxel isotropic resampling and intensity normalization to make the MR image aspects comparable across various patients.However, these variations in MR images need to be examined more systematically with both MST-and GRLM-based features.Second, we focused on the slice with the maximum tumor area, which was due to a couple of factors such as balancing between accuracy, efficiency, and clinical relevance.As a preliminary study, we intend to pinpoint the region that potentially has the most significant impact on diagnosis and treatment decisions.However, the analysis with the entire tumor volume should be explored in the future study.Next, in this study, we performed a binary classification approach rather than a prediction of survival time.The reason for choosing the 12-month survival cutoff rather than utilizing a more continuous prediction of survival time is a matter that warrants clarification within the clinical context.Also, binary classification simplifies the outcome into distinct categories, which can be easier to interpret and communicate to both medical professionals and patients.However, using continuous prediction instead of binary classification in data analysis can offer several advantages.For example, continuous prediction provides a more detailed and nuanced understanding of outcomes.Also, continuous prediction retains the full range of survival times, preserving information that might be lost in binary categorization.Also, it can provide individualized predictions of survival time for each patient.This personalized information can guide treatment decisions, support patient counseling, and help tailor interventions to the unique circumstances of each individual.However, the choice should align with the research objectives and the utility of the outcomes in the given context.Since continuous prediction provides several advantages and valuable information, it needs to be examined in future studies.Lastly, for predicting immune signature status, we identified 34 patients based on the availability of information regarding immune signature status.Although we used cross-validation as the primary validation method, which is a valid and often preferred approach, especially when the dataset is limited in size, the future study should be explored with a larger dataset.Glioblastoma is a highly aggressive and malignant type of brain tumor and researchers have been increasingly studying the differences in disease patterns and progression for brain cancers.Recently, several studies investigated sex-specific difference in progression and impact on survival for the glioblastoma patients 38,39 .They suggested that gaining a better understanding of the implications and interplay of sex differences in brain cancers will be informative for clinical practice and biological research, which could be considered in our future study with our proposed model.
In this study, we presented graph-based methods for characterizing the spatial proximity of radiologicallydefined habitats in glioblastoma tumors, using a minimal spanning tree and graph run-length matrix construction.According to our results, the features from both MST-and GRLM-based methods provided quantitative metrics of image heterogeneity that have prognostic value for patient survival with high accuracy for both methods.We surmise that the spatial proximity features of habitats based on the MST and GRLM approaches may offer a promising method as a clinical prognostic tool in glioblastomas.However, this approach needs to be validated further in an independent patient cohort to confirm its predictive potential.

Figure 1 .
Figure 1.(a) T1 post-contrast image, (b) T2 FLAIR image with arrows point to the enhanced tumor areas, and (c) the segmentation results on T2 FLAIR image using MITK.

Figure 2 .
Figure 2. MSTs of the four habitats within the tumor ROI on an overlapped T1 and T2 ROI map.(a) T1 postcontrast high voxel intensity group, (b) T1 post-contrast low voxel intensity group, (c) T2 FLAIR high intensity voxel group, and (d) T2 FLAIR low intensity voxel group, respectively.Different gray levels within the ROI represent different habitats and their overlapping areas (scale 0-10).

Figure 3 .
Figure 3.An example of ROI map with a Delaunay triangulation across all vertices.Different gray levels within the ROI represent different habitats and their overlapping areas.

2 Figure 4 .
Figure 4. (a) Illustration of a single initial node located at the center and (b) a graph run-length matrix for this single initial node.This method was adapted from Tosun et al.

Figure 5 .
Figure 5. ROC curve for prediction of 12-month survival status.The x-axis is the false positive rate (or 1 -specificity); the y-axis is the true positive rate (or sensitivity).The area under the ROC curve is 0.832 for MST and 0.773 for GRLM.The optimal cut off points are (0.23, 0.73) and (0.21, 0.82) for GRLM and MST, respectively.

Table 1 .
Patient demographics.The data for age and overall survival are means with standard deviations.OS overall survival.

Table 2 .
The results of the ROC analysis for the 12-month survival prediction.AUC area under ROC curve, MAE mean absolute error, ACC accuracy.