An interpretable biparametric MRI habitat radiomics model for predicting bone metastasis in prostate cancer: a dual-center study
Original Article

An interpretable biparametric MRI habitat radiomics model for predicting bone metastasis in prostate cancer: a dual-center study

Juntao Gong1, Feixiang Li1, Yun Sun1, Jinghuan Huang1, Yingying Zhang1, Ou Yang1, Jing Li2, Ze Song3, Kai Ai4, Gang Huang5

1Gansu University of Chinese Medicine, Lanzhou, China; 2Zhangye People’s Hospital Affiliated to Hexi University, Zhangye, China; 3Department of Medical Imaging, School of Medicine, Hexi University, Zhangye, China; 4Philips Healthcare (Suzhou) Co., Ltd., Suzhou, China; 5Department of Radiology, Gansu Provincial Hospital, Lanzhou, China

Contributions: (I) Conception and design: J Gong, F Li, G Huang; (II) Administrative support: K Ai, G Huang; (III) Provision of study materials or patients: J Li, Z Song; (IV) Collection and assembly of data: Y Sun, J Huang, Y Zhang, O Yang; (V) Data analysis and interpretation: J Gong, F Li; (VI) Manuscript writing: All authors; (VII) Final approval of manuscript: All authors.

Correspondence to: Gang Huang, MD. Department of Radiology, Gansu Provincial Hospital, 204 Donggang West Road, Chengguan District, Lanzhou 730000, China. Email: huang_g2024@163.com.

Background: Bone metastasis (BM) critically determines prognosis in prostate cancer (PCa), but its noninvasive prediction remains challenging. Habitat imaging may better capture intratumoral heterogeneity and improve BM risk assessment. This study aimed to develop and validate an interpretable habitat radiomics model based on biparametric magnetic resonance imaging (bp-MRI) for predicting BM in PCa.

Methods: This retrospective dual-center study included 238 PCa patients who underwent preoperative bp-MRI [3.0T, T2‑weighted imaging (T2WI) and apparent diffusion coefficient (ADC) maps]. Patients from the primary center were split into training (n=133) and internal test (n=57) sets; 48 patients from a second center formed an external validation cohort. Tumors were segmented into three habitat subregions via K-means clustering on ADC maps. Radiomic features were extracted from T2WI for each subregion and the whole tumor. Four models (clinical, conventional radiomics, habitat radiomics, and combined) were built using a Gaussian process (GP) classifier. Performance was evaluated using area under the receiver operating characteristic curve (AUC), DeLong test, net reclassification improvement (NRI), and calibration curves. SHapley Additive exPlanations (SHAP) analysis assessed interpretability.

Results: The combined model achieved the highest AUCs: 0.927 (training), 0.841 (internal test), and 0.821 (external validation). The habitat model achieved AUCs of 0.870, 0.838, and 0.760 in the training, internal test, and external validation sets, respectively, outperforming the conventional and clinical models in all three cohorts. In the training and internal test sets, the DeLong test showed that the habitat and combined models had significantly better discrimination than the conventional and clinical models (P<0.05); in the external validation set, the combined model remained significantly superior, whereas the habitat model’s superiority did not reach statistical significance. NRI demonstrated improved risk stratification by the habitat model. SHAP identified key predictive features from habitat subregion 2.

Conclusions: An interpretable habitat radiomics model based on bp-MRI shows potential for predicting BM risk in PCa. The combined model provides favorable performance for clinical decision support. External generalizability requires confirmation in larger, multicenter cohorts.

Keywords: Biparametric magnetic resonance imaging (bp-MRI); bone metastasis (BM); habitat radiomics; prostate cancer (PCa); SHapley Additive exPlanations (SHAP)


Submitted Apr 20, 2026. Accepted for publication Jun 29, 2026. Published online Jul 30, 2026.

doi: 10.21037/tau-2026-0380


Highlight box

Key findings

• The combined biparametric magnetic resonance imaging (bp-MRI) habitat radiomics model achieved areas under the receiver operating characteristic curve (AUCs) of 0.927 (training), 0.841 (internal test), and 0.821 (external validation) for predicting bone metastasis (BM) in prostate cancer (PCa).

• The habitat model alone (AUC 0.760 externally) outperformed conventional radiomics and clinical models, and SHapley Additive exPlanations analysis identified Habitat 2 features (first-order skewness and dependence variance) as the most influential predictors, suggesting a “high-energy, low-heterogeneity” signature associated with BM risk.

What is known and what is new?

• BM predicts poor prognosis in PCa, but noninvasive prediction remains challenging. Conventional whole-tumor radiomics overlooks intratumoral heterogeneity.

• This study presents an interpretable bp-MRI habitat radiomics model that partitions tumors into three functional subregions. The combined model achieved robust BM prediction (AUC 0.821) and identified a “high-energy, low-heterogeneity” signature in Habitat 2 as the key predictive feature.

What is the implication, and what should change now?

• Habitat radiomics captures intratumoral heterogeneity, offering superior BM risk stratification. Prospective multicenter studies are needed to integrate this approach into routine bp-MRI workflows for noninvasive, cost-effective BM risk assessment at initial PCa diagnosis.


Introduction

Prostate cancer (PCa) is one of the most common malignancies in men, with an estimated nearly 1.5 million new cases annually worldwide, and is the second leading cause of cancer-related mortality, accounting for nearly 400,000 deaths per year (1). In China, a substantial number of patients still present with metastatic disease at initial diagnosis, particularly bone metastasis (BM) (2). BM is a common complication in advanced PCa, ultimately occurring in approximately 80% of patients (3), and serves as a key determinant of treatment choices and prognosis (4).

