Dosiomic predictors of biochemical failure in patients with localized prostate cancer treated with Iodine-125 low-dose-rate brachytherapy

Abstract

Background

This study aimed to identify dosiomic features that have a significant impact on biochemical failure (BCF) following low-dose rate (LDR) brachytherapy treatment using Iodine-125 seeds for prostate cancer and to provide insights into LDR brachytherapy treatment efficacy using a dosiomic approach.

Methods

Between January 2005 and February 2015, 1,205 patients with localized prostate cancer underwent Iodine-125 seed implantation without combined external irradiation. A total of 96 patients were selected for this study, including 48 with BCF and 48 without BCF. The patients were divided into two cohorts: derivation and validation. Dose distribution images (DDs) were calculated from computed tomography (CT) images taken one month after implantation. A total of 1,130 dosiomic features, including shape-and-size, histogram, and texture features, were extracted from these DDs, their wavelet-transformed images, and Laplacian-of-Gaussian (LoG)-filtered images. The features obtained were categorized into three groups: shape-and-size (S), histogram (H), and texture (T). The Boruta algorithm was used to eliminate less important features. Two analyses were performed: Analysis A performed a multivariate logistic regression analysis using data from the validation cohort to identify significant features. Analysis B generated logistic regression models using derivation cohort data. The accuracy of BCF prediction was assessed using the validation cohort, with performance measured using the area under the receiver operating characteristic curve (AUC).

Results

After the feature reduction process, two, two, and four features remained in the S, H, and T feature groups, respectively. In analysis A, the multivariate logistic regression identified four dominant features, two from each of the S and T groups. In analysis B, the AUC of the logistic regression prediction models using S, H, and all four features were 0.81, 0.77, and 0.86, respectively.

Conclusions

Four significant dosiomic features were identified. Notably, three features—elongation, Maximum2DDiameterRow, and wavelet-HHL_Skewness—strongly distinguished patients with favorable prognoses from others. These findings suggest that dosiomic features from postimplant CT and dose distribution may serve as effective factors for evaluating LDR brachytherapy outcomes in patients with prostate cancer.

Introduction

Low-dose-rate (LDR) brachytherapy for localized prostate cancer has a long history, beginning in 1917 with the transperineal implantation of radium needles by Barringer, and is even today located in the standard treatment option with good clinical performance for low- and intermediate-risk localized prostate cancer [1,2,3,4]. LDR brachytherapy using Iodine-125 seeds showed an outcome equivalent to that of high-dose-rate [HDR) brachytherapy for low-, intermediate-, and selected high-risk patients and comparable oncological outcomes with radical prostatectomy [5,6,7].

Recently, radiomics has been increasingly gaining importance in cancer research, as it extracts a large number of quantitative imaging features from medical images such as computed tomography (CT), magnetic resonance imaging (MRI), and positron emission tomography (PET), contributing to improving diagnostic, prognostic, and predictive accuracy. Radiomics is expected to become a driving force of data-driven decision-making in cancer treatment strategies. [8,9,10,11] Regarding prostate cancer, CT-based and MR-based radiomics studies have demonstrated their potential to classify risk groups of the patients and prostate lesion aggressiveness based on Gleason score. [12, 13] Dosiomics, or dose-based radiomics, is the application of a radiomics strategy to dose distributions, extracting quantitative features from these dose distributions for use in prognostic prediction models [14, 15]. Previous studies have demonstrated that dosiomic prediction models outperform those with conventional dose indices. Regarding prostate cancer, the previous research indicated that the dosiomic features extracted from the clinical target volume of external beam radiotherapy were significantly correlated with biochemical failure (BCF) [14]. These reports demonstrated its powerful prediction performance; however, its areas of application were limited to the dose distribution of external beam radiotherapy.

Prior studies have indicated that CT-based postimplant dose evaluation is essential and strongly recommended [4, 16]. This study aimed to identify dosiomic features regarding CT-based dose evaluation at approximately 30 days postimplantation. These features were assessed for their impact on BCF following LDR brachytherapy using Iodine-125 seeds and providing insights into LDR brachytherapy treatment efficacy through a dosiomic approach.

Materials and methods

Patients

Between January 2005 and February 2015, a total of 1,205 patients with localized prostate cancer underwent Iodine-125 seed implantation without combined external irradiation. The inclusion criteria of the treatment were basically low- and intermediate-risk prostate cancer according to the risk classification in the National Comprehensive Cancer Network (NCCN) 2019 guidelines. The medical records of 133 patients were randomly selected and reviewed. There were 85 and 48 patients with non-BCF and BCF, respectively. BCF was defined according to the Phoenix ASTRO consensus [17, 18]. To eliminate the inclination of patients’ clinical background, propensity score matching with the nearest matching method using age, prostate-specific antigen at diagnosis, clinical T-stage, ISUP grade group, and follow-up duration was applied, and the number of patients with non-BCF and BCF included in this study were 48 and 48, respectively [19]. Details of the patient characteristics are shown in Table 1. A total of 96 patients were randomly divided into two cohorts: derivation (48) and validation (48). Both cohorts included 24 patients with BCF and 24 patients with non-BCF.

