Prediction of Pituitary Adenoma Volumetric Response to Gamma Knife Radiosurgery Using Machine Learning-Supported MRI Radiomics

1. Introduction

Pituitary adenomas (PAs) are predominantly benign brain tumors, constituting 10%-20% of all intracranial neoplasms, and account for most sellar and parasellar tumors. Gamma Knife radiosurgery (GKRS) has emerged as an effective and minimally invasive treatment modality with high therapeutic efficacy for various intracranial pathologies, including non-functioning and functioning PAs. Radiosurgery (SRS), involving hypofractionated SRS delivered in up to five fractions, achieves tumor control rates ranging from 80% to over 90% in published reports [1,2,3]. Although GKRS achieves favorable outcomes, the individual treatment response varies substantially among patients, necessitating more reliable predictors of individual treatment outcomes [4,5,6].

A common approach in PA treatment is the transsphenoidal removal of the central part of the PA, frequently leaving the remnants within the cavernous sinus for GKRS treatment. Volume response, commonly assessed through post-treatment imaging, is an essential indicator of GKRS efficacy and can inform subsequent management decisions. Identifying reliable predictors of volume response would enable early intervention for poor- and non-responders, potentially improving overall patient outcomes. If reliable pre-treatment outcome prediction existed, surgical treatment could be avoided in selected PAs predicted to respond significantly to GKRS. Conversely, in cases with a less favorable prognosis to GKRS, where higher radiation doses may be required, upfront surgical removal of the central parts of the PA can be preferable to avoid radiation injury to adjacent critical structures, particularly the optic nerves, chiasm, and tracts.

Radiomics focuses on extracting and analyzing intensity, morphological and texture features from medical images. These features capture complex tumor characteristics that are often invisible to the human eye, enabling the potential identification of imaging biomarkers for treatment response prediction. Radiomics has shown promise in various oncological contexts, allowing for non-invasive prediction of treatment response or prognosis of disease outcome [7,8].

We hypothesized that radiomics prediction signatures could outperform the benchmark signature relying only on routine clinicopathological (CP) features for predicting PA volumetric response.

The objective of this study was to develop several single-modality and multi-modality radiomic prediction signatures, both with and without CP parameters, to predict the volumetric response of PAs to GKRS, using pre-treatment MRI. Additionally, the study aimed to compare the predictive performance of these signatures in the test folds. Improving the prediction of treatment outcomes is highly clinically relevant because it enables personalized treatment planning, optimizes therapeutic strategies and thereby enhances patient care.

2. Methods

This study adheres to the guidelines set out in the STROBE statement for cohort studies.

2.1. Ethics Approval Statement

The study received approval (BLINDED FOR REVIEW) from our Institutional Review Board and adheres to The Code of Ethics of the World Medical Association (Declaration of Helsinki), as published in the British Medical Journal (July 18, 1964) and its 7th revised edition in 2013. The need for written informed consent was waived by the ethics committee due to the retrospective nature of the analysis.

2.2. Patients

Included were 81 patients between 11.4 and 83.2 years of age (mean 45.6) with imaging-diagnosed PA treated at our Gamma Knife center. MRI had been performed within four weeks before their radiosurgery, with available follow-up data after an interval of six months or longer (range 6.7–105.5 months, mean 40.4; Table 1). Tumor volumes ranged between 0.18 and 40.33 cm³ (mean 6.30).

2.3. Sample Size Calculation

The prospective sample size calculation was based on a pilot experiment with 30 patients and required 36 patients for alpha=0.05, beta=0.20, and AUC=0.76 (Medcalc 14.8.1; MedCalc Software Ltd., Ostend, Belgium). These 30 patients were included in the final cohort. The actually obtained AUCs for the calculated scores integrating radiomics and CP features ranged between 0.759 and 0.928, with a final sample size of 81 patients in T1w, 81 in CE-T1w, 48 in T2w, and 41 in FLAIR.

2.4. Gamma Knife Treatment

The Gamma Knife technique was previously described.⁷ On the day of the treatment, after placing a stereotactic G frame (Elekta AB), under sedation and local anesthesia, contrast-enhanced 3D computed tomography (CT) imaging was obtained. Pretreatment MRI sequences, acquired less than four weeks before, were then coregistered to the stereotactic CT. Treatments were planned on a Leksell GammaPlan 10.1 workstation (Elekta Instrument AB, Stockholm, Sweden) carefully respecting dose constraints, particularly of the optic apparatus.