Currently, several imaging modalities are available for assessing BM. The 99mTc-methylene diphosphonate whole-body bone scan is recommended as the preferred initial screening method due to its high sensitivity, with a reported pooled sensitivity of approximately 79% (5,6). Positron emission tomography-computed tomography (PET/CT) offers higher diagnostic accuracy (6), and computed tomography (CT) assesses structural bone destruction. However, these modalities involve radiation exposure, have limited specificity leading to false-positives, and offer restricted quantitative assessment (7). Conventional magnetic resonance imaging (MRI) provides excellent soft-tissue resolution but its performance is highly dependent on subjective experience and yields few reproducible quantitative features (8,9). Whole-body MRI (WB-MRI) has emerged as a powerful, radiation-free alternative for comprehensive metastatic assessment, offering superior soft-tissue contrast for detecting bone marrow involvement. However, WB-MRI protocols are not yet universally standardized or routinely integrated into initial PCa staging in all centers, and they typically extend examination time and cost (10).

Biparametric magnetic resonance imaging (bp-MRI), which includes T2-weighted imaging (T2WI) and diffusion-weighted imaging (DWI), has become a key tool for detecting and grading PCa (11,12). The apparent diffusion coefficient (ADC), a quantitative parameter derived from DWI, enhances diagnostic objectivity and reproducibility (12-14). However, routine prostate bp-MRI protocols are confined to the pelvis and do not systematically assess common distant metastatic sites (15).

Radiomics extracts high-throughput features to capture tumor heterogeneity and has shown promise in predicting PCa BM, with one study reporting an area under the receiver operating characteristic curve (AUC) of 0.898 (16-18). However, conventional radiomics often treats the tumor as a homogeneous entity, overlooking the profound spatial intratumoral heterogeneity that is a known driver of aggressive behavior and treatment outcomes (19,20). To overcome this limitation, a method that captures the full spatial heterogeneity by analyzing distinct biologic subregions is required, offering a more nuanced view of the tumor’s complexity. Recent studies have indicated that combining intratumoral and peri-tumoral features improves prediction models (21), leading to the concept of habitat analysis. This method partitions tumors into subregions with distinct microenvironmental characteristics, providing a refined framework for analyzing tumor heterogeneity (22,23). This approach is predicated on the rationale that spatially distinct subregions may harbor unique biological information which is lost when averaging features across the entire tumor volume. In other malignancies, such as glioblastoma and lung cancer, habitat or subregion-based radiomic analyses have demonstrated superior performance over conventional whole-tumor approaches in tasks like survival prediction and lesion characterization (24,25). Despite advances in other cancers, studies applying habitat imaging to predict BM in PCa remain scarce.