Table 1 Patient characteristics

Treatment planning, CT images, and dose distributions

Treatment planning was performed concurrently with real-time ultrasound-guided Iodine-125 seed implantation using the modified peripheral loading method, with the intention of achieving a minimum dose of 145 Gy and a dose ranging from 195 to 200 Gy covering 90% of the prostate volume (D90% = 195–200 Gy) [20]. A postimplant CT scan with a 2-mm slice thickness was acquired approximately 30 days after implantation, and dose distributions (DDs) were calculated on the CT scan using Variseed software version 9.0.2 (Varian Medical Systems, Palo Alto, CA, USA). The field-of-view of CT images varied from 100 to 200 mm², but the calculation grid size for all DDs was uniformly set to 2.0 × 2.0 × 2.0 mm³ due to the specification of Variseed software.

Dosiomic feature extraction

Two experienced radiation oncologists contoured a region-of-interest (ROI) of the prostate gland using Variseed software, and the contoured ROIs were checked by the radiation oncologists against each other. The extracted dosiomic features within the ROI of the prostate gland included 14 shape-and-size features (S), 18 histogram (first-order) features (H), and 75 texture features (T). The S features were the ones that reflect the information of the geometrical shape and size of the prostate gland. Here, the S features were not excluded in this study, although they were not directly related to DDs. The H and T features reflect the information of differential dose-volume histograms and the textural patterns of DDs within the ROI, respectively. Dosiomic features can be determined not only on the original DDs but also on those to which image filters were applied. Image filters may enhance the characteristics that are difficult to find in the original DD images and help find meaningful features [21]. In this study, S, H, and T features were extracted from the original DDs, and H and T features were extracted from eight types of wavelet-transformed and Laplacian of Gaussian (LoG)-filtered DDs (Fig. 1) using PyRadiomics Ver. 3.0.1 [22]. Wavelet transform was applied by three-dimensional discrete wavelet transform with a total of eight filters, and its frequency components were HHH, HHL, HLH, HLL, LHH, LHL, LLH, and LLL, where “H” and “L” denote high-pass and low-pass filters, respectively, and the order of three capital letters represents x-, y- and z-direction, respectively [23]. The Laplacians of Gaussian filters (sigma = 1.0, 1.5, and 2.0) were also applied to the original DDs. Consequently, 1,130 dosiomic feature values were obtained (Fig. 2) and normalized using the z-score standardization method [24, 25].

Fig. 1

Dose distribution, wavelet transformed sub-bands, and Laplacian of Gaussian-filtered images to be analyzed. The area circled by the red line indicates the prostate ROI

Fig. 2

Analysis workflow of the study

Study design

This study contained two analyses: (A) searching for significant dosiomic features to predict BCF and (B) evaluating the significance of dosiomic features found in analysis A using BCF/non-BCF prediction with a logistic regression model and log-rank test. The workflow of this study is shown in Fig. 3.

Fig. 3

Dosiomic feature groups and feature reduction using Boruta

Feature reduction

Before analysis, a feature-reduction process was implemented for each feature group in the derivation cohort using the Boruta algorithm [26]. Boruta is a novel feature selection algorithm that identifies all potentially relevant features using an importance measure in a Random Forest classification algorithm. It operates by creating shadow features through duplication and random shuffling of the original feature set and then training a random forest classifier using both the original and shadow features. The algorithm evaluates the feature’s importance by comparing the performance of the original features to that of the best-performing shadow feature. Through an iterative process, features significantly outperforming shadow features are deemed “important,” while those underperforming are considered “unimportant.” This process continues until all features are classified or the maximum number of runs is reached. Boruta employs statistical tests to validate the feature’s importance and reduce false positives. The Boruta package running in Python (boruta_py) was used in the current study [27].

A total of 1,130 features were extracted from the DDs (14 S, 18 H, and 75 T), wavelet-transformed DDs (144 H and 600 T), and LoG-filtered DDs (54 H and 225 T). These features were sorted into three groups: S, H, and T. The three groups were individually processed using the Boruta algorithm to eliminate less important features (Fig. 3). The remaining features after the Boruta feature-reduction process were used in analyses A and B.

Analysis A: significant dosiomic features

