Prognostic prediction in nasopharyngeal carcinoma using multi-region interaction features based on MRI habitat analysis

Introduction

Nasopharyngeal carcinoma (NPC) is a common malignant tumor of the head and neck. According to data from the International Agency for Research on Cancer, NPC demonstrates significant regional disparities, with a high incidence in Southeast Asia and East Asia (1). Due to the insidious onset of NPC, nonspecific clinical symptoms, and the low positive predictive value (PPV) of routine screening methods, over 70% of NPC patients are diagnosed at locally advanced stages with poor therapeutic efficacy, significantly compromising treatment outcomes (2,3). Furthermore, modern tumor biology studies (4) have confirmed that significant spatial heterogeneity exists within tumors, and cells in distinct intratumoral regions may harbor different genetic mutations. The distribution of these mutations exhibits marked spatial variation, which influences tumor proliferation rate, invasiveness, and sensitivity to drug therapy and radiotherapy. Therefore, intratumoral spatial heterogeneity not only affects treatment selection and efficacy evaluation but also directly determines patient prognosis and survival. Consequently, more precise and noninvasive identification of spatial heterogeneity within NPC tumors is critical for prognostic prediction in NPC patients.

In clinical practice, magnetic resonance imaging (MRI), computed tomography (CT), and positron emission tomography (PET) are commonly used radiological imaging modalities for NPC evaluation, providing morphological and functional information of lesion regions (5,6). Radiomics, through its standardized workflow of image acquisition, feature extraction, feature selection, and model construction (7), has demonstrated significant potential in predicting treatment outcomes and cancer genetics. Recent studies indicate that MRI radiomics has been widely applied in the early diagnosis, clinical staging, treatment efficacy evaluation, prediction of post-radiotherapy adverse reactions, and prognosis of NPC (8-13), and has achieved good results. Among these, Jiang et al. (13) developed and validated a radiomics model based on pre-treatment MRI, which performed well in predicting the overall survival (OS) of NPC patients with local residual lesions after radiotherapy, with a C-index of 0.788. Raghavan Nair et al. (14) extracted radiomics features from T1WI, T2WI, and Contrast-enhanced T1-weighted imaging (CE-T1WI) images of 58 NPC patients, confirming that radiomics features have high predictive value for local recurrence of NPC. Cao et al. (15) predicted distant metastasis in early-stage NPC using multiparametric MRI radiomics combined with clinical features. The results showed that the radiomics combined model could effectively distinguish between patients with distant metastasis and non-metastatic patients. These studies all demonstrate that radiomics models perform well in prognostic prediction for NPC patients.

However, traditional radiomics methods often extract a large number of high-throughput quantitative features from conventional images based on the entire tumor region, including information on tumor size, shape, intensity, and texture, but they cannot reflect the information of fine subregions within the tumor (16). Current studies have considered the significant intratumoral heterogeneity in NPC (17), but traditional radiomics commonly summarizes tumor heterogeneity using global whole-tumor features, but lacks explicit spatial sub-regional or localized heterogeneity modeling, limiting its ability to characterize fine-grained intra-tumoral spatial complexity. Therefore, we need to explore new non-invasive methods to quantify intratumoral heterogeneity in NPC.

In recent years, habitat imaging uses imaging biomarkers and multiparametric voxel-level clustering to delineate distinct intratumoral sub-regions with similar imaging and putative biological characteristics (16). Multiple studies have demonstrated that in various malignancies including gliomas, brain metastases, breast cancer, cervical cancer, and hepatocellular carcinoma, MRI-based habitat imaging identifies distinct habitat regions through tumor clustering, and the established prognostic prediction models effectively predict tumor outcomes (18-25). Currently, studies have explored habitat analysis for predicting treatment response to neoadjuvant chemotherapy in NPC (26). However, research on habitat analysis-based models for prognostic prediction in NPC patients remains limited. Additionally, Cho et al. (22) identified five distinct habitats from early and delayed phases of DCE-MRI, constructing a habitat model to predict disease-free survival (DFS). The results showed that the habitat model served as an independent risk factor for poor DFS outcomes and outperformed radiomics models. Wu et al. (27) partitioned each breast cancer tumor into multiple spatially separated subregions based on enhanced DCE-MRI images, defined a multi-regional spatial interaction (MSI) matrix, extracted 22 imaging features, and stratified patients into different risk groups using network analysis. The prognostic value of imaging stratification alongside clinicopathological factors was assessed using multivariate Cox regression.