Therefore, this study aimed to develop and validate an interpretable habitat radiomics model based on biparametric MRI to noninvasively predict the risk of BM in PCa patients. We sought to assess whether this habitat-based approach could outperform conventional radiomics and clinical models, and to evaluate its generalizability through a dual-center study design. We present this article in accordance with the TRIPOD reporting checklist (available at https://tau.amegroups.com/article/view/10.21037/tau-2026-0380/rc).


Methods

The overall radiomics analysis followed a standardized workflow encompassing image acquisition, tumor segmentation, preprocessing, habitat subregion clustering, feature extraction, selection, and model development, as visually summarized in Figure 1.

Figure 1 Flowchart of the predictive model construction based on medical imaging habitat subregion analysis and radiomic features. ADC, apparent diffusion coefficient; AUC, area under the receiver operating characteristic curve; DCA, decision curve analysis; LASSO, least absolute shrinkage and selection operator; ROI, region of interest; SHAP, SHapley Additive exPlanations.

Study population and data collection

The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. This study was approved by the Ethics Committee of Gansu Provincial Hospital (approval No. 2025-489) and the Ethics Committee of Zhangye People’s Hospital Affiliated to Hexi University (approval No. HFYER-2023.07). Informed consent was waived in this retrospective study.

Patients were sourced from two institutions: Gansu Provincial Hospital (Center 1, primary center) and Zhangye People’s Hospital Affiliated to Hexi University (Center 2, collaborating center). We consecutively enrolled patients with pathologically confirmed PCa who underwent preoperative bp-MRI between January 2021 and December 2023.

Inclusion criteria were: (I) histopathologically confirmed PCa; (II) available preoperative prostate bp-MRI.

Exclusion criteria included: (I) prior neoadjuvant therapy; (II) incomplete clinical, pathological, or imaging data; (III) inadequate MRI quality for tumor segmentation; (IV) absence of a baseline bone scan.

Confirmation of metastatic status and cohort division

BM status was determined by 99mTc-methylene diphosphonate whole-body bone scan. To enhance the reliability of the non-metastatic cohort, bone scan-negative patients were further confirmed by pelvic or whole-body MRI to exclude early marrow-only metastases.

A total of 238 eligible patients were included. Among them, 190 patients from the primary center were randomly split into a training set (n=133, 70%) and an internal test set (n=57, 30%). An independent external validation set comprised 48 patients from a second center, all meeting the same criteria.

Clinical variables

The following baseline clinical and pathological parameters were collected: age, body mass index (BMI), total prostate-specific antigen (TPSA), free prostate-specific antigen (FPSA), alkaline phosphatase (ALP), number of positive biopsy cores (PNC) and International Society of Urological Pathology (ISUP) grade group. Prostate volume (PV) was calculated from T2-weighted imaging using the ellipsoid formula (width × length × height × 0.52) (26). Prostate-specific antigen density (PSAD) was computed as TPSA / PV.

Image acquisition

All MRI examinations were performed on 3.0T scanners (Siemens) using a phased-array body coil with patients in the supine position. The imaging protocol included axial T1-weighted imaging (T1WI), axial and sagittal T2-weighted imaging (T2WI), and axial DWI, with ADC maps generated automatically, using a mono-exponential decay model. For Center 1, ADC maps were calculated using the b-values of 50 and 1,500 s/mm2. For Center 2, ADC maps were calculated using the b-values of 50 and 800 s/mm2. Detailed acquisition parameters are provided in Appendix 1.

Tumor segmentation and reproducibility assessment

Tumor segmentation was performed independently by two radiologists (each with 10 years of experience) using ITK-SNAP software (version 3.6.0) to delineate the volume of interest (VOI) slice-by-slice. Any discrepancies were resolved by consensus with a third senior radiologist (with over 20 years of experience) to ensure accuracy.

To evaluate the reproducibility of the manual segmentations, a rigorous consistency analysis was conducted. A second radiologist (Reader B) independently segmented tumors in a randomly selected 30% subset (n=57) of cases, blinded to the initial results. The first reader (Reader A) repeated the segmentation on the same subset after a 4-week interval to minimize recall bias. Inter- and intra-reader agreements were quantified using the Dice similarity coefficient (DSC) and the intraclass correlation coefficient (ICC, based on a two-way random-effects model for absolute agreement).

Image preprocessing

The MRI data underwent preprocessing to standardize images and ensure model stability. To ensure spatial consistency across centers, all T2WI and ADC images were resampled to an isotropic voxel size of 1.0×1.0×1.0 mm3. The steps included: (I) intensity inhomogeneity correction of T2WI maps using the N4 bias field correction algorithm; (II) spatial registration of sequences with Advanced Normalization Tools (ANTs) to align ADC with T2WI; (III) outlier removal by cropping data to the 0.01–0.99 quantile range; and (IV) intensity normalization to [0, 1] via zero-centering and unit variance scaling.

Subregion clustering

In the patient-population level, the entire tumor volume was segmented into distinct habitat subregions using an unsupervised K-means clustering algorithm based on grayscale values from the ADC maps. The number of clusters (k) was tested from 2 to 10, with the optimal value (k=3) determined solely on the training cohort based on the highest Calinski-Harabasz (CH) score (Figure S1). All tumor voxels from the training cohort were pooled, and K-means was applied to define global cluster centroids. For each patient, voxels were then assigned to the nearest centroid. To ensure consistent labeling across patients, the three resulting clusters were sorted by their mean ADC value (lowest to highest) and labeled as Habitat 1 (lowest ADC, presumably representing high-cellularity tumor core), Habitat 2 (intermediate ADC), and Habitat 3 (highest ADC, potentially reflecting necrosis or stroma). A three-cluster solution was chosen because it provided the best separation of ADC-defined microenvironments (low, intermediate, and high cellularity), which is biologically plausible for PCa. No spatial constraints (e.g., connectivity or adjacency) were applied during clustering; voxels were treated independently based solely on their ADC values. Cross-center harmonization was addressed through image preprocessing (isotropic resampling, bias field correction, and intensity normalization) as described above, but no additional ComBat or batch-effect adjustment was performed. Further methodological details are provided in Appendix 1.

Representative cases of the study cohort are illustrated in Figure 2, which compares the bp-MRI findings and habitat clustering results between a patient without BM and one with BM.

Figure 2 Representative MRI images and habitat subregion analysis results from two prostate cancer patients. (A-D) Bone metastasis-negative patient: (A) T2-weighted image demonstrates a focal hypointense lesion within the prostate; (B) ADC map shows marked hypointensity; (C,D) Habitat clustering analysis reveals low intratumoral heterogeneity, dominated by a single predominant subregion. (E-H) Bone metastasis-positive patient: (E) T2-weighted image displays an irregular hypointense lesion in the prostate; (F) ADC map exhibits significant hypointensity, suggesting elevated cellular density; (G,H) Habitat clustering analysis indicates high intratumoral heterogeneity, comprising three distinct subregions (blue, green, and red). ADC, apparent diffusion coefficient; MRI, magnetic resonance imaging.

Habitat and conventional radiomics feature extraction

Low-order radiomic features were extracted from both the clustered subregions and the original whole tumor using the PyRadiomics package (version 3.0.1) in Python (version 3.9.13) (available at: https://pyradiomics.readthedocs.io/en/latest/index.html). These features were categorized into 3 groups: (I) shape; (II) first-order statistics; and (III) texture. Habitat radiomic features were extracted from each subregion, with missing values imputed using the median. Features from the ADC map of habitat subregion 1 were labeled as ADC_featurename_h1, whereas those from the T2WI image were labeled featurename_h1. The same naming convention was applied to subsequent subregions. Non-clustered conventional radiomic features from the whole tumor were prefixed with intra_.

Feature selection and model construction

All preprocessing and modeling steps were strictly performed on the training set and applied to the test sets. Radiomic features were standardized using Z-score normalization with mean and standard deviation derived from the training set prior to selection. Feature selection followed a three-step sequential process: initial screening with SelectKBest, redundancy reduction via minimum redundancy maximum relevance (mRMR), and final selection using least absolute shrinkage and selection operator (LASSO) regression (alpha =0.05). Clinical variables associated with BM, identified through univariate and multivariate logistic regression, formed the clinical model.

To prevent data leakage, the normalization parameters (mean and standard deviation) and the final selected feature subsets (from SelectBest, mRMR, and LASSO) were derived exclusively from the training set and then frozen. These frozen parameters and feature sets were subsequently applied to the internal test set and external validation set without any refitting or further feature selection.

Using the selected habitat, whole-tumor, and clinical features, classification models were developed on the training cohort. Ten classical machine learning algorithms were evaluated, such as random forest, support vector machine (SVM), logistic regression, and Gaussian process (GP), among others, using fivefold cross-validation. The optimal model was subsequently validated on the test set.

Statistical analysis

Statistical analyses were conducted using Python (v3.9.13), R (v4.2.2), and SPSS 27. Normally distributed data are presented as mean ± standard deviation, non-normally distributed data as median [interquartile range (IQR)], and categorical data as counts (percentages). Univariate and multivariate logistic regression were used to identify clinical predictors of BM (P<0.05 significant).

Model performance was evaluated using receiver operating characteristic (ROC) analysis (AUC, accuracy, sensitivity, specificity), decision curve analysis (DCA), calibration curves, and the DeLong test. All models were trained on the training set (n=133), then validated on an internal test set (n=57) and an external validation set (n=48). SHAP analysis was applied to interpret feature contributions.


Results

Clinical characteristics

The external validation cohort consisted of 48 patients, of whom 16 (33.3%) had BM. This prevalence provided a clinically relevant scenario for testing the model’s discriminative ability. Patient baseline characteristics are summarized in Table 1. With the exception of the PNC, ALP, TPSA, PSAD, ISUP and FPSA, no significant differences were observed in other clinical parameters (including age, BMI, volume and DRE) between patients with and without BM (P>0.05). The clinical model, constructed using PNC and FPSA identified as independent predictors by multivariate analysis (Table 2), achieved AUCs of 0.787 [95% confidence interval (CI): 0.711–0.863] in the training set and 0.704 (95% CI: 0.566–0.842) in the test set.

Table 1

Clinical and pathological characteristics of 190 PCa patients in different cohorts

Variables Train vs. Test Non-BM vs. BM
Train (n=133) Test (n=57) P Non-BM (n=108) BM (n=82) P
Age (years) 74.00 (67.00, 77.00) 72.00 (68.00, 80.00) 0.57 72.50 (68.00, 77.25) 74.00 (68.00, 78.00) 0.52
BMI (kg/m2) 23.17 (21.11, 25.46) 21.97 (19.94, 24.34) 0.02 23.046 (21.14, 25.16) 22.32 (20.58, 24.80) 0.52
ALP (U/L) 87.00 (68.00, 144.00) 90.00 (69.56, 166.68) 0.49 76.00 (64.00, 97.13) 132.52 (86.61, 225.75) <0.01
PNC (cores) 10.00 (6.00, 12.00) 10.00 (6.00, 12.00) 0.34 8.00 (4.00, 11.00) 11.50 (9.25, 12.00) <0.01
TPSA (ng/mL) 72.89 (22.60, 100.00) 74.31 (37.97, 100.00) 0.65 37.77 (15.30, 77.85) 100.00 (87.39, 100.00) <0.01
Volume (cm3) 43.11 (31.72, 61.78) 42.43 (30.02, 65.52) 0.96 43.75 (30.90, 62.42) 40.94 (30.58, 61.23) 0.39
PSAD (ng/mL/cm3) 1.23 (0.58, 2.10) 1.28 (0.63, 2.34) 0.42 0.80 (0.36, 1.60) 1.84 (1.17, 2.77) <0.01
FPSA (ng/mL) 8.93 (2.84, 30.00) 12.44 (3.79, 30.00) 0.32 4.46 (1.92, 11.69) 30.00 (14.26, 30.00) <0.01
DRE 0.63 0.85
   0 82 (61.65) 33 (57.89) 66 (61.11) 49 (59.76)
   1 51 (38.35) 24 (42.11) 42 (38.89) 33 (40.24)
ISUP 0.17 0.044
   1 8 (6.02) 1 (1.75) 8 (7.41) 1 (1.22)
   2 14 (10.53) 4 (7.02) 13 (12.04) 5 (6.09)
   3 13 (9.77) 10 (17.54) 15 (13.89) 8 (9.76)
   4 49 (36.84) 15 (26.32) 37 (34.26) 27 (32.92)
   5 49 (36.84) 27 (47.37) 35 (32.41) 41 (50.00)

Data are presented as median (IQR) or n (%). ALP, alkaline phosphatase; BM, bone metastasis; BMI, body mass index; DRE, digital rectal examination; FPSA, free prostate-specific antigen; IQR, interquartile range; ISUP, International Society of Urological Pathology; PCa, prostate cancer; PNC, positive core number; PSAD, prostate-specific antigen density; TPSA, total prostate-specific antigen.

Table 2

Univariate and multivariate logistic regression analyses for clinical and pathological characteristics

Variables Univariate logistic regression Multivariate logistic regression
B OR (95% CI) P B OR (95% CI) P
Age 0.02 1.02 (0.98–1.06) 0.38
ALP 0.01 1.01 (1.01–1.02) <0.01*
PNC 0.26 1.29 (1.17–1.44) <0.01* 0.14 1.15 (1.01–1.31) 0.04*
TPSA 0.03 1.03 (1.02–1.04) <0.01*
BMI 0.01 1.01 (0.94–1.08) 0.77
Volume 0 1.00 (0.99–1.01) 0.57
PSAD 0.76 2.14 (1.58–2.88) <0.01*
FPSA 0.11 1.12 (1.09–1.15) <0.01* 0.12 1.12 (1.05–1.20) <0.01*
DRE
   0 1.00 (reference)
   1 0.06 1.06 (0.59–1.90) 0.85
ISUP
   1 1.00 (reference)
   2 1.12 3.08 (0.30–31.33) 0.34
   3 1.45 4.27 (0.45–40.44) 0.21
   4 1.76 5.84 (0.69–49.48) 0.11
   5 2.24 9.37 (1.12–78.64) 0.04*

*, statistical significant value of P<0.05. ALP, alkaline phosphatase; BMI, body mass index; CI, confidence interval; DRE, digital rectal examination; FPSA, free prostate-specific antigen; ISUP, International Society of Urological Pathology; OR, odds ratio; PNC, positive core number; PSAD, prostate-specific antigen density; TPSA, total prostate-specific antigen.

Segmentation reproducibility

The manual segmentation demonstrated excellent reproducibility. The inter-reader analysis yielded a mean DSC of 0.78 (95% CI: 0.70–0.84) and an ICC of 0.77. The intra-reader analysis showed a mean DSC of 0.84 (95% CI: 0.79–0.90) and an ICC of 0.85.

Selection of habitat subregion and conventional radiomic features

All radiomic features used for model construction were extracted from T2WI images. Initially, 408 habitat subregion features and 136 whole-tumor radiomic features were extracted from each ROI. First, the SelectKBest method identified that 30 key features were highly correlated with the BM status. Subsequently, the mRMR method further reduced this set to 22 features having high correlation with the BM status and low redundancy. Finally, the LASSO method with an alpha value of 0.05 was used to select the optimal radiomic features. Using these 3 rounds of feature selection and dimensionality reduction, 6 habitat radiomic features and 5 whole-tumor radiomic features were ultimately chosen to construct the habitat model and the conventional radiomics model, respectively. Finally, the combined model was constructed using a late-fusion method integrating the habitat, whole-tumor, and clinical models.

Model development and validation

We evaluated ten machine learning algorithms—SVM, Autoencoder, Linear Discriminant Analysis (LDA), Random Forest, Logistic Regression, LASSO regression, Adaptive Boosting, Decision Tree, GP, and Naïve Bayes—using five-fold cross-validation on the training set. For each algorithm, hyperparameter optimization was performed via grid search with the following ranges: for GP, we tested RBF, Matern, and Rational Quadratic kernels with length scales ranging from 0.1 to 2.0; for SVM, we varied the regularization parameter C (0.1, 1, 10) and kernel type (linear, RBF); for Random Forest, we adjusted the number of trees (50, 100, 200) and maximum depth (5, 10, none); and for other algorithms, commonly used parameter grids were applied. The final model selection was based on the highest mean cross-validated AUC across the five folds. The GP classifier with an RBF kernel and a length scale of 0.8 achieved the highest mean cross-validated AUC (0.892) and was therefore selected as the optimal model. The optimized hyperparameters were then fixed, and the model was retrained on the full training set before evaluation on the internal test and external validation sets.

Figure 3A-3C displays the optimal ROC curves for the clinical, conventional radiomics, and habitat models, with detailed performance metrics provided in Table 3. The habitat and combined models consistently achieved higher AUCs across all cohorts. Specifically, in the training set, the habitat model attained an AUC of 0.870 (95% CI: 0.806–0.934) and the combined model reached 0.927 (0.883–0.971), outperforming the conventional radiomics (0.738, 0.654–0.822) and clinical models (0.787, 0.711–0.863). In the internal test set, the habitat and combined models achieved AUCs of 0.838 (0.733–0.942) and 0.841 (0.731–0.951), respectively, compared to 0.691 (0.553–0.829) for the conventional and 0.704 (0.566–0.842) for the clinical model. In the external validation set, the habitat model yielded an AUC of 0.760 (0.594–0.926), while the combined model achieved 0.821 (0.675–0.966), with corresponding accuracies of 0.787 and 0.851, sensitivities of 0.813 and 0.688, and specificities of 0.774 and 0.936.

Figure 3 ROC curves and DeLong test results. (A-C) ROC curves of the four models in the training, internal test, and external validation cohorts. (D-F) Bar plots showing pairwise DeLong test comparisons. Asterisks (*) indicate P<0.05. AUC, area under the receiver operating characteristic curve; ROC, receiver operating characteristic.

Table 3

Prediction performance of different models

Model Cohort AUC 95% CI Cutoff ACC Sen Spe PPV NPV
Clinical Train 0.787 0.711–0.863 0.576 0.737 0.649 0.803 0.712 0.753
Test 0.703 0.566–0.842 0.659 0.684 0.480 0.844 0.706 0.680
External 0.667 0.498–0.836 0.482 0.681 0.875 0.581 0.519 0.900
Radiomics Train 0.738 0.654–0.822 0.497 0.684 0.772 0.618 0.602 0.783
Test 0.691 0.553–0.829 0.531 0.667 0.680 0.656 0.607 0.724
External 0.688 0.515–0.860 0.500 0.702 0.687 0.709 0.550 0.814
Habitat Train 0.870 0.806–0.934 0.468 0.835 0.895 0.789 0.761 0.909
Test 0.838 0.733–0.942 0.502 0.789 0.840 0.750 0.724 0.857
External 0.760 0.594–0.926 0.513 0.787 0.813 0.774 0.650 0.889
Combined Train 0.927 0.883–0.971 0.502 0.879 0.842 0.908 0.873 0.885
Test 0.841 0.731–0.951 0.459 0.825 0.880 0.781 0.759 0.893
External 0.821 0.675–0.966 0.545 0.851 0.688 0.936 0.846 0.853

ACC, accuracy; AUC, area under the receiver operating characteristic curve; CI, confidence interval; NPV, negative predictive value; PPV, positive predictive value; Sen, sensitivity; Spe, specificity.

The DeLong test confirmed that in the training and test sets, the combined and habitat models significantly outperformed the conventional and clinical models (P<0.05, Figure 3D-3F). Although the habitat model’s superiority in the external validation set was not statistically significant, net reclassification improvement (NRI) analysis showed positive values versus the clinical (0.3333) and conventional (0.1795) models. The combined model exhibited the highest NRI values against all baseline models (Figure 4). DCA revealed greater net clinical benefit for the combined model (Figure 5A-5C), and calibration curves confirmed its superior precision (Figure 5D-5F). SHAP analysis identified original_firstorder_Skewness_h2 and original_gldm_DependenceVariance_h2 as the most influential features across all datasets, both exerting a negative effect on the model output (Figure 6).

Figure 4 Heatmap of NRI results between models. The color represents the magnitude of the NRI value: red indicates positive improvement (the model in the row outperforms that in the column), and blue indicates negative improvement. **, P<0.01; ***, P<0.001. NRI, net reclassification improvement.
Figure 5 Decision curve and calibration curve analysis. (A-C) DCA showing net clinical benefit. (D-F) Calibration curves showing agreement between predicted probabilities and observed outcomes. DCA, decision curve analysis.
Figure 6 SHAP analysis of the habitat model. (A) Summary plot showing the importance and directional effect of the top features. (B) Feature importance plot showing the relative contribution of the top features to the model’s predictions. SHAP, SHapley Additive exPlanations.

Habitat subregion feature analysis

This study conducted the ADC value distribution for the 3 clustered subregions (Habitat 1, Habitat 2, and Habitat 3) and plotted the corresponding histograms to further characterize the imaging properties of each habitat subregion (Figure S2). The analysis demonstrated that the three subregions exhibited significantly distinct ADC value (P<0.01) distribution patterns, supporting the validity of the clustering segmentation. Specifically, Habitat 2 demonstrated unique radiomics characteristics: its first-order energy was significantly higher than that of the other two subregions, whereas its global statistics (such as variance) were at a moderate level (Figure 7A-7C).

Figure 7 Distribution histograms of key radiomic features across habitat subregions. (A) First-order energy. (B) Skewness. (C) Variance. ADC, apparent diffusion coefficient.

Discussion

In this study, we developed and validated four predictive models for BM risk in PCa using bp-MRI. Both the habitat radiomics model and the combined model (integrating habitat, whole-tumor, and clinical features) significantly outperformed conventional radiomics and clinical models. Decision curve and calibration curve analyses confirmed the superior clinical utility and predictive accuracy of the combined model.

Across all datasets, the DeLong test showed that both the combined and habitat models significantly outperformed the clinical and conventional models, while no significant difference between the combined and habitat models themselves. This likely reflects the dominant predictive contribution of the habitat features, which may already capture the tumor heterogeneity most relevant to BM, limiting the incremental value of other features—consistent with reports that a small subset of discriminative features can suffice for optimal performance (27,28). Supporting this, a habitat-like analysis of prostate multiparametric magnetic resonance imaging (mp-MRI) identified subregions linked to high-grade disease (29), suggesting inherent suitability for capturing aggressive tumor phenotypes.

In the external validation set, the habitat model was not statistically superior to the conventional model by DeLong test, but the combined model showed significantly better discrimination. NRI analysis revealed that the habitat model improved risk stratification over the clinical (NRI =0.333) and conventional (NRI =0.179) models. The discrepancy between AUC comparison and positive NRI may be due to the limited sample size (n=48), reducing power to detect modest effect, as well as potential cross-center variability affecting habitat features stability (30,31). Nonetheless, the habitat model maintained clinically relevant discrimination (AUC =0.760) in this independent cohort, and its integration into the combined model yielded the most robust performance.

Our “functional-structural” strategy leverages complementary bp-MRI information. ADC-based subregion delineation provides a functionally informed tumor partition reflecting cellular density and microstructural heterogeneity (14), while texture analysis of T2WI within habitats characterizes subtle structural organization (17). This approach aligns with habitat imaging rationale linking distinct phenotypes to tumor biology and offers a clinically intuitive tool for visualizing functionally distinct subregions (19,22,32).

SHAP analysis identified original_firstorder_Skewness_h2 and original_gldm_DependenceVariance_h2 as the most influential features, with negative SHAP values indicating an inverse association with BM risk. Habitat 2 showed a “high energy-low heterogeneity” pattern: its first-order energy was higher than other subregions while variance was moderate. Higher energy reflects concentrated intensity distribution, suggesting a uniform, metabolically active cell population (17). Lower skewness and lower dependence variance in Habitat 2—key predictors of BM risk—corresponding to symmetrical intensity distribution and reduced local texture heterogeneity. This pattern aligns with homogeneous, densely cellular tumor tissue, which may reflect the presence of a dominant tumor clone (33,34) associated with increased invasive and metastatic potential (32,35). Conversely, higher skewness and variance could be related to greater tumor microenvironment heterogeneity (e.g., immune or stromal components), which might restrain progression, but this interpretation remains speculative (36,37). These findings suggest that within specific habitats, spatial homogeneity rather than heterogeneity could be more closely associated with aggressive tumor behavior. This observation offers a hypothesis-generating perspective on the clinical implications of tumor heterogeneity, warranting further investigation (32).

The model’s generalizability was supported by the external validation cohort (n=48), where the combined model achieved an AUC of 0.821 and the habitat model 0.760. SHAP analysis confirmed the consistent influence of Habitat 2 features across centers, supporting cross‑center reliability.

This study has several limitations. First, its retrospective design and limited primary cohort (n=190) may introduce selection bias. The small external validation set (n=48) reduced statistical power, potentially explaining the non-significant difference between the habitat and conventional models externally (30). Moreover, despite rigorous feature selection, the sample size relative to the number of extracted radiomic features may carry a residual overfitting risk. Future large-scale, multicenter prospective studies are needed for further validation. Second, despite standardized preprocessing, variability in image acquisition parameters across centers remains a challenge despite preprocessing. Third, tumor ROI delineation, though performed by experienced radiologists in consensus, remains subjective, and inter-observer variability was not formally assessed. Additionally, radiomic feature extraction and model building were performed after outcome labels were known; formal blinding of predictors was not implemented. Fourth, BM diagnosis relied primarily on bone scintigraphy, which has false-positive and sensitivity limitations. While MRI confirmed negative cases, positive cases lacked prostate-specific membrane antigen (PSMA) PET/CT or histopathologic confirmation—a major limitation we acknowledge. Future studies should use a composite reference standard (PSMA PET/CT, follow-up imaging, or histopathology). Finally, biological interpretations are preliminary and hypothesis-generating; proposed links between imaging features and tumor subclones or microenvironment require direct histopathologic or radiogenomic validation. Additionally, several limitations concerning study methodology should be acknowledged. No a priori sample size calculation was performed; the sample size was determined by consecutively available eligible patients. Missing data were handled by median imputation, but the number and pattern of missing values were not reported. Outcome assessment was not blinded to predictor information. The full prediction model (coefficients or formula) is not provided, which limits independent reproducibility, and instructions for individual prediction are not available.


Conclusions

In conclusion, our bp-MRI-based habitat radiomics model shows promise for noninvasive BM risk assessment in PCa. The combined model, integrating habitat, whole-tumor, and clinical features, provides favorable predictive performance and may assist in clinical decision-making. Prospective multicenter validation is warranted before clinical adoption.


Acknowledgments

The authors gratefully acknowledge the technical support provided by the Department of Medical Imaging at Gansu Provincial Hospital. We also thank all the patients who participated in this study.


Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tau.amegroups.com/article/view/10.21037/tau-2026-0380/rc

Data Sharing Statement: Available at https://tau.amegroups.com/article/view/10.21037/tau-2026-0380/dss

Peer Review File: Available at https://tau.amegroups.com/article/view/10.21037/tau-2026-0380/prf

Funding: This work was supported by the Gansu Provincial Major Project for Sci-Tech Innovation in the Health Industry (grant No. GSWSZD2024-03), Zhangye City Science and Technology Plan Project (grant No. ZY2025BJ45), and the Gansu University of Chinese Medicine Graduate Innovation and Entrepreneurship Project (grant No. 2026CXCY-062).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tau.amegroups.com/article/view/10.21037/tau-2026-0380/coif). K.A. is an employee of Philips Healthcare (Suzhou) Co., Ltd. The other authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. This study was approved by the Ethics Committee of Gansu Provincial Hospital (approval No. 2025-489) and the Ethics Committee of Zhangye People’s Hospital Affiliated to Hexi University (approval No. HFYER-2023.07). Informed consent was waived in this retrospective study.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.


References

  1. Bray F, Laversanne M, Sung H, et al. Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2024;74:229-63. [Crossref] [PubMed]
  2. Temporal patterns of cancer burden in Asia, 1990-2019: a systematic examination for the Global Burden of Disease 2019 study. Lancet Reg Health Southeast Asia 2024;21:100333.
  3. Halabi S, Kelly WK, Ma H, et al. Meta-Analysis Evaluating the Impact of Site of Metastasis on Overall Survival in Men With Castration-Resistant Prostate Cancer. J Clin Oncol 2016;34:1652-9. [Crossref] [PubMed]
  4. Sartor O, de Bono JS. Metastatic Prostate Cancer. N Engl J Med 2018;378:645-57. [Crossref] [PubMed]
  5. Tilki D, van den Bergh RCN, Briers E, et al. EAU-EANM-ESTRO-ESUR-ISUP-SIOG Guidelines on Prostate Cancer. Part II-2024 Update: Treatment of Relapsing and Metastatic Prostate Cancer. Eur Urol 2024;86:164-82.
  6. Hupe MC, Philippi C, Roth D, et al. Expression of Prostate-Specific Membrane Antigen (PSMA) on Biopsies Is an Independent Risk Stratifier of Prostate Cancer Patients at Time of Initial Diagnosis. Front Oncol 2018;8:623. [Crossref] [PubMed]
  7. Lin Y, Mao Q, Chen B, et al. When to perform bone scintigraphy in patients with newly diagnosed prostate cancer? a retrospective study. BMC Urol 2017;17:41.
  8. Wang S, Burtt K, Turkbey B, et al. Computer aided-diagnosis of prostate cancer on multiparametric MRI: a technical review of current research. Biomed Res Int 2014;2014:789561. [Crossref] [PubMed]
  9. Ruprecht O, Weisser P, Bodelle B, et al. MRI of the prostate: interobserver agreement compared with histopathologic outcome after radical prostatectomy. Eur J Radiol 2012;81:456-60. [Crossref] [PubMed]
  10. Padhani AR, Lecouvet FE, Tunariu N, et al. METastasis Reporting and Data System for Prostate Cancer: Practical Guidelines for Acquisition, Interpretation, and Reporting of Whole-body Magnetic Resonance Imaging-based Evaluations of Multiorgan Involvement in Advanced Prostate Cancer. Eur Urol 2017;71:81-92. [Crossref] [PubMed]
  11. Harvey H, deSouza NM. The role of imaging in the diagnosis of primary prostate cancer. J Clin Urol 2016;9:11-7. [Crossref] [PubMed]
  12. Turkbey B, Rosenkrantz AB, Haider MA, et al. Prostate Imaging Reporting and Data System Version 2.1: 2019 Update of Prostate Imaging Reporting and Data System Version 2. Eur Urol 2019;76:340-51. [Crossref] [PubMed]
  13. van der Leest M, Cornel E, Israël B, et al. Head-to-head Comparison of Transrectal Ultrasound-guided Prostate Biopsy Versus Multiparametric Prostate Resonance Imaging with Subsequent Magnetic Resonance-guided Biopsy in Biopsy-naïve Men with Elevated Prostate-specific Antigen: A Large Prospective Multicenter Clinical Study. Eur Urol 2019;75:570-8. [Crossref] [PubMed]
  14. Woo S, Suh CH, Kim SY, et al. Head-to-Head Comparison Between Biparametric and Multiparametric MRI for the Diagnosis of Prostate Cancer: A Systematic Review and Meta-Analysis. AJR Am J Roentgenol 2018;211:W226. [Crossref] [PubMed]
  15. Weinreb JC, Barentsz JO, Choyke PL, et al. PI-RADS Prostate Imaging - Reporting and Data System: 2015, Version 2. Eur Urol 2016;69:16-40. [Crossref] [PubMed]
  16. Aerts HJ, Velazquez ER, Leijenaar RT, et al. Decoding tumour phenotype by noninvasive imaging using a quantitative radiomics approach. Nat Commun 2014;5:4006. [Crossref] [PubMed]
  17. van Griethuysen JJM, Fedorov A, Parmar C, et al. Computational Radiomics System to Decode the Radiographic Phenotype. Cancer Res 2017;77:e104-7. [Crossref] [PubMed]
  18. Wang Y, Yu B, Zhong F, et al. MRI-based texture analysis of the primary tumor for pre-treatment prediction of bone metastases in prostate cancer. Magn Reson Imaging 2019;60:76-84. [Crossref] [PubMed]
  19. Lambin P, Leijenaar RTH, Deist TM, et al. Radiomics: the bridge between medical imaging and personalized medicine. Nat Rev Clin Oncol 2017;14:749-62. [Crossref] [PubMed]
  20. Bedard PL, Hansen AR, Ratain MJ, et al. Tumour heterogeneity in the clinic. Nature 2013;501:355-64. [Crossref] [PubMed]
  21. Lin S, He P, You R. Prediction of bone metastasis of prostate cancer based on intratumoral and peritumoral radiomics of MRI T2WI combined with ADC images. Front Oncol 2025;15:1555315. [Crossref] [PubMed]
  22. Braman NM, Etesami M, Prasanna P, et al. Intratumoral and peritumoral radiomics for the pretreatment prediction of pathological complete response to neoadjuvant chemotherapy based on breast DCE-MRI. Breast Cancer Res 2017;19:57. [Crossref] [PubMed]
  23. Zhou M, Hall L, Goldgof D, et al. Radiologically defined ecological dynamics and clinical outcomes in glioblastoma multiforme: preliminary results. Transl Oncol 2014;7:5-13. [Crossref] [PubMed]
  24. Kickingereder P, Burth S, Wick A, et al. Radiomic Profiling of Glioblastoma: Identifying an Imaging Predictor of Patient Survival with Improved Performance over Established Clinical and Radiologic Risk Models. Radiology 2016;280:880-9. [Crossref] [PubMed]
  25. Beig N, Khorrami M, Alilou M, et al. Perinodular and Intranodular Radiomic Features on Lung CT Images Distinguish Adenocarcinomas from Granulomas. Radiology 2019;290:783-92. [Crossref] [PubMed]
  26. Rosenkrantz AB, Neil J, Kong X, et al. Prostate cancer: Comparison of 3D T2-weighted with conventional 2D T2-weighted imaging for image quality and tumor detection. AJR Am J Roentgenol 2010;194:446-52. [Crossref] [PubMed]
  27. Liu J, Guo W, Zeng P, et al. Vertebral MRI-based radiomics model to differentiate multiple myeloma from metastases: influence of features number on logistic regression model performance. Eur Radiol 2022;32:572-81. [Crossref] [PubMed]
  28. Zhang W, Liang F, Zhao Y, et al. Multiparametric MR-based feature fusion radiomics combined with ADC maps-based tumor proliferative burden in distinguishing TNBC versus non-TNBC. Phys Med Biol 2024; [Crossref]
  29. Fehr D, Veeraraghavan H, Wibmer A, et al. Automatic classification of prostate cancer Gleason scores from multiparametric magnetic resonance images. Proc Natl Acad Sci U S A 2015;112:E6265-E6273. [Crossref] [PubMed]
  30. van Smeden M, Moons KG, de Groot JA, et al. Sample size for binary logistic prediction models: Beyond events per variable criteria. Stat Methods Med Res 2019;28:2455-74. [Crossref] [PubMed]
  31. Traverso A, Wee L, Dekker A, et al. Repeatability and Reproducibility of Radiomic Features: A Systematic Review. Int J Radiat Oncol Biol Phys 2018;102:1143-58. [Crossref] [PubMed]
  32. Junttila MR, de Sauvage FJ. Influence of tumour micro-environment heterogeneity on therapeutic response. Nature 2013;501:346-54. [Crossref] [PubMed]
  33. Nowell PC. The clonal evolution of tumor cell populations. Science 1976;194:23-8. [Crossref] [PubMed]
  34. Gerlinger M, Rowan AJ, Horswell S, et al. Intratumor heterogeneity and branched evolution revealed by multiregion sequencing. N Engl J Med 2012;366:883-92. [Crossref] [PubMed]
  35. Migliozzi S, Adabbo B, Garofano L, et al. Restraint of cancer cell plasticity by spatial homotypic clustering. Cancer Cell 2025;43:2206-2223.e10. [Crossref] [PubMed]
  36. Ganeshan B, Panayiotou E, Burnand K, et al. Tumour heterogeneity in non-small cell lung carcinoma assessed by CT texture analysis: a potential marker of survival. Eur Radiol 2012;22:796-802. [Crossref] [PubMed]
  37. Galon J, Costes A, Sanchez-Cabo F, et al. Type, density, and location of immune cells within human colorectal tumors predict clinical outcome. Science 2006;313:1960-4. [Crossref] [PubMed]
Cite this article as: Gong J, Li F, Sun Y, Huang J, Zhang Y, Yang O, Li J, Song Z, Ai K, Huang G. An interpretable biparametric MRI habitat radiomics model for predicting bone metastasis in prostate cancer: a dual-center study. Transl Androl Urol 2026;15(8):273. doi: 10.21037/tau-2026-0380

Download Citation