Multivariate logistic regression analyses were performed using the remaining features after Boruta to identify significant features in each feature group (S, H, and T) in the validation cohort. The p-values were adjusted using the Benjamini and Hochberg method [28].

Analysis B: biochemical failure prediction using logistic regression model

To evaluate the significance of the dosiomic features found in Analysis A, logistic regression models were fitted to the remaining features from each of the S, H, and T feature groups in the derivation cohort, and the prediction performance of the models was evaluated using data from the validation cohort. The prediction performances were compared using the area under the curve (AUC) values from the receiver operator characteristic (ROC) analysis. An evaluation of whether these features could distinguish between the BCF and non-BCF groups was also performed using the log-rank test.

Statistical analysis

The clinical backgrounds of patients with BCF and non-BCF were compared using the Mann-Whitney U test, performed using R software, version 4.2.2 (R Project for Statistical Computing). Propensity score matching between patients with BCF and non-BCF, log-rank tests, and Kaplan-Meier curves were calculated using the R software. Multivariate logistic regression analysis, prediction model generation, and validation were performed using the Statsmodels module (version 0.14.0), running on Python version 3.11.4.

Results

Feature reduction

The extracted 1,130 features were sorted into three groups: S, H, and T. Features in the three groups were individually processed using the Boruta algorithm to eliminate the less important features. Table 2 lists the dosiomic features left after the Boruta feature reduction process. The remaining features of each feature group were two from S, two from H, and four from T.

Table 2 Remaining dosiomic features after Boruta feature reduction process, and their p-values in analysis A

Analysis A

Multivariate logistic regression analyses were performed using the remaining dosiomic features in each group after Boruta feature reduction. Table 2 shows the p-values; the asterisks indicate significant differences. The features that showed significant differences were Elongation and Maximum2DDiameterRow from the S group and wavelet-HHL_Skewness and LoG_Sigma_1.5 mm_Energy from the H group, but none of the significant features were found from the T group. Figure 4 (a)–(d) show boxplots of the four features, with the vertical axes in their original scales without z-scoring normalization. Table 2 also includes the mean values and standard deviation (SD) of these features in their original scales.

Fig. 4

Boxplots of significant features found in analysis A and related features with their own scale: (a) original_shape_Maximum2DDiameterRow, (b) original_shape_Elongation, (c) wavelet_HHL_histogram_Skewness, (d) log_sigma_1.5mm_3D_histogram_Energy, (e) wavelet_HHL_histogram_Mean, (f) original_shape_MeshVolume, (g) number_of_seeds

Analysis B

Three logistic regression models were generated from the remaining dosiomic features: two features from S, two features from H, and all four features in the derivation cohort, and their distinguishing performance was evaluated using the same features in the validation cohort. The resulting ROC curves and AUC are shown in Fig. 5 (a)–(c). The obtained AUC values were 0.81, 0.77, and 0.86 for the models generated by two from S, two from H, and all four features, respectively. The logistic regression model that used these four features yielded the highest AUC values.

Fig. 5

Receiver operating characteristics curves and area-under-the-curve values of biochemical failure prediction by logistic regression models using features found in analysis A

The log-rank test was performed to evaluate the accuracy of the identified dosiomic features. Figure 6 (a)–(d) show the Kaplan-Meier curves of patient groups separated by the median of each dosiomic feature found and p-values of the log-rank tests. The obtained p-values are 0.003, 0.7 × 10^− 8, 0.004, and 0.1, from Elongation (S), Maximum2DDiameterRow (S), wavelet-HHL_Skewness (H) and LoG_Sigma_1.5 mm_Energy (H), respectively. The first three features significantly distinguished between the two groups.

Fig. 6

Kaplan-Meier curves of two patient groups separated by the median of the selected features. The p-values were resulting from the log-rank test

Discussion

LDR brachytherapy is a very efficient and cost-effective treatment, although it is also a high-complexity procedure according to the complexity index of interventional radiotherapy index (COMIRI) classification [29, 30]. The results of the present study indicate that dosiomic analysis might be an effective tool for CT-based postimplant dose evaluation.

In this study, two significant features from each of the S and H groups were found in analysis A, and three features—Elongation and Maximum2DDiameterRow from the S group and wavelet-HHL_Skewness from the H group—were finally identified as well-differentiating features.

Features from the S group primarily represented the geometrical shape and volume of the treatment target, that is, the patient’s prostate, rather than the dose distribution. Elongation, which is defined as the square root of the ratio between the major and second major axis lengths of the prostate, characterizes the elliptical nature of its shape and indicates how closely it resembles a spherical or oval sphere [11]. The mean value and standard deviation of elongation were 0.8368 ± 0.0788 and 0.8922 ± 0.0659 for BCF and non-BCF groups, respectively, suggesting that the non-BCF group’s prostates were more spherical and a more spheroidal shape may be advantageous.