The margin dose varied from 12 to 40 Gy (mean 20.5 Gy), depending on tumor size, location, and hormone status (Table 1). Forty-nine PAs were treated in a single session, with margin doses ranging from 12 to 35 Gy (mean 17.8 Gy). According to our institutional protocol, lesions abutting the anterior optic pathway are treated using hypofractionated radiosurgery (HFSRS). Thirty-two PAs were treated with 5.0 to 9.0 Gy (mean 6.1 Gy) for 3 to 5 days (mean 4.03). Biologically effective dose (BED) is routinely used to compare doses of different dose-fraction regimens relying on the widely accepted Linear Quadratic model, with its known limitations for high doses [9] HFSRS doses can be converted to single fraction equivalent doses (SFED) to intuitively compare radiation effects to conventional physical doses.[10] Margin SFED varied from 11.1 – 35.0 Gy (mean 16.2 Gy), applying an α/β ratio of 2.47 Gy for nonfunctional PAs or 4.91 Gy for functioning Pas [11,12].

2.5. MRI

MRI was performed on a 3-Tesla scanner (Achieva; Philips, Eindhoven, Netherlands). 3D T1-weighted non-contrast (T1w), contrast-enhanced (CE-T1w), T2w, FLAIR sequences were acquired with the following sequence parameters:

• 3D-T1w magnetization-prepared rapid acquisition (MPRAGE) sequence: gradient echo, TR/TE/TI 6.8/3.2/900ms, flip angle 8°, measured voxel size 0.6*0.6*1.0 mm, before and after intravenous injection of contrast medium.
• T2w sequence: TR/TE 3693.8/80ms, 150 transversal slices, thickness 1mm, matrix 512×512.
• Fluid-attenuated inversion recovery (FLAIR) sequence: TR/TE/TI 11,000/120/2800ms, 90 transversal slices, thickness 2mm, matrix 512×512.

2.6. Postprocessing

In all 81 patients, PA volumes were delineated and measured from CE-T1w images on the Leksell GammaPlan workstation. MRI sequences were coregistered to the stereotactic 3D contrast-enhanced CT acquired on treatment day. Image sets were verified for high image quality. Sequences with artifacts were excluded.

2.7. Follow-up

Imaging and clinical follow-up were performed every six months for the first two years after GKRS and annually thereafter.

2.8. Feature Extraction

The radiomics analysis was conducted using the open-source Pyradiomics plugin, integrated into 3D Slicer (version 5.6.1) [13]. The Pyradiomics parameter file (Supplementary file 1), instructed the computation of all available image transformations and feature types, resulting in the total of 2156 features per MRI scan.

Feature extraction was performed on the original images (107 features) and on images transformed by: wavelet, square, square root, logarithm, gradient, exponential, Laplacian of Gaussian (LoG) and Local Binary Pattern (LBP 2D) filters. For detailed descriptions of the extracted radiomics features please see: https://pyradiomics.readthedocs.io/en/latest/features.html.

2.9. Model Selection

Features were pre-selected by discarding those without a significant correlation with the PA volume change outcome, based on Pearson’s correlation test (Statistica 12.0, StatSoft, Hamburg, Germany). Predictive models were constructed using the features selected by LASSO regression with integrated 10-fold cross-validation and the following classifiers: random forest, naïve Bayes, kNN, logistic regression, neural network and SVM (Orange Data Mining, University of Ljubljana, Slovenia).

2.10. Evaluation of Predictive Performance

The PA volume change per natural logarithm of time was chosen as the endpoint to account for the near-exponential decrease of tumor volume over time. [15] Beyond an initial 6-months phase, our data indicates that most PA volumes follow an exponential trend. To compare volume changes during long FU periods (FUP), tumor volume change per natural logarithm of time is more accurate than volume change per time. This was calculated as the difference in tumor volume before SRS and at the last follow-up, divided by the initial volume and the natural logarithm of the time since SRS to the last FU: [15]

                         Volume at last FUP - Initial volume at SRS

Relative volume change per ln(month) = ---------------------------------------------------------

                                  Initial Volume * ln(FUP)

Predictive performance was assessed using ROC analysis (IBM SPSS v28, Armonk, NY, USA; Orange Data Mining), with statistical analyses considering both continuous independent variables and continuous or categorized dependent outcomes.

2.11. Validation

Validation involved bootstrap for ROC analysis and split-sample cross-validation via LASSO. Besides LASSO cross-validation in the training folds (51 patients), models were validated on eight random test folds (30 patients).

3. Results