In this study, we aim to utilize MRI-based habitat analysis to partition NPC tumors into multiple subregions and to quantify the spatial interactions among these subregions by extracting MSI features that characterize intratumoral heterogeneity. The primary objective is to develop and validate an MSI-based predictive model for mortality risk prediction and subsequent risk stratification in patients with NPC. As a secondary objective, survival analysis using Cox proportional hazards models is performed to further evaluate the prognostic value and clinical utility of MSI features. By providing a more explicit and interpretable representation of intratumoral spatial heterogeneity, this framework aims to support more informed individualized treatment decision-making. We present this article in accordance with the TRIPOD reporting checklist (available at https://qims.amegroups.com/article/view/10.21037/qims-2025-1-2551/rc).

Methods

Study design

This study was a retrospective, imaging-based prognostic study in patients with NPC. Eligible patients with pre-treatment MRI data and complete follow-up information were included in the analysis. The primary endpoint was defined as all-cause mortality during the entire available follow-up period. Patients were classified as ‘death’ if death occurred at any time during follow-up, and as ‘survival’ if they were alive at the last follow-up, in which case the outcome was treated as censored. MSI features were extracted from pre-treatment T2-weighted MRI (T2W-MRI) images (detailed acquisition parameters are provided in Appendix 1). Based on these features, feature selection and machine learning model construction were performed to generate individualized predicted probabilities for mortality risk, which were used to estimate each patient’s risk of death. The binary classification task was designed to model overall mortality risk rather than survival at a predefined time horizon. In addition, Cox proportional hazards regression was performed as a complementary time-to-event analysis to evaluate the association between model-derived risk scores and time-to-event outcomes, thereby assessing the prognostic relevance of the proposed model. An overview of the study workflow is shown in Figure 1.

Figure 1 Experimental workflow. 2D, two-dimensional; 3D, three-dimensional; MRI, magnetic resonance imaging; MSI, multi-regional spatial interaction; ROC, receiver operating characteristic.

Study patients

This retrospective study enrolled 1,297 patients with NPC who were consecutively admitted to Sun Yat-sen University Cancer Center (SYSUCC) between January 2010 and March 2014 (Table 1). The study was approved by the Ethics Committee of Sun Yat-sen University Cancer Center (No. B2019-222-01) and was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. Given the retrospective nature of the study, the requirement for written informed consent was waived by the ethics committee. The inclusion and exclusion criteria were consistent with those defined in a previous SYSUCC study on NPC prognosis (28).

Treatment

All patients underwent intensity-modulated radiotherapy (IMRT). Treatment strategies, including induction chemotherapy or concurrent chemoradiotherapy (CCRT), were based on the American Joint Committee on Cancer (AJCC) 8th edition staging system and institutional guidelines. Patients with stage II–IV disease received CCRT, with induction chemotherapy. Specifically, 11.7% of patients were treated with RT alone, 38.3% with CCRT alone, and 50.0% with induction chemotherapy plus CCRT. The supplementary material provides detailed technical parameters for IMRT, CCRT, and induction chemotherapy.

Data splitting

Among the cohort, 1,140 cases were non-deceased patients, and 157 cases were deceased patients. To address significant class imbalance in the dataset, we proposed a stratified sampling-based data partitioning scheme (Figure 2A). Stratified sampling (29), a widely used statistical method, aims to enhance sample representativeness and estimation accuracy. Its core principle involves dividing the population into homogeneous subgroups (termed “strata”) and independently sampling from each stratum. First, 50 deceased and 50 non-deceased patients were randomly selected to form an independent test set. The remaining 1,090 non-deceased patients were then randomly divided into 10 equal subsets, each of which was combined with the same set of 107 remaining deceased patients to construct 10 balanced training datasets (216 cases per set). Although deceased samples were repeatedly used across training datasets, this design was intentionally chosen to reduce training instability under extreme class imbalance and to preserve statistical power given the limited number of outcome events. Each training dataset was processed independently for feature selection and model training, and final predictions were obtained through an ensemble of models. All performance evaluations were conducted exclusively on the independent test set to ensure an unbiased assessment of model generalization performance.

Figure 2 Dataset partitioning and model ensemble workflow. (A) The imbalanced dataset containing death prognosis information is partitioned into 10 training subsets and 1 testing subset. (B) Individual models are trained on each training subset, followed by performing model ensemble through weighted averaging based on AUC scores to generate the final predictive model.

Partitioning of intra-tumoral subregions

After acquiring the original T2W-MRI images and physician-delineated three-dimensional (3D) NPC tumor regions, we developed a tumor intraregional subpartitioning method. Detailed T2W-MRI acquisition parameters are provided in the supplementary materials. As shown in Figure S1A, voxel-wise intensity values extracted from the 3D region of interest (ROI) across all patients in the training cohort were normalized using z-score normalization, in order to eliminate inter-patient differences in gray-level distribution. Subsequently, for each voxel within the tumor volume, the gray intensity of the central voxel within each tumor region and the gray-level differences between this voxel and its four in-plane neighboring voxels were extracted to form a global feature matrix (Figure S1B), which was used as the input for K-Means clustering (30). Although tumor delineation was performed on the 3D volume, in-plane neighboring voxels were considered to capture local spatial intensity variations while maintaining robustness to inter-slice spacing differences. The global label matrix obtained from clustering enabled the target region tumor volume to be divided into subregions with relatively homogeneous imaging characteristics and local spatial consistency, thereby achieving intratumoral subregion partitioning for individual tumors based on population voxel-level characteristics. Furthermore, the optimal number of clusters was evaluated within the range of 2–10, and the corresponding Calinski-Harabasz (CH) index (31) was calculated using the scikit-learn package in Python to assess clustering quality by measuring intra-cluster compactness and inter-cluster separation. A higher CH index indicates better clustering performance; therefore, the optimal cluster number k was determined based on the training cohort data. After determining the optimal k, both the training and test cohorts employed an identical clustering procedure, and the same optimal k was uniformly applied for K-Means clustering to generate tumor habitat maps for each cohort, thereby ensuring consistency and comparability in subregion partitioning between the two datasets.

Quantitative MSI features

Based on the generated habitat maps, a MSI matrix (27) was independently constructed for each tumor to quantify intratumoral spatial heterogeneity. Specifically, each tumor voxel was assigned a habitat label after clustering. Subsequently, for each image, a local neighborhood search was performed for every tumor voxel to identify the “cluster label pairs” between that voxel and its neighboring voxels. Each time a label pair (i,j) occurred, the corresponding entry (i,j) in the MSI matrix was incremented by one. This process was repeated for all voxels until the entire ROI was traversed (as shown in Figure 1). In the MSI matrix, the diagonal elements represent the spatial continuity within the same subregion, i.e., the volume of each habitat, whereas the off-diagonal elements reflect the interaction strength between different subregions, indicating the extent of spatial boundaries where different habitats are adjacent. Finally, 45 MSI features were extracted from each tumor’s MSI matrix, including 38 first-order statistical features and 7 second-order texture features (Table S1), to quantify intratumoral spatial heterogeneity from multiple perspectives. As comparative features, radiomics features were simultaneously extracted from each patient’s tumor ROI using the open-source PyRadiomics package, including shape, gray-level histogram, and multiple texture categories (gray-level co-occurrence matrix, gray-level run length matrix, gray-level size zone matrix, neighboring gray-tone difference matrix, gray-level dependence matrix). A total of 1,710-dimensional radiomics features were obtained and used for performance comparison with the model constructed using MSI features.

Feature selection

To ensure a fair comparison between MSI features and radiomics features, an identical feature selection workflow was applied to both feature types and independently executed across ten balanced training subsets. First, univariate statistical testing (independent-samples t-test) was performed to remove features without significant differences between the survival and death groups. Subsequently, highly collinear features were eliminated using Pearson correlation, with a correlation threshold of |r|>0.8, to reduce redundancy and improve model stability. The chi-squared test was then applied to the remaining features to further select those significantly associated with mortality.

Feature selection was performed independently within each training subset. Features that were consistently retained across the training subsets were considered stable and were used as inputs for model construction. Radiomics-only and MSI-only models were built using their respective stable feature sets, while the combined model was constructed by concatenating the selected radiomics and MSI features.

Model construction and ensemble integration

After feature selection, all retained features were normalized to the range of 0–1, with normalization parameters computed from the training set and subsequently applied to the test set. The dataset was divided into ten sub-training sets and one independent test set according to the partitioning strategy described in this study, and a base classifier was trained separately on each sub-training set.

To comprehensively evaluate the predictive value and robustness of the proposed MSI features, multiple commonly used machine learning algorithms representing different modeling paradigms—including decision tree, logistic regression, support vector machine, random forest, gradient boosting, AdaBoost, Naïve Bayes, XGBoost, LightGBM, and CatBoost—were implemented and compared. All models were trained and evaluated under identical experimental settings, including the same data splits, feature selection workflow, and evaluation metrics, to ensure a fair comparison focused on feature effectiveness rather than classifier-specific advantages.

Decision tree classifiers are capable of handling heterogeneous data and modeling nonlinear relationships ten balanced sub-training sets and have been applied in various medical prognostic prediction studies (32), particularly in settings where model interpretability is required. In this study, a decision tree classifier was employed to construct a mortality prediction model based on MSI features. For comparison, additional machine learning classifiers were trained using the same feature sets and experimental protocol, and their predictive performances were systematically evaluated. Radiomics-based prediction models were constructed following the same modeling and evaluation framework, and the classifier with the best performance was selected for subsequent ensemble modeling. For each type of feature, an individual machine learning model (referred to as a base model) was independently trained on each of the ten sub-training sets using the same classifier, identical feature sets, and the same training protocol. To enhance the robustness of the ensemble model, ten decision tree base learners were trained on ten balanced sub-training sets. Five-fold cross-validation was performed on each sub-training set, and the average AUC was calculated as the performance metric. The predictive outputs of the ten decision tree base models were then integrated using a weighted averaging strategy based on the five-fold average AUC of the training set, where five-fold cross-validation was adopted to provide a stable performance estimate under limited and imbalanced outcome events (33). The weight w_i was computed as follows:

wi=max(0,AUCi−0.5)∑j=110max(0,AUCj−0.5)[1]

where AUC_i represents the mean AUC of the i-th base model obtained from five-fold cross-validation, and the denominator sums the truncated AUC contributions of all base models (j=1,2,...,10). The max function is used to truncate models performing below the random level (0.5). The final ensemble predicted probability was obtained by weighted summation (Figure 2B):

P=∑i=110wipi[2]

The primary task of this study was binary classification of mortality outcome. Patients with a predicted probability of death greater than 0.5 were classified as high-risk (predicted death), while those with a predicted probability ≤0.5 were classified as low-risk (predicted survival).

In addition to the classification task, Cox proportional hazards regression was performed as a supplementary survival analysis to further evaluate the prognostic relevance of the imaging-derived features. For Cox modeling, the selected features, together with survival time and survival status, were used to train a Cox proportional hazards model within each training subset. Five-fold cross-validation was conducted for each subset, and the concordance index (C-index) was calculated to assess survival prediction performance.

The C-index obtained from cross-validation was used to determine the weight of each Cox model during ensemble integration, ensuring that models with superior prognostic performance contributed more to the final survival model. The ensemble Cox model was constructed by weighted aggregation of the sub-models based on their respective C-index values.

In the test set, each patient received a final survival risk score derived from the ensemble Cox model. Patients were stratified into high-risk and low-risk groups using a predefined cutoff determined by the median risk score in the training data. This fixed cutoff was applied to the test set for Kaplan-Meier survival analysis and log-rank testing to assess the prognostic stratification capability of the model.

Model performance evaluation and statistical analysis

The performance of classifiers in predicting death among NPC patients was evaluated using the AUC. Confusion matrices were plotted to calculate model accuracy, recall (sensitivity), specificity, precision (PPV), and negative predictive value (NPV). Calibration curves were generated in the test set to assess the agreement between predicted probabilities and observed event frequencies. Decision curve analysis (DCA) quantified the net benefit across threshold probabilities to evaluate the clinical utility of the ensemble models. Log-rank tests were applied to Kaplan-Meier survival curves to compare survival differences between high- and low-risk groups. The discriminative ability of survival models was measured using the C-index, which reflects the alignment between predicted risk scores and actual survival outcomes. Categorical variables were analyzed with chi-square or Fisher’s exact tests, while continuous variables were assessed via the Mann-Whitney U test. All statistical tests were two-sided, with P<0.05 indicating statistical significance.

In addition to overall performance evaluation, a predefined subgroup analysis was conducted to assess the robustness of the mortality prediction model across different tumor stages. Specifically, patients in the independent test set were stratified according to T stage (T1–T4), and the predictive performance of the final trained model was evaluated separately within each subgroup. This analysis was prespecified as a secondary, exploratory evaluation and was not used for model training, parameter tuning, or feature selection.

Results

Patient demographics

A total of 1,297 patients with NPC were included in the analysis. The median age was 45 years, and 73.4% of the patients were male. According to the AJCC 8th edition staging system, the majority of patients had stage II–IVa disease, with 286 patients (22.0%) classified as stage II and 900 patients (69.4%) as stage III–IVa. Treatment distributions and other baseline characteristics are summarized in Table 1.

Determining the number of clusters

By applying the K-Means clustering algorithm to the global feature matrix of tumor voxels from patients in the training set, the CH index was calculated for cluster numbers ranging from 2 to 10. As shown in Figure S2, the CH index reached its maximum when the number of clusters was 4, indicating that intra-cluster compactness and inter-cluster separation were optimal at this cluster number. Therefore, four clusters were selected as the number of habitats for segmentation. After performing the overall clustering, each tumor voxel was assigned to a habitat region, labeled as 1, 2, 3, or 4. Mapping the cluster labels of tumor voxels back to the original images for each patient generated color-coded habitat maps (Figure S1): red represents habitat 1, green represents habitat 2, blue represents habitat 3, and yellow represents habitat 4. These results demonstrate that the K-Means clustering method can segment tumor subregions at the population level with spatially uniform distribution and consistent signal intensity.

Selected optimal features

Selected MSI features

Based on the habitat maps generated through clustering, we constructed MSI matrices and extracted 45 MSI features to quantify intratumoral spatial heterogeneity. As shown in Figure S3, each patient’s MSI matrix displays 10 first-order features, where the diagonal elements represent the internal volume of each individual habitat region, and the off-diagonal elements represent the volumes of the boundaries between corresponding pairs of habitat regions.

During the feature selection stage, the 10 sub-training sets independently underwent feature screening, with each sub-training set evaluating the association between features and mortality outcomes using the Chi-squared test. After comprehensive analysis, a final set of 8 selected features was determined for subsequent model construction, including Boundary Volume Class 3–4, Contrast, Correlation, Diag Mean, Energy, Entropy, Kurtosis, and Volume Class 4. The AUC-weighted ensemble strategy balances the contributions of each sub-model to the final prediction, thereby improving prediction stability and robustness. Specifically, Boundary Volume Class 3–4 represents the boundary volume between habitat regions 3 and 4, Diag Mean represents the average of the total volumes across all habitat regions, and Volume Class 4 represents the overall volume of habitat region 4. As shown in Figure 3, the original tumor images and clustered habitat maps of two NPC patients illustrate differences in tumor heterogeneity: Patient 1 (deceased) exhibited higher Boundary Volume Class 3–4 and Volume Class 4 values compared to Patient 2 (survived). Further analysis in the training set (Figure S3) confirmed that these two features significantly differed between high- and low-mortality risk groups (P<0.05), suggesting that larger Volume Class 4 and Boundary Volume Class 3–4 may indicate higher mortality risk and possess potential prognostic value.

Figure 3 Tumor imaging and clustered habitat maps for two patients with NPC. (A) Schematic illustration of the whole tumor region for Patient 1. (B) Habitat clustering map of the entire tumor for Patient 1. (C) Visualization of the feature Boundary Volume Class 3–4 for Patient 1. (D) Visualization of the feature Volume Class 4 for Patient 1. (E) Schematic illustration of the whole tumor region for Patient 2. (F) Habitat clustering map of the entire tumor for Patient 2. (G) Visualization of the feature Boundary Volume Class 3–4 for Patient 2. (H) Visualization of the feature Volume Class 4 for Patient 2. The lower-right corner of each panel shows a 3D rendering, while the feature values derived in this study are indicated in parentheses in the lower-left corner. Each row represents different imaging views from the same patient, and the columns correspond to the original tumor image, the 3D habitat map (all regions), the combined visualization of Habitats 3 and 4, and Habitat 4 in isolation. Volume values shown in the figure represent the number of voxels (voxel count). NPC, nasopharyngeal carcinoma.

Selected radiomic features

Meanwhile, during the radiomics feature selection process, we employed a consistent feature selection workflow, integrating the results from the 10 sub-training sets through a union operation. This yielded a final set of 10 radiomics features, including the number of voxels in the original mask, first-order minimum value after exponential transformation, non-uniformity and dependence emphasis features from the gray-level dependence matrix in the gradient/wavelet domains (reflecting large dependence–high gray-level structures and small dependence–low gray-level patterns across multiple sub-bands), the total energy of the two-dimensional local binary pattern, and intensity and low-gray-region distribution features from the neighborhood gray-tone difference matrix in the wavelet domain.

Performance evaluation of MSI model

Across the 10 sub-training sets, all three feature-based models (MSI, radiomics, and combined features) demonstrated stable predictive performance, indicating robustness to sample perturbations. Among them, the MSI feature–based model consistently achieved the highest AUC values, with ROC curves showing a more concentrated distribution across training subsets (Figure 4), suggesting stronger and more stable discriminative ability for mortality risk compared with radiomics and combined features.

Figure 4 Training ROC curves of the three feature sets across 10 sub-training subsets: (A) MSI, (B) radiomics, and (C) combined features. Each curve represents a model trained on one subset, and the corresponding legend reports the AUC along with its 95% CI. AUC, area under the curve; CI, confidence interval; MSI, multi-regional spatial interaction; ROC, receiver operating characteristic.

In the independent test set, the MSI-based model achieved an AUC of 0.79, with an accuracy of 0.75, recall of 0.80, and negative predictive value of 0.78 (Table 2). When compared with models constructed using radiomics features alone and combined features, the MSI-based habitat model exhibited superior overall predictive performance (Figure 5A). The best-performing radiomics-based model achieved an AUC of 0.77 (Table S2), while the combined-feature model achieved an AUC of 0.77 (Table S3). The addition of radiomics features to MSI features did not further improve predictive performance, suggesting that MSI features capture complementary spatial heterogeneity information that may be diluted by redundant or noisy features when simple feature concatenation is applied.

Figure 5 Comparative visualization of MSI model, radiomics model and combined feature model. (A) ROC curve comparison; (B) calibration curve comparison; (C) decision curve comparison; (D) survival curves of the MSI model. AUC, area under the curve; CI, confidence interval; HR, hazard ratio, MSI, multi-regional spatial interaction; ROC, receiver operating characteristic.

Further evaluation using confusion-matrix–derived metrics showed that the MSI-based model outperformed radiomics and combined-feature models in terms of accuracy, sensitivity, specificity, PPV, and NPV (Table 3; Figure S4). Calibration curves demonstrated good agreement between predicted and observed mortality probabilities for both MSI- and radiomics-based models (Figure 5B). Decision curve analysis indicated that the MSI-based model provided a higher net clinical benefit across a wide range of threshold probabilities in the test set (Figure 5C).

Based on the predicted risk scores derived from the MSI-based model, patients were stratified into high- and low-risk groups. Kaplan-Meier survival analysis demonstrated a significant separation between these two groups (log-rank P<0.05), with a hazard ratio of 5.13 (95% confidence interval: 2.71–9.71), indicating a substantially increased mortality risk in the high-risk group (Figure 5D). In addition, survival modeling further supported the prognostic value of MSI features. The MSI-based survival model achieved a concordance index of 0.71 in the test set, outperforming the radiomics-based model (C-index =0.68) and the combined-feature model (C-index =0.67). These results indicate that MSI-based habitat features provide more informative and robust prognostic information for mortality prediction in NPC patients.

Predictive performance of the MSI model across different T-stage subgroups

We evaluated the performance of the MSI-based mortality prediction model in the test set by stratifying patients according to T stage. As shown in Table 4, the model achieved accuracy values of 0.73, 0.90, 0.66, and 0.81 for mortality prediction in T1-, T2-, T3-, and T4-stage NPC patients, respectively.

Discussion

Accurate mortality risk prediction is essential for effective risk stratification and treatment planning in patients with NPC. In this study, we demonstrate that an MRI-based habitat-driven MSI model provides superior performance for mortality prediction compared with conventional whole-tumor radiomics. The proposed MSI model achieved a higher AUC on the independent test set and consistently greater clinical net benefit on decision curve analysis, indicating that explicitly modeling intratumoral spatial interactions yields prognostically meaningful information beyond global radiomic descriptors.

Habitat-based imaging has been widely investigated as a strategy to characterize intratumoral heterogeneity by partitioning tumors into subregions with distinct imaging phenotypes. Previous studies have demonstrated the prognostic relevance of voxel-level habitat analysis across multiple tumor types using both superpixel-based and intensity-driven clustering approaches (28,32,33). These studies collectively support the premise that subregional imaging phenotypes capture clinically meaningful heterogeneity that is not adequately reflected by whole-tumor analysis alone. However, most existing approaches primarily focus on habitat delineation or subregional feature extraction, without explicitly considering the spatial interactions between different habitats.

By leveraging spatial intensity variation and local gradient information from T2W-MRI, the proposed habitat maps capture both global signal heterogeneity and local spatial transitions within NPC tumors. This enables an intuitive and interpretable visualization of intratumoral spatial heterogeneity, allowing distinct tumor subregions and their interfaces to be explicitly identified. Such spatial information provides insights into heterogeneous tumor composition and organization, which may reflect differences in tumor aggressiveness, invasiveness, and microenvironmental interactions that are relevant to patient prognosis.

Notably, MSI features describing both subregional volume and inter-regional boundary interactions demonstrated significant associations with mortality risk. In particular, larger Volume Class 4 and increased Boundary Volume Class 3–4 were associated with a higher risk of death. These findings suggest that not only the extent of aggressive tumor habitats, but also the degree of spatial interaction and interface between heterogeneous subregions, may reflect invasive tumor behavior and adverse biological processes. Importantly, such spatial interaction patterns are inherently overlooked by conventional radiomics approaches, which aggregate features across the entire tumor volume without accounting for inter-regional relationships.

Previous studies have explored diverse post-clustering analytical strategies to elevate habitat information to patient-level biomarkers, including correlation with histopathology (21), survival modeling using Cox proportional hazards analysis (34), and extraction of ecological diversity or subregional radiomics features for outcome prediction (35). These works demonstrate that habitat-derived subregional information can be meaningfully linked to clinical outcomes. Building upon this foundation, our study extends prior work by explicitly modeling spatial interactions among habitat subregions, thereby providing a complementary and clinically relevant perspective on intratumoral heterogeneity.

Compared with conventional radiomics, which summarizes imaging features at the whole-tumor level, the MSI framework captures both intra-habitat characteristics and inter-habitat spatial relationships. From a clinical perspective, the MSI-based model exhibited a recall-oriented performance profile, suggesting an enhanced ability to identify patients at high risk of death. In prognostic applications, particularly for aggressive malignancies such as NPC, higher sensitivity is often considered clinically advantageous, as it may reduce the likelihood of underestimating patients who could benefit from closer surveillance or intensified treatment strategies. Although specificity and precision were slightly lower than those of certain baseline classifiers, the improved recall highlights the potential utility of MSI features for risk stratification.

Furthermore, subgroup analysis across different T-stage categories demonstrated that the MSI-based model maintained consistent discriminatory performance in patients with varying local tumor extents. This finding suggests that MSI features may capture prognostically relevant tumor heterogeneity that is not fully reflected by conventional anatomic staging alone. Accordingly, the proposed MSI framework may provide complementary prognostic information beyond traditional T-stage classification. Such a more comprehensive representation of intratumoral heterogeneity may partly account for the superior performance of the MSI model in mortality prediction and risk stratification observed in this study. By enabling more accurate identification of high-risk patients, this approach may help inform individualized clinical risk assessment and management strategies.

Beyond NPC, the proposed habitat-based MSI framework may be applicable to other solid tumors characterized by pronounced intratumoral heterogeneity, such as glioblastoma. By partitioning tumors into phenotypically distinct subregions and explicitly modeling their spatial interactions, MSI features provide a structured representation of intratumoral organization that extends beyond conventional whole-tumor radiomic descriptors. Importantly, such region-level spatial interaction information may serve as a meaningful interface between imaging-derived biomarkers and geometry-aware or feature-informed computational tumor modeling frameworks, where spatial priors or inter-regional coupling are required to constrain tumor growth or invasion dynamics. Recent studies (36) have highlighted the value of integrating imaging-informed regional features into computational modeling of tumor behavior, supporting the feasibility of such cross-scale integration. Although computational or biomechanical modeling was beyond the scope of the present study, the proposed framework offers a scalable, imaging-driven foundation for future integration with mechanistic tumor modeling approaches.

This study has several limitations. First, the analysis was conducted using T2-weighted NPC MRI data from a single center, which may limit the generalizability of the findings. To more comprehensively validate the prognostic capability of the proposed models, future work should incorporate larger multi-center cohorts. Second, due to the relatively limited number of deceased cases, a stratified resampling strategy was adopted during model development, which represents an inherent trade-off in retrospective studies with rare outcome events and warrants further validation in larger datasets. Third, the absence of corresponding clinical variables constrained the performance of the combined-feature models. Integrating relevant clinical features in future studies may substantially enhance predictive accuracy. Lastly, due to the lack of pathological data and high-resolution anatomical imaging, the biological basis and clinical relevance of key MSI features—such as Boundary Volume Class 3–4 and Volume Class 4—could not be further elucidated. Future investigations incorporating pathological and structural imaging information will be necessary to clarify their underlying mechanisms and applicability.

Conclusions

This study investigated an MRI-based habitat analysis approach for quantifying intratumoral spatial heterogeneity in NPC and systematically compared it with conventional radiomics features. Our findings demonstrate that partitioning the tumor into subregions and extracting MSI features provides a more precise characterization of intratumoral complexity, thereby enhancing the predictive performance for mortality risk and enabling more accurate stratification of high- and low-risk patients. Among the selected features, Volume Class 4 and Boundary Volume Class 3–4 showed particularly strong prognostic value and may serve as potential imaging biomarkers for mortality risk assessment in NPC. Beyond this cancer type, the proposed habitat-analysis modeling framework has potential applicability across different malignancies and merits further validation in other cancer cohorts.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://qims.amegroups.com/article/view/10.21037/qims-2025-1-2551/rc

Data Sharing Statement: Available at https://qims.amegroups.com/article/view/10.21037/qims-2025-1-2551/dss

Funding: The study was partially supported by project grants from the National Natural Science Foundation of China (Nos. 82260358 and 82171906) and Guangdong Basic and Applied Basic Research Foundation (No. 2025A1515011590).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://qims.amegroups.com/article/view/10.21037/qims-2025-1-2551/coif). All authors report that the present work was partially supported by the National Natural Science Foundation of China (Nos. 82260358 and 82171906) and the Guangdong Basic and Applied Basic Research Foundation (No. 2025A1515011590). P.L. has a pending Chinese patent related to the research topic, with no commercial licensing or financial benefit. Shuchao Chen, Shu Chen, M.L. and H.C. are co-inventors of a pending Chinese patent related to the research topic, with no commercial licensing or financial benefit. The authors have no other 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. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the Ethics Committee of Sun Yat-sen University Cancer Center (No. B2019-222-01). The requirement for informed consent was waived due to the retrospective design of the 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. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, Bray F. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin 2021;71:209-49. [Crossref] [PubMed]
2. Sun XS, Liang YJ, Liu SL, Li XY, Chen QY, Guo SS, Wen YF, Liu LT, Xie HJ, Tang QN, Yan JJ, Guo L, Ma J, Tang LQ, Mai HQ. Establishment and validation of a nomogram for predicting survival in patients with de novo metastatic nasopharyngeal carcinoma. Oral Oncol 2019;94:73-9. [Crossref] [PubMed]
3. Lian CL, Zhou R, Zhou Y, Zhou P, Wu SG. Assessment of Response to Different Induction Chemotherapy Regimens in Locally Advanced Nasopharyngeal Carcinoma. Drug Des Devel Ther 2023;17:551-62. [Crossref] [PubMed]
4. Eisenstein M. Why tumour geography matters - and how to map it. Nature 2024;635:1031-3. [Crossref] [PubMed]
5. Amin MB, Greene FL, Edge SB, Compton CC, Gershenwald JE, Brookland RK, Meyer L, Gress DM, Byrd DR, Winchester DP. The Eighth Edition AJCC Cancer Staging Manual: Continuing to build a bridge from a population-based to a more “personalized” approach to cancer staging. CA Cancer J Clin 2017;67:93-9.
6. Liao XB, Mao YP, Liu LZ, Tang LL, Sun Y, Wang Y, Lin AH, Cui CY, Li L, Ma J. How does magnetic resonance imaging influence staging according to AJCC staging system for nasopharyngeal carcinoma compared with computed tomography? Int J Radiat Oncol Biol Phys 2008;72:1368-77. [Crossref] [PubMed]
7. Junn JC, Soderlund KA, Glastonbury CM. Imaging of Head and Neck Cancer With CT, MRI, and US. Semin Nucl Med 2021;51:3-12. [Crossref] [PubMed]
8. Cheng J, Su W, Wang Y, Zhan Y, Wang Y, Yan S, Yuan Y, Chen L, Wei Z, Zhang S, Gao X, Tang Z. Magnetic resonance imaging based on radiomics for differentiating T1-category nasopharyngeal carcinoma from nasopharyngeal lymphoid hyperplasia: a multicenter study. Jpn J Radiol 2024;42:709-19. [Crossref] [PubMed]
9. Wong LM, Ai QYH, Zhang R, Mo F, King AD. Radiomics for Discrimination between Early-Stage Nasopharyngeal Carcinoma and Benign Hyperplasia with Stable Feature Selection on MRI. Cancers (Basel) 2022;14:3433. [Crossref] [PubMed]
10. Li Q, Yu Q, Gong B, Ning Y, Chen X, Gu J, Lv F, Peng J, Luo T. The Effect of Magnetic Resonance Imaging Based Radiomics Models in Discriminating stage I-II and III-IVa Nasopharyngeal Carcinoma. Diagnostics (Basel) 2023;13:300. [Crossref] [PubMed]
11. Zhao L, Gong J, Xi Y, Xu M, Li C, Kang X, Yin Y, Qin W, Yin H, Shi M. MRI-based radiomics nomogram may predict the response to induction chemotherapy and survival in locally advanced nasopharyngeal carcinoma. Eur Radiol 2020;30:537-46. [Crossref] [PubMed]
12. Hou J, Li H, Zeng B, Pang P, Ai Z, Li F, Lu Q, Yu X. MRI-based radiomics nomogram for predicting temporal lobe injury after radiotherapy in nasopharyngeal carcinoma. Eur Radiol 2022;32:1106-14. [Crossref] [PubMed]
13. Jiang S, Han L, Liang L, Long L. Development and validation of an MRI-based radiomic model for predicting overall survival in nasopharyngeal carcinoma patients with local residual tumors after intensity-modulated radiotherapy. BMC Med Imaging 2022;22:174. [Crossref] [PubMed]
14. Raghavan Nair JK, Vallières M, Mascarella MA, El Sabbagh N, Duchatellier CF, Zeitouni A, Shenouda G, Chankowsky J. Magnetic Resonance Imaging Texture Analysis Predicts Recurrence in Patients With Nasopharyngeal Carcinoma. Can Assoc Radiol J 2019;70:394-402. [Crossref] [PubMed]
15. Cao X, Wang X, Song J, Su Y, Wang L, Yin Y. Pretreatment multiparametric MRI radiomics-integrated clinical hematological biomarkers can predict early rapid metastasis in patients with nasopharyngeal carcinoma. BMC Cancer 2024;24:435. [Crossref] [PubMed]
16. Wu LX, Ding N, Ji YD, Zhang YC, Li MJ, Shen JC, Hu HT, Jin L, Yin SN. Habitat Analysis in Tumor Imaging: Advancing Precision Medicine Through Radiomic Subregion Segmentation. Cancer Manag Res 2025;17:731-41. [Crossref] [PubMed]
17. Liu Y, He S, Wang XL, Peng W, Chen QY, Chi DM, et al. Tumour heterogeneity and intercellular networks of nasopharyngeal carcinoma at single cell resolution. Nat Commun 2021;12:741. [Crossref] [PubMed]
18. Tan AC, Ashley DM, López GY, Malinzak M, Friedman HS, Khasraw M. Management of glioblastoma: State of the art and future directions. CA Cancer J Clin 2020;70:299-312. [Crossref] [PubMed]
19. Qiao J, Kang H, Ran Q, Tong H, Ma Q, Wang S, Zhang W, Wu H. Metabolic habitat imaging with hemodynamic heterogeneity predicts individual progression-free survival in high-grade glioma. Clin Radiol 2024;79:e842-53. [Crossref] [PubMed]
20. Jeong SY, Park JE, Kim N, Kim HS. Hypovascular Cellular Tumor in Primary Central Nervous System Lymphoma is Associated with Treatment Resistance: Tumor Habitat Analysis Using Physiologic MRI. AJNR Am J Neuroradiol 2022;43:40-7. [Crossref] [PubMed]
21. Lee DH, Park JE, Kim N, Park SY, Kim YH, Cho YH, Kim HS. Tumor habitat analysis by magnetic resonance imaging distinguishes tumor progression from radiation necrosis in brain metastases after stereotactic radiosurgery. Eur Radiol 2022;32:497-507. [Crossref] [PubMed]
22. Cho HH, Kim H, Nam SY, Lee JE, Han BK, Ko EY, Choi JS, Park H, Ko ES. Measurement of Perfusion Heterogeneity within Tumor Habitats on Magnetic Resonance Imaging and Its Association with Prognosis in Breast Cancer Patients. Cancers (Basel) 2022;14:1858. [Crossref] [PubMed]
23. Xu C, Wang Z, Wang A, Zheng Y, Song Y, Wang C, Yang G, Ma M, He M. Breast Cancer: Multi-b-Value Diffusion Weighted Habitat Imaging in Predicting Pathologic Complete Response to Neoadjuvant Chemotherapy. Acad Radiol 2024;31:4733-42. [Crossref] [PubMed]
24. Wang C, Wu F, Wang F, Chong HH, Sun H, Huang P, Xiao Y, Yang C, Zeng M. The Association Between Tumor Radiomic Analysis and Peritumor Habitat-Derived Radiomic Analysis on Gadoxetate Disodium-Enhanced MRI With Microvascular Invasion in Hepatocellular Carcinoma. J Magn Reson Imaging 2025;61:1428-39. [Crossref] [PubMed]
25. Huang F, Huang Q, Liao X, Gao Y. Prediction of high-risk prostate cancer based on the habitat features of biparametric magnetic resonance and the omics features of contrast-enhanced ultrasound. Heliyon 2024;10:e37955. [Crossref] [PubMed]
26. Zhu Y, Zheng D, Xu S, Chen J, Wen L, Zhang Z, Ruan H. Intratumoral habitat radiomics based on magnetic resonance imaging for preoperative prediction treatment response to neoadjuvant chemotherapy in nasopharyngeal carcinoma. Jpn J Radiol 2024;42:1413-24. [Crossref] [PubMed]
27. Wu J, Cao G, Sun X, Lee J, Rubin DL, Napel S, Kurian AW, Daniel BL, Li R. Intratumoral Spatial Heterogeneity at Perfusion MR Imaging Predicts Recurrence-free Survival in Locally Advanced Breast Cancer Treated with Neoadjuvant Chemotherapy. Radiology 2018;288:26-35. [Crossref] [PubMed]
28. Zhu S, Li S, Cao D, Luo C, Liang Z, Liang S, Zhang G, Zhao Q, Ruan G, Liu L, Fu G, Li H. Prognostic Value of Skull Base Foramen Invasion Subclassification in T Category Modification and Induction Chemotherapy Management for Nasopharyngeal Carcinoma: Post-Hoc Analysis of a Dual-Center Retrospective Cohort Study. Adv Sci (Weinh) 2025;12:e2408182. [Crossref] [PubMed]
29. Sadaiyandi J, Arumugam P, Sangaiah AK, Zhang C. Stratified Sampling-Based Deep Learning Approach to Increase Prediction Accuracy of Unbalanced Dataset. Electronics 2023;12:4423.
30. Likas A, Vlassis N, Verbeek JJ. The global k-means clustering algorithm. Pattern Recognition 2003;36:451-61.
31. Liu Y, Li Z, Xiong H, Gao X, Wu J, Wu S. Understanding and enhancement of internal clustering validation measures. IEEE Trans Cybern 2013;43:982-94. [Crossref] [PubMed]
32. Kononenko I. Machine learning for medical diagnosis: history, state of the art and perspective. Artif Intell Med 2001;23:89-109. [Crossref] [PubMed]
33. Liang JD, Ping XO, Tseng YJ, Huang GT, Lai F, Yang PM. Recurrence predictive models for patients with hepatocellular carcinoma after radiofrequency ablation using support vector machines with feature selection methods. Comput Methods Programs Biomed 2014;117:425-34. [Crossref] [PubMed]
34. Xie C, Yang P, Zhang X, Xu L, Wang X, Li X, Zhang L, Xie R, Yang L, Jing Z, Zhang H, Ding L, Kuang Y, Niu T, Wu S. Sub-region based radiomics analysis for survival prediction in oesophageal tumours treated by definitive concurrent chemoradiotherapy. EBioMedicine 2019;44:289-97. [Crossref] [PubMed]
35. Bailo M, Pecco N, Callea M, Scifo P, Gagliardi F, Presotto L, Bettinardi V, Fallanca F, Mapelli P, Gianolli L, Doglioni C, Anzalone N, Picchio M, Mortini P, Falini A, Castellano A. Decoding the Heterogeneity of Malignant Gliomas by PET and MRI for Spatial Habitat Analysis of Hypoxia, Perfusion, and Diffusion Imaging: A Preliminary Study. Front Neurosci 2022;16:885291. [Crossref] [PubMed]
36. Ghahramani MR, Bavi O. Heterogeneous biomechanical/mathematical modeling of spatial prediction of glioblastoma progression using magnetic resonance imaging-based finite element method. Comput Methods Programs Biomed 2024;257:108441. [Crossref] [PubMed]