Maximum2DDiameterRow is defined as the largest pairwise distance between any two points in the sagittal slices; therefore, it directly represents the size of the prostate itself. The mean value and SD of Maximum2DDiameterRow were 39.13 ± 4.537 and 44.78 ± 5.565 for BCF and non-BCF groups, respectively. The prostate size in confirmed patients with BCF was smaller than that in confirmed patients with non-BCF. Not only elongation but Maximum2DDiameterRow also strongly distinguished the benign prognosis cohort from the others, as shown in Fig. 6 (a) and (b). McNeely et al. indicated that patients with smaller prostates were at higher risk of having a lower D90% (dose covering 90% of the prostate volume) than those with larger prostates, which might cause an increase in BCF, and it is an assertion that might overlap with the results of the current study [31].

One significant and one non-significant feature was found in the H group: skewness from wavelet-HHL images and energy from LoG-filtered images with 1.5 mm of sigma. Skewness measures the asymmetry of a histogram; a higher value indicates that the tail of the histogram is inclined toward the side with high values.

$$\:skewness=\frac{{\mu\:}_{3}}{{\sigma\:}^{3}}=\frac{\frac{1}{{N}_{p}}{\sum\:}_{i=1}^{{N}_{p}}{\left(\varvec{X}\left(i\right)-\stackrel{-}{X}\right)}^{3}}{{\left(\sqrt{\frac{1}{{N}_{p}}{\sum\:}_{i=1}^{{N}_{p}}{\left(\varvec{X}\left(i\right)-\stackrel{-}{X}\right)}^{2}}\right)}^{3}}$$

(1)

where $\:{\mu\:}_{3}$ is the 3rd central moment, $\:\sigma\:$ is standard deviation, $\:{N}_{p}$ is the number of voxels in the ROI, $\:\varvec{X}\left(i\right)$ is a voxel value of i^th voxel in the ROI, and $\bar X$ is an average value of $\:\varvec{X}\left(i\right)$ in the ROI. The wavelet-HHL sub-band is derived by applying high-pass filtering in the x-y plane (axial slices) and low-pass filtering in the z-direction; typical examples of both patients with BCF and non-BCF are shown in Fig. 7. While the mean values of the wavelet-HHL voxel distributions were similar between the BCF and non-BCF groups, as shown in Fig. 4 (e), the difference in skewness indicated a variation in tail elongation in the higher value region. This trend can be inferred as follows: the wavelet-HHL sub-band enhances the voxel values around the iodine sources, and patients with non-BCF, who have a larger prostate volume and need more iodine sources, have more voxels with higher values in wavelet-HHL, as shown in Fig. 4 (f), (g), and Fig. 7.

Fig. 7

Sample images of wavelet-HHL sub-bands, (a) patient with non-BCR and (b) patient with BCR. Window width (WW) and level (WL) of those images were 100,000 / 0 for images in the upper row and 1000 / 40,000 for the ones in the lower row

Energy is the sum of the squared doses of each voxel, that is, a larger irradiated dose causes a higher energy value.

$$\:energy=\sum\:_{i=1}^{{N}_{p}}{\varvec{X}\left(i\right)}^{2}$$

(2)

The mean value and SD of LoG-1.5 mm_Energy were 3.374 × 10¹³ ± 1.036 × 10¹³ and 3.755 × 10¹³ ± 1.451 × 10¹³ for BCF and non-BCF groups, respectively. This feature is understandable because patients with non-BCF showed higher energy values. However, this feature did not separate the BCF and non-BCF cohorts very well, as shown in Fig. 6 (d) (p = 0.1).

Consequently, three features out of four significant features effectively distinguished patients with benign prognoses from others (p < 0.05). The present study indicated that S and H features distinguished patients with benign prognosis from others after LDR brachytherapy, in contrast to the earlier study by Murakami et al. on dosiomic analysis of external beam therapy [14].

The limitations of this study include the relatively small number of patients involved in the analysis and its retrospective nature. Multiple institutional studies with large amounts of patient data may strengthen the findings of the present study.

Conclusion

The present study identified two shape-and-size and two histogram features as significant prognostic factors using a dosiomic analysis approach. Three features, elongation and Maximum2DDiameterRow from shape-and-size feature group and wavelet-HHL_Skewness from histogram feature group, out of four strongly distinguished patients with benign prognosis from the others. This study demonstrates that dosiomic analysis using postimplant CT and dose distribution may serve as effective factors for postimplant evaluation of LDR brachytherapy outcomes in patients with localized prostate cancer.