3.1. Patients’ Characteristics

The treatment results are presented in Table 1. After a mean FUP of 40.4 months, the control rate was 98.8%. As there are no standardized radiographic criteria for assessing PA treatment response, we used the RANO criteria [17] stated by Imber et al. for PA response characterization.[18] Partial response was achieved in 60 [74.1%] PAs, with stable disease found in 20 [24.7%] of patients, while one [1.2%] progressive PA was observed at last follow-up. The volumetric outcome was notably better in functioning PAs, with a volume reduction of 2.21% per month, compared to a lower volume reduction per month of 1.22% for non-functional PA. This difference in volumetric response is probably attributed to the significantly higher doses used in treating functioning PAs, required to obtain hormone remission.

3.2. Experimental Design

The workflow of the study is shown in Figure 1. Models were developed for individual modalities (CP, T1w, CE-T1w, T2w, FLAIR) and combined modalities (T1w+CE-T1w, CP+T1w+CE-T1w), focusing on T1w and CE-T1w MRI sequences available for all 81 patients. This facilitated multi-sequence model construction and validation across eight test folds. Tumor volume response prediction used R² from linear LASSO (Table 2) for continuous outcomes, and logit LASSO with six classifiers for categorized outcomes, using AUC and accuracy for evaluation (Table 3 and Table 4).

3.3. The Predictive Models for PA Response to Radiosurgery

The constructed predictive models were initially tested against the continuous, uncategorized tumor volume response outcomes (Table 2). Thereby, T2w showed the best association with the outcomes (R² =0.665), followed by the most comprehensive models: CP+T1w+CE-T1w (R² =0.584) and T1w+CE-T1w (R²=0.502). The optimal performance of T2w was notable, but because of its availability in a smaller number of patients, this result was interpreted as preliminary. Our focus was thus directed towards CP+T1w+CE-T1w, identified as the second-best performer (Table 2).

The optimal cutpoint(s) for the comprehensive model, combining CP and radiomics features (CP+T1w+CE-T1w), were identified by testing against six different cutpoints in the outcome (Figure 2). The cutpoint of -0.25% tumor volumetric response per natural logarithm of time was identified as the most suitable for prognosis using the available CP and radiomics features (Figure 2).

After determining the optimal cutpoint for outcome categorization, the predictive performance of individual CP parameters was evaluated using the t-test for continuous and the chi-square test for binary values (Table 3). Age, fraction number, accumulated dose, SFED, and hormone secretion status showed significant associations with the -0.25% outcome (Table 3). Models using only CP or CE-T1w features had the weakest link to the -0.25% outcome, while those with T1w features performed best. Notably, there was no significant difference in predictive performance among T1w, T1w+CE-T1w, and CP+T1w+CE-T1w models, as confirmed by an independent-samples t-test. To reduce the feature burden on LASSO, features were pre-selected based on their Pearson correlation with the outcome cutpoint at -0.25% to lighten the LASSO's feature load, enhancing the models' prognostic performance. These features were refined through L1 LASSO selection, with the top eight used in classification by random forest, naïve Bayes, kNN, logistic regression, neural network, and SVM classifiers (Table 4).Comparing classifier performance using the CP+T1w+CE-T1w model for the -0.25% outcome cutpoint, revealed that logistic regression, neural network, and SVM showed similar classification performance in test folds, with random forest, naïve Bayes, and kNN performing less effectively (Table 4).

Table 5 details selected features and coefficients for the model incorporating CP, T1w, and CE-T1w features, highlighting the inclusion of only one CE-T1w feature, likely due to its inferior predictive performance (Table 5).

4. Discussion

This initial study demonstrates that radiomic signatures built on pre-GKRS MRI data can be used to predict tumor volumetric response to radiosurgery. The computed radiomic models outperformed the benchmark model, which included only clinicopathological features. The T1w and T2w models were the best predictive performers, while the dual sequence T1w+CE-T1w and the merged CP+T1w+T2w models did not provide a statistically significant improvement over the single modality models, as judged by the t-test. All models, except for T2w and FLAIR (due to their low numbers), underwent additional evaluation for generalizability using cross-validation. The results demonstrated that the predictive performance in the training folds was largely retained in unseen test folds, with the remaining AUCs consistently exceeding 0.90.

