Development and application of a novel tumor habitat analysis technique based on dynamical modeling

Oropharyngeal cancer (OPC) exhibits varying responses to chemoradiation therapy, making treatment outcome prediction challenging. Traditional imaging-based methods often fail to capture the spatial heterogeneity within tumors, which influences treatment resistance and disease progression. Advances in modeling techniques allow for more nuanced analysis of this heterogeneity, identifying distinct tumor regions, or habitats, that drive patient outcomes.

To interrogate the association between treatment-induced changes in spatial heterogeneity and chemoradiation resistance of oropharyngeal cancer (OPC) based on a novel tumor habitat analysis.

A mathematical model was used to estimate tumor time dynamics of patients with OPC based on the applied analysis of partial differential equations. The position and momentum of each voxel was propagated according to Fokker-Planck dynamics, that is, a common model in statistical mechanics. The boundary conditions of the Fokker-Planck equation were solved based on pre- and intra-treatment (i.e., after 2 weeks of therapy) ¹⁸F-FDG-PET SUV images of patients (n = 56) undergoing definitive (chemo)radiation for OPC as part of a previously conducted prospective clinical trial. Tumor-specific time dynamics, measured based on the solution of the Fokker-Planck equation, were generated for each patient. Tumor habitats (i.e., non-overlapping subregions of the primary tumor) were identified by measuring vector similarity in voxel-level time dynamics through a fuzzy c-means clustering algorithm. The robustness of our habitat construction method was quantified using a mean silhouette metric to measure intra-habitat variability. Fifty-four habitat-specific radiomic texture features were extracted from pre-treatment SUV images and normalized by habitat volume. Univariate Kaplan-Meier analyses were implemented as a feature selection method, where statistically significant features (p < 0.05, log-rank) were used to construct a multivariate Cox proportional-hazards model. Parameters from the resulting Cox model were then used to construct a risk score for each patient, based on habitat-specific radiomic expression. The patient cohort was stratified by median risk score value and association with recurrence-free survival (RFS) was evaluated via log-rank tests.

Dynamic tumor habitat analysis partitioned the gross disease of each patient into three spatial subregions. Voxels within each habitat suggested differential response rates in different compartments of the tumor. The minimum mean silhouette value was 0.57 and maximum mean silhouette value was 0.8, where values above 0.7 indicated strong intra-habitat consistency and values between 0.5 and 0.7 indicated reasonable intra-habitat consistency. Nine radiomic texture features (three GLRLM, two GLCOM, and three GLSZM) and SUVmax were found to be prognostically significant and were used to build the multivariate Cox model. The resulting risk score was associated with RFS (p = 0.032). By contrast, potential confounding factors (primary tumor volume and mean SUV) were not significantly associated with RFS (p = 0.286 and p = 0.231, respectively).

We interrogated spatial heterogeneity of oropharyngeal tumors through the application of a novel algorithm to identify spatial habitats on SUV images. Our habitat construction technique was shown to be robust and habitat-specific feature spaces revealed distinct underlying radiomic expression patterns. Radiomic features were extracted from dynamic habitats and used to build a risk score which demonstrated prognostic value.