Radiomics analysis treats MRI as minable data by extracting quantitative computational features. It gains importance as clinical imaging becomes increasingly widespread. To the best of our knowledge, no prior studies have employed texture or radiomics analyses on pituitary adenomas. Our group and others have performed texture analyses on meningiomas, [19] vestibular schwannomas [20,21,22] and brain metastases.[23,24]. In a previous study investigating the use of radiomics features to develop a treatment prognosis following SRS, our group applied radiomics to pre-radiosurgical MRI to predict the long-term outcomes of WHO grade 1 meningiomas after SRS, achieving satisfactory predictive performance, as evidenced by an AUC reaching 0.88.[7] Other groups used radiomics to predict outcome and pseudoprogression of vestibular schwannoma treated with SRS.[25,26] Several studies explored radiomics to predict local control of brain metastases [27,28,29,30] or arteriovenous malformations [31,32,33] after SRS using radiomics.

Of the four sequences investigated, we found that T2w provided the best predictive performance for volumetric PA change (Table 2). This result agrees with our previous findings applying texture analysis to predict volumetric outcomes of SRS in benign meningiomas and vestibular schwannomas. For benign meningiomas, we found that the histogram parameter standard deviation of voxel intensities of T2w images correlated best with volumetric change after SRS.[19] Increased T2w intensity values have been reported to relate to a soft consistency of meningiomas, increased vascularity, cellular atypia, and angioblastic or melanocytic components as well as cystic degeneration and ischemic necrosis.[34,35,36] For vestibular schwannomas, kurtosis of T2w image intensity values predicted progression best, with a sensitivity and specificity of 71% and 78%, while the minimum of the T2w voxel intensity values correlated significantly with the final regression of tumor volume per month.[20]

Several studies investigated the use of MRI to preoperatively assess tumor consistency in PAs, as tumor consistency is a critical factor in surgical planning. Most existing studies support the ability of T2w to predict intraoperative consistency, while T1w has not been shown to offer any predictive value.[37] Hypointensity on T2w likely correlates with firmer tumors, possibly attributable to their increased collagen content and vascularity. In comparison, softer tumors tend to be hyperintense on T2w, which may relate to higher water content and/or cystic components.[38,39,40] Although multiple preoperative PA consistency assessment methods have been studied, none demonstrated sufficient accuracy and reliability in clinical use.[41] More recently, radiomics and machine learning-based models have achieved high precision and good AUC values.[41] Wan et al. developed a radiomics model built on combined T1w/CE-T1w/T2w images, with eleven imaging features exhibiting statistically significant differentiation between soft and hard PAs, providing an excellent performance with an AUC of 0.90.[42]

In functional PAs, in addition to volumetric tumor control, hormone remission is a mandatory treatment objective. As more genes must be silenced to achieve hormone remission, substantially higher PA margin doses are required. Because of higher margin doses, we observed a more favorable volumetric response in functional PAs (mean volume change/month -2.21%), compared to non-functioning PAs (mean volume change/month -1.22%). A preliminary analysis revealed enhanced associations within subdivided PA cohorts, achieving an AUC close to 1. Due to the small sample size of functional PA treatments, we were unable to separately perform a reliable statistical analysis in the subset of functional PA tumors. However, our findings warrant further investigation.

The high dimensionality of radiomics analysis is both its strength and its weakness, as its high dimensionality has been widely criticized. To address this issue, we pre-selected features by excluding those that were not significantly associated with the outcome, resulting in an extensive reduction in dimensionality. Additionally, our study showed that adding CP features to multi-sequence imaging data (CP+T1w+CE-T1w) did not improve predictive performance over single-sequence radiomics models.

5. Limitations

Although the patient group was highly homogeneous and vastly exceeded the sample size requirement, its size still posed a limitation. To enhance the clinical validity of the reported feature association with PA volume response to GKRS, further studies in larger patient groups are warranted. The relatively short follow-up time and the retrospective design of the prognostic models were additional limitations. Moreover, predictive studies necessitate confirming the generalizability of the acquired findings through external validation, while internal prognostic validation in unseen test folds, that were not part of the development cohort, was already carried out within this study. Therefore, further validation in external cohorts is needed to establish the prognostic clinical validity of the obtained predictive models. Additionally, despite the objective nature of the computational analysis technique, the workflow included residual subjectivity during the tumor VOI segmentation.

6. Conclusions

We demonstrated that a radiomics-based model using conventional MR imaging achieved excellent predictive classification and generalization performance, surpassing a model based on clinicopathological parameters. By achieving reliable predictions of the volumetric outcomes of SRS, radiomics might enable individualized treatment strategies, ultimately contributing to improved treatment outcomes for patients undergoing SRS for pituitary adenomas.