Methods and systems for predicting neoadjuvant chemoradiation of tumors based on multi-dimensional habitat imaging features

Abstract

This invention discloses a method and system for predicting neoadjuvant chemoradiotherapy (nCRT) in tumors based on multidimensional habitat imaging features. The method includes: acquiring DWI-MRI images of LARC patients after nCRT; segmenting the tumor region and dividing it into multiple habitat subregions with different imaging phenotypes using an unsupervised clustering algorithm; extracting high-dimensional radiomics features from the entire tumor and each habitat subregion; constructing a joint prediction model integrating clinicopathological features, overall radiomics features, and intratumoral ecological diversity features of habitat subregions; and using this model to predict whether the patient will achieve total tumor density (TDS) after nCRT, outputting the prediction probability and classification results. This invention provides a non-invasive, high-precision efficacy prediction tool by visualizing and quantifying the heterogeneity within the tumor, which can be used to guide individualized treatment decisions, such as screening patients suitable for a "wait-and-see" strategy.

Description

Technical Field

[0001] This invention belongs to the fields of medical image processing, artificial intelligence, and precision oncology, and relates to a method and system for predicting neoadjuvant chemoradiotherapy based on multi-dimensional habitat imaging features. Specifically, it involves a method for extracting and analyzing visualization features of intratumoral heterogeneity (ITH) based on diffusion-weighted magnetic resonance imaging (DWI-MRI) habitat imaging, and applying it to build a predictive model to evaluate the tumor downstaging (TSD) effect in patients with locally advanced rectal cancer (LARC) after receiving neoadjuvant chemoradiotherapy (nCRT). Background Technology

[0002] Rectal cancer is one of the most common malignant tumors worldwide. For locally advanced rectal cancer (LARC), the standard treatment is neoadjuvant chemoradiotherapy (nCRT) combined with total mesorectal excision. However, patient responses to nCRT are highly heterogeneous; some patients achieve significant pathological complete remission or tumor downstaging, thus having the opportunity to adopt a "wait and see" strategy or undergo local excision to preserve anal function and significantly improve their quality of life. Therefore, accurately predicting the efficacy of nCRT preoperatively is crucial for developing individualized treatment decisions.

[0003] Currently, clinical practice primarily relies on postoperative pathological examinations (such as tumor regression grading, TRG) to assess treatment efficacy, which is an invasive and time-consuming method. Preoperative prediction is mainly based on clinical characteristics and routine imaging assessments, but its accuracy is limited. MRI, especially DWI sequences, can non-invasively reflect the cell density and water molecule diffusion restriction of tumor tissue, and has become an important imaging tool for assessing the response to rectal cancer treatment.

[0004] Radiomics provides a new dimension for tumor phenotyping by extracting quantitative features from medical images in high throughput. However, traditional radiomics methods typically analyze the entire tumor as a homogeneous region, ignoring the spatial heterogeneity that exists within the tumor. This intratumoral heterogeneity (ITH) is an important biological basis for treatment resistance and prognostic differences.

[0005] Habitat imaging is an emerging image analysis method that draws on ecological concepts to divide subregions within a tumor that exhibit different phenotypes on imaging into distinct "habitats," thereby visualizing and quantifying intratumoral hematoxylin and eosinophil (ITH). This method can reveal more nuanced differences in biological behavior within the tumor.

[0006] Currently, there is a lack of a comprehensive technical solution that can effectively integrate DWI-MRI imaging information, quantify intratumoral spatial heterogeneity, and construct a high-precision predictive model for preoperative non-invasive assessment of TDS after nCRT in LARC patients. Therefore, developing an ITH analysis and prediction method based on DWI-MRI habitat imaging has significant clinical needs and application value.

[0007] The aforementioned existing technologies have the following drawbacks in predicting the efficacy of nCRT in LARC patients: 1. Insufficient utilization of spatial heterogeneity information within tumors; traditional holistic analysis methods cannot capture key local features that affect treatment response.

[0008] 2. There is a lack of predictive model frameworks that can effectively integrate advanced radiomics features with clinicopathological information.

[0009] 3. The accuracy, robustness, and clinical interpretability of existing predictive models need further improvement, making it difficult to reliably guide individualized treatment decisions. Summary of the Invention

[0010] The present invention aims to overcome the above-mentioned shortcomings of the prior art and provides a method for predicting tumor downstaging after neoadjuvant chemoradiotherapy for rectal cancer based on DWI-MRI habitat imaging.

[0011] The core of this invention lies in visualizing and quantifying the heterogeneity within tumors through habitat imaging technology, and constructing a joint prediction model that integrates multi-dimensional features.

[0012] The first aspect of the present invention relates to a method for predicting neoadjuvant chemoradiotherapy for tumors based on multidimensional habitat imaging features, comprising the following steps: Step 1. Data Acquisition and Preprocessing: Acquire diffusion-weighted magnetic resonance imaging (DWI-MRI) data and corresponding clinicopathological data of patients with locally advanced rectal cancer after neoadjuvant chemoradiotherapy (nCRT); preprocess the DWI-MRI data. Step 2. Build a prediction model. The model architecture consists of the following modules: tumor segmentation and habitat subregion division module, multi-dimensional feature extraction module, and feature selection and training module.

[0013] The tumor segmentation and habitat subregion division module refers to manually or automatically delineating the tumor region as the region of interest (ROI) on the preprocessed DWI-MRI image, and using unsupervised clustering to divide the tumor into K visually distinguishable habitat subregions based on the ROI. The multi-dimensional feature extraction module refers to extracting not only image omics features from the divided habitat sub-regions, but also advanced features such as sub-region topological structure features and sub-region spatial invasiveness features. The feature selection and training module involves initially screening the features extracted by the multi-dimensional feature extraction module using univariate analysis, then using the mRMR algorithm for dimensionality reduction to form a joint feature set. Subsequently, a genetic algorithm is used to further screen and optimize the joint feature set and clinical features. First, a chromosome population containing fixed clinical features is initialized. The XGBoost classifier is trained on the training set and its performance is evaluated on the validation set, with the area under the receiver operating characteristic curve (AUC) used as the primary indicator to calculate fitness. Next, genetic operations are performed, updating the population through tournament selection, two-point crossover, and position flip mutation, while protecting the fixed clinical features from being altered. Then, convergence is monitored, and optimization is terminated when the maximum number of iterations or the fitness plateau is reached. Finally, the feature subset with the best overall performance on the validation set is output, and XGBoost is retrained using the best subset to obtain the final model.

[0014] Step 3. Prediction and Output: Input the joint feature set of the patient to be predicted into the trained prediction model to obtain the predicted probability or classification label of the patient reaching TDS, and output the prediction result. The trained model is then transformed into a stable and interpretable clinical decision support tool to achieve automated prediction and result output for new patients.

[0015] Furthermore, step 1 specifically includes: Step 101: Data Acquisition. Collect multiparametric MRI images of the pelvis obtained from patients with pathologically confirmed locally advanced rectal cancer (LARC) after neoadjuvant chemoradiotherapy (nCRT) and before surgery (usually 6-8 weeks after radiotherapy). Step 102: Bias field correction, commonly using the N4 algorithm, which iteratively solves for the bias field through maximum likelihood estimation. B(x)B(x) Step 103: Image registration. Spatially align images from different sequences (e.g., T2WI and DWI) or different time points from the same patient to ensure consistency in subsequent analysis areas. Use rigid or affine transformations, with the high-resolution T2WI image as a fixed reference, to register the DWI image with it.

[0016] Step 104: Resampling. Standardize the image voxel size to eliminate anisotropic resolution differences caused by different scanning protocols. Use linear interpolation or B-spline interpolation to resample all images to isotropic voxels. The interpolation formula is as follows: (1,1) in The grayscale value of the target image. The voxel values ​​of the original image Step 105: Intensity standardization, to make the grayscale values ​​of images from different scanning devices and different patients comparable. Z-score standardization is calculated using the following formula: (1,2) in m and s It is the mean and standard deviation of an image or a specific region. Furthermore, step 2 specifically includes: Step 201: First, a simple linear iterative clustering algorithm is used to segment the tumor region into multiple compact, uniform three-dimensional hypervoxels; then, a K-means clustering algorithm is used to cluster the hypervoxel feature vector sets to discover different image phenotype "habitats". N Each hypervoxon is divided into K Within each cluster, the intra-class squared error is minimized; the number of habitat subregions K is optimized using at least one evaluation index among the silhouette coefficient, Calinski-Harabasz index, and Davies-Bouldin index, and the optimal value is determined by combining the elbow rule with the silhouette coefficient. K value.

[0017] Step 202: Extract radiomics features from each of the identified habitat subregions. These features include at least the following three categories: Shape characteristics: Describe the three-dimensional geometry of the tumor or its subregions, such as volume, surface area, sphericity, surface area to volume ratio, and maximum three-dimensional diameter.

[0018] First-order statistical characteristics: describe the distribution of voxel intensity within a region, such as mean, median, standard deviation, skewness, kurtosis, energy, entropy, etc.

[0019] Texture features: describe the spatial arrangement and interrelationship of voxel intensity, calculated by methods such as gray-level co-occurrence matrix, gray-level run-length matrix, gray-level size region matrix, and neighborhood gray-level difference matrix, such as contrast, correlation, homogeneity, entropy, and short run-length advantage.

[0020] Step 203: To further quantify the spatial configuration and invasive potential of intratumoral heterogeneity, subregional spatial topological features are extracted, including: 1. Dominant Subregion Analysis Characteristics: Calculate the volume percentage of each subregion relative to the total tumor volume to identify dominant subregions; 2. Spatial distribution characteristics of subregions: By calculating the spatial contact surface area and voxel adjacency relationship between different subregions, or by using spatial correlation analysis, it can be determined whether a specific subregion tends to be adjacent or separated in space; The formula for calculating the spatial relationship of subregions is as follows: (2,1) in and grayscale value of a pixel These features are designed to depict the full picture of response heterogeneity within a tumor, quantifying whether a feature is spatially clustered, discrete, or randomly distributed, rather than viewing a single subregion in isolation.

[0021] Step 204: To focus on the morphological complexity of the tumor-normal tissue boundary or subregion margins to assess their invasiveness, subregion spatial invasiveness features are extracted, including: 1. Edge wrinkling degree characteristics, calculated using the following formula: (2,2) It is the size of pixels. N Minimum number required for coverage 2. Spatial regularity characteristics, calculated using the following formula: (2,3) V For volume, A This represents the surface area. Values ​​closer to 1 indicate a more regular shape, while smaller values ​​indicate a more irregular shape. 3. Interface sharpness characteristics: Assess the steepness of the gradient change in image signal intensity between a specific subregion and its adjacent normal tissue or other subregions; a blurred interface suggests invasive growth. At the boundary between the subregion and the adjacent normal tissue, extract a narrow band along the normal direction and calculate the average gradient amplitude of the DWI signal within this narrow band. The calculation formula is as follows: (2,4) Band is a set of pixels, and p is a single pixel. N The total number of pixels, the lower G The value indicates blurred boundaries, suggesting invasive growth under a microscope. Step 205: Extraction of clinicopathological features, including at least one of the following: modified tumor regression grade (mTRG), postoperative N stage, carcinoembryonic antigen (CEA) level, and carbohydrate antigen (CA) 19-9 level; Step 206: Standardize and select features from the extracted high-dimensional features. First, univariate analysis (t-test / Mann-Whitney U test) is used to screen for features that show significant differences between the TDS and non-TDS groups (p<0.05). Then, the mRMR algorithm is used for further dimensionality reduction, and finally, about 200 of the most predictive features are retained to form a joint feature set.

[0022] Step 207: Subsequently, a genetic algorithm is used to screen and optimize the joint feature set and clinical features. A binary chromosome population containing fixed clinical features is initialized, with a population size of 200. For each individual's feature subset, an XGBoost classifier is trained on the training set and its performance is evaluated on the validation set. The fitness value is calculated using the area under the receiver operating characteristic curve as the core indicator, while a feature quantity penalty term is introduced to control model complexity. Then, genetic operations are performed, using a tournament selection strategy to select superior individuals. A new generation of population is generated using two-point crossover and position flip mutation operators, ensuring that the clinical feature gene positions are not changed during the crossover and mutation process, and supplemented by an elite retention mechanism to maintain the optimal solution of the population. The convergence status of the evolutionary process is continuously monitored, and optimization is terminated when the preset maximum number of iterations is reached or the fitness does not significantly improve for several consecutive generations. Finally, the feature subset with the best overall performance on the validation set is output, providing optimized feature input for subsequent model construction.

[0023] Furthermore, step 3 specifically includes: Step 301: Use the area under the receiver operating characteristic curve (AUC), calibration curve, and decision curve analysis (DCA) to comprehensively evaluate the model's discrimination, calibration, and clinical usefulness.

[0024] Step 302: Save the finally trained optimal model (such as the XGBoost model) in a standardized format (such as a Python pickle file, ONNX format, or PMML format) to ensure that it can be stably loaded and used in a production environment.

[0025] Based on the same inventive concept, the present invention also provides a system for implementing the above method, comprising: Data Input and Management Module: Used to import and store patients' DWI-MRI images and clinical data.

[0026] Image processing and analysis module: used to perform image preprocessing, tumor segmentation, habitat classification and radiomics feature extraction.

[0027] Model computation module: Used to store pre-trained prediction models and perform feature computation and classification prediction on new patient data.

[0028] Results visualization and reporting module: used to display tumor habitat segmentation maps, feature importance maps, model prediction probabilities, and clinical recommendation reports.

[0029] User interface: Provides doctors with access to the operating system, view results, and adjust parameters.

[0030] Compared with the prior art, the present invention has the following significant advantages: 1. Innovative quantification of intratumoral heterogeneity: For the first time, the concept of habitat imaging was systematically applied to DWI-MRI analysis of LARC. Through unsupervised clustering, the tumor was visualized and segmented into different phenotypic subregions, realizing a quantitative description of intratumoral spatial heterogeneity and capturing key biological information affecting the treatment effect.

[0031] 2. Constructing a high-precision prediction model: By integrating habitat ITH features, overall radiomics features, and clinicopathological features, if the constructed joint prediction model shows excellent prediction performance on the independent test set, significantly outperforming the single feature model, then this will provide a more reliable basis for clinical decision-making.

[0032] 3. Enhancing Clinical Translational Value: This invention provides a complete end-to-end solution, from raw images to clinical prediction results. The model, validated through decision curve analysis, demonstrates clear clinical net benefits and effectively assists physicians in selecting advantageous patients suitable for a "wait-and-see" strategy, promoting individualized and precise treatment of rectal cancer.

[0033] 4. Good interpretability: Through feature importance analysis (such as SHAP), the imaging features that contribute the most to the prediction and their corresponding tumor subregions can be identified, which enhances the credibility of the model and helps to understand the biological significance behind the imaging features.

[0034] 5. Achieved in-depth quantification of tumor spatial biological behavior: By extracting subregional spatial topology and edge invasiveness features, this invention not only describes what a tumor "is", but also reveals how its different internal components "are spatially arranged" and "how they grow invasively", establishing a more direct and interpretable correlation between imaging features and the malignant biological behavior of tumors (such as invasion and treatment resistance), surpassing the dimensions of traditional radiomics analysis. Attached Figure Description

[0035] Figure 1 This is a flowchart illustrating the overall process of the method described in the embodiments of the present invention, showing the overall processing flow of steps 1 to 3. Detailed Implementation

[0036] The present invention will be further described in detail below with reference to the accompanying drawings and embodiments. The following embodiments are only used to explain the present invention and do not constitute a limitation on the scope of protection of the present invention.

[0037] Example 1

[0038] A method for predicting neoadjuvant chemoradiotherapy for tumors based on multidimensional habitat imaging features includes the following steps: Step 1. Data Acquisition and Preprocessing: Acquire diffusion-weighted magnetic resonance imaging (DWI-MRI) data and corresponding clinicopathological data of patients with locally advanced rectal cancer after neoadjuvant chemoradiotherapy (nCRT); preprocess the DWI-MRI data. Step 2. Build the prediction model. The model architecture consists of the following modules: tumor segmentation and habitat subregion division module, multi-dimensional feature extraction module, and feature selection and training module. The tumor segmentation and habitat subregion division module refers to manually or automatically delineating the tumor region as the region of interest (ROI) on the preprocessed DWI-MRI image, and using unsupervised clustering to divide the tumor into K visually distinguishable habitat subregions based on the ROI. The multi-dimensional feature extraction module refers to extracting not only image omics features from the divided habitat sub-regions, but also advanced features such as sub-region topological structure features and sub-region spatial invasiveness features. The feature selection and training module involves initially screening the features extracted by the multi-dimensional feature extraction module using univariate analysis, then using the mRMR algorithm for dimensionality reduction to form a joint feature set. Subsequently, a genetic algorithm is used to further screen and optimize the joint feature set and clinical features. First, a chromosome population containing fixed clinical features is initialized. The XGBoost classifier is trained on the training set and its performance is evaluated on the validation set, with the area under the receiver operating characteristic curve (AUC) used as the primary indicator to calculate fitness. Next, genetic operations are performed, updating the population through tournament selection, two-point crossover, and position flip mutation, while protecting the fixed clinical features from being altered. Then, convergence is monitored, and optimization is terminated when the maximum number of iterations or the fitness plateau is reached. Finally, the feature subset with the best overall performance on the validation set is output, and XGBoost is retrained using the best subset to obtain the final model.

[0039] Step 3. Prediction and Output: Input the joint feature set of the patient to be predicted into the trained prediction model to obtain the predicted probability or classification label of the patient reaching TDS, and output the prediction result. The trained model is then transformed into a stable and interpretable clinical decision support tool to achieve automated prediction and result output for new patients.

[0040] Furthermore, step 1 specifically includes: Step 101: Data Acquisition. Collect multiparametric MRI images of the pelvis obtained from patients with pathologically confirmed locally advanced rectal cancer (LARC) after neoadjuvant chemoradiotherapy (nCRT) and before surgery (usually 6-8 weeks after radiotherapy). Step 102: Bias field correction, commonly using the N4 algorithm, which iteratively solves for the bias field through maximum likelihood estimation. B(x)B(x) Step 103: Image registration. Spatially align images from different sequences (e.g., T2WI and DWI) or different time points from the same patient to ensure consistency in subsequent analysis areas. Use rigid or affine transformations, with the high-resolution T2WI image as a fixed reference, to register the DWI image with it.

[0041] Step 104: Resampling. Standardize the image voxel size to eliminate anisotropic resolution differences caused by different scanning protocols. Use linear interpolation or B-spline interpolation to resample all images to isotropic voxels. The interpolation formula is as follows: (1,1) in The grayscale value of the target image. The voxel values ​​of the original image Step 105: Intensity standardization, to make the grayscale values ​​of images from different scanning devices and different patients comparable. Z-score standardization is calculated using the following formula: (1,2) in m and s It is the mean and standard deviation of an image or a specific region. Furthermore, step 2 specifically includes: Step 201: First, a simple linear iterative clustering algorithm is used to segment the tumor region into multiple compact, uniform three-dimensional hypervoxels; then, a K-means clustering algorithm is used to cluster the hypervoxel feature vector sets to discover different image phenotype "habitats". N Each hypervoxon is divided into K Within each cluster, the intra-class squared error is minimized; the number of habitat subregions K is optimized using at least one evaluation index among the silhouette coefficient, Calinski-Harabasz index, and Davies-Bouldin index, and the optimal value is determined by combining the elbow rule with the silhouette coefficient. K value.

[0042] Step 202: Extract radiomics features from each of the identified habitat subregions. These features include at least the following three categories: Shape characteristics: Describe the three-dimensional geometry of the tumor or its subregions, such as volume, surface area, sphericity, surface area to volume ratio, and maximum three-dimensional diameter.

[0043] First-order statistical characteristics: describe the distribution of voxel intensity within a region, such as mean, median, standard deviation, skewness, kurtosis, energy, entropy, etc.

[0044] Texture features: describe the spatial arrangement and interrelationship of voxel intensity, calculated by methods such as gray-level co-occurrence matrix, gray-level run-length matrix, gray-level size region matrix, and neighborhood gray-level difference matrix, such as contrast, correlation, homogeneity, entropy, and short run-length advantage.

[0045] Step 203: To further quantify the spatial configuration and invasive potential of intratumoral heterogeneity, subregional spatial topological features are extracted, including: 1. Dominant Subregion Analysis Characteristics: Calculate the volume percentage of each subregion relative to the total tumor volume to identify dominant subregions; 2. Spatial distribution characteristics of subregions: By calculating the spatial contact surface area and voxel adjacency relationship between different subregions, or by using spatial correlation analysis, it can be determined whether a specific subregion tends to be adjacent or separated in space; The formula for calculating the spatial relationship of subregions is as follows: (2,1) in and grayscale value of a pixel These features are designed to depict the full picture of response heterogeneity within a tumor, quantifying whether a feature is spatially clustered, discrete, or randomly distributed, rather than viewing a single subregion in isolation.

[0046] Step 204: To focus on the morphological complexity of the tumor-normal tissue boundary or subregion margins to assess their invasiveness, subregion spatial invasiveness features are extracted, including: 1. Edge wrinkling degree characteristics, calculated using the following formula: (2,2) It is the size of pixels. N Minimum number required for coverage 2. Spatial regularity characteristics, calculated using the following formula: (2,3) V For volume, A This represents the surface area. Values ​​closer to 1 indicate a more regular shape, while smaller values ​​indicate a more irregular shape. 3. Interface sharpness characteristics: Assess the steepness of the gradient change in image signal intensity between a specific subregion and its adjacent normal tissue or other subregions; a blurred interface suggests invasive growth. At the boundary between the subregion and the adjacent normal tissue, extract a narrow band along the normal direction and calculate the average gradient amplitude of the DWI signal within this narrow band. The calculation formula is as follows: (2,4) Band is a set of pixels, and p is a single pixel. N The total number of pixels, the lower G The value indicates blurred boundaries, suggesting invasive growth under a microscope. Step 205: Extraction of clinicopathological features, including at least one of the following: modified tumor regression grade (mTRG), postoperative N stage, carcinoembryonic antigen (CEA) level, and carbohydrate antigen (CA) 19-9 level; Step 206: Standardize and select features from the extracted high-dimensional features. First, univariate analysis (t-test / Mann-Whitney U test) is used to screen for features that show significant differences between the TDS and non-TDS groups (p<0.05). Then, the mRMR algorithm is used for further dimensionality reduction, and finally, about 200 of the most predictive features are retained to form a joint feature set.

[0047] Step 207: Subsequently, a genetic algorithm is used to screen and optimize the joint feature set and clinical features. A binary chromosome population containing fixed clinical features is initialized, with a population size of 200. For each individual's feature subset, an XGBoost classifier is trained on the training set and its performance is evaluated on the validation set. The fitness value is calculated using the area under the receiver operating characteristic curve as the core indicator, while a feature quantity penalty term is introduced to control model complexity. Then, genetic operations are performed, using a tournament selection strategy to select superior individuals. A new generation of population is generated using two-point crossover and position flip mutation operators, ensuring that the clinical feature gene positions are not changed during the crossover and mutation process, and supplemented by an elite retention mechanism to maintain the optimal solution of the population. The convergence status of the evolutionary process is continuously monitored, and optimization is terminated when the preset maximum number of iterations is reached or the fitness does not significantly improve for several consecutive generations. Finally, the feature subset with the best overall performance on the validation set is output, providing optimized feature input for subsequent model construction.

[0048] Furthermore, step 3 specifically includes: Step 301: Use the area under the receiver operating characteristic curve (AUC), calibration curve, and decision curve analysis (DCA) to comprehensively evaluate the model's discrimination, calibration, and clinical usefulness.

[0049] Step 302: Save the final trained optimal model (such as the XGBoost model) in a standardized format (such as a Python pickle file, ONNX format, or PMML format) to ensure that it can be stably loaded and called in a production environment.

[0050] Example 2

[0051] This embodiment relates to a system for implementing the method of predicting neoadjuvant chemoradiotherapy for tumors based on multidimensional habitat imaging features as described in Embodiment 1, comprising: Data Input and Management Module: Used to import and store patients' DWI-MRI images and clinical data.

[0052] Image processing and analysis module: used to perform image preprocessing, tumor segmentation, habitat classification and radiomics feature extraction.

[0053] Model computation module: Used to store pre-trained prediction models and perform feature computation and classification prediction on new patient data.

[0054] Results visualization and reporting module: used to display tumor habitat segmentation maps, feature importance maps, model prediction probabilities, and clinical recommendation reports.

[0055] User interface: Provides doctors with access to the operating system, view results, and adjust parameters.

[0056] Example 3

[0057] This embodiment relates to an application case of the method for predicting tumor downstaging after neoadjuvant chemoradiotherapy for rectal cancer based on multidimensional habitat imaging features, as described in Embodiment 1.

[0058] 1. Data Preparation (corresponding to step 1): Collect data from a group of pathologically confirmed LARC patients from the hospital information system. All patients underwent standard nCRT and subsequently DWI-MRI and radical surgery. Collect postoperative pathological T-staging to determine TDS status (gold standard). Divide patients into training and independent test sets either chronologically or randomly.

[0059] 2. Image Processing and Feature Extraction: DWI-MRI images were preprocessed using open-source tools such as 3D Slicer and PyRadiomics.

[0060] The tumor ROI was manually delineated on DWI (b=1000 s / mm²) images by two experienced radiologists under blinded conditions, and consensus was reached when there was disagreement.

[0061] In the Python environment, the K-means algorithm from the `scikit-learn` library was used to cluster the grayscale and texture features of voxels within tumor ROIs. The optimal number of clusters, K=3, was determined using the elbow rule and contour coefficient, thus dividing each tumor into three habitat subregions.

[0062] Radiomics features, subregion spatial topology features, and subregion spatial invasiveness features were extracted from the entire tumor ROI and three subregions using the PyRadiomics library.

[0063] Extract clinicopathological features from medical records, including mTRG (0-1 vs 2-3) and postoperative N staging (N0 vs N+).

[0064] 3. Feature selection and model training: The extracted high-dimensional features were standardized and feature selected. First, univariate analysis (t-test / Mann-Whitney U test) was used to screen for features that showed significant differences between the TDS and non-TDS groups (p<0.05). Then, the mRMR algorithm was used for further dimensionality reduction, and finally, about 200 of the most predictive features were retained to form a joint feature set.

[0065] Subsequently, a genetic algorithm was used to screen and optimize the joint feature set and clinical features. A binary chromosome population containing fixed clinical features was initialized, with a population size of 200. For the feature subset represented by each individual, the XGBoost classifier was trained on the training set and its performance was evaluated on the validation set. The fitness value was calculated using the area under the receiver operating characteristic curve as the core indicator, while a feature quantity penalty term was introduced to control the model complexity.

[0066] Then, genetic operations are performed, and dominant individuals are selected through a tournament selection strategy. A new generation of population is generated by two-point crossover and position flip mutation operators. During the crossover and mutation process, it is ensured that the clinical characteristic gene positions are not changed, and an elite preservation mechanism is used to maintain the optimal solution of the population.

[0067] The convergence status of the evolution process is continuously monitored, and optimization is terminated when the preset maximum number of iterations is reached or the fitness does not improve significantly for several consecutive generations.

[0068] The final output is the feature subset with the best overall performance on the validation set, which provides optimized feature input for subsequent model construction. XGBoost is then retrained using the best subset to obtain the optimal model.

[0069] 4. Clinical predictive applications: Integrate the trained optimal model into the prediction system.

[0070] When a new LARC patient completes nCRT, their DWI-MRI images are imported into the system. The system automatically performs preprocessing, segmentation, habitat classification, and feature extraction, and combines the patient's mTRG and N-staging information to generate a joint feature vector.

[0071] The system inputs the feature vector into the optimal model, calculates and outputs the predicted probability of the patient reaching TDS in real time.

[0072] The system generates a report, including a habitat segmentation visualization, predicted probabilities, and clinical recommendations based on preset thresholds (e.g., >0.7), such as "suitable for in-depth evaluation pending observation strategy." Physicians can use this report to discuss individualized treatment plans with patients.

[0073] Those skilled in the art will understand that, without departing from the principles and spirit of this invention, the image segmentation methods (such as those using deep learning for automatic segmentation), clustering algorithms, feature types, machine learning models, etc., in the above embodiments can be replaced or modified. For example, U-Net can be used for automatic tumor segmentation, or a deep convolutional neural network can be used for end-to-end feature learning and classification. All such modifications should be considered to fall within the scope of protection defined by the claims of this invention.

Claims

1. A method for predicting neoadjuvant chemoradiotherapy for tumors based on multidimensional habitat imaging features, characterized in that, Includes the following steps: Step 1. Data Acquisition and Preprocessing: Acquire diffusion-weighted magnetic resonance imaging (DWI-MRI) data and corresponding clinicopathological data of patients with locally advanced rectal cancer after neoadjuvant chemoradiotherapy (nCRT); preprocess the DWI-MRI data. Step 2. Build a prediction model, which includes a tumor segmentation and habitat subregion division module, a multi-dimensional feature extraction module, and a feature screening and training module connected in sequence. Among them, the tumor segmentation and habitat subregion division module, based on the region of interest (ROI) on the preprocessed DWI-MRI image, uses unsupervised clustering to divide the tumor into K visually distinguishable habitat subregions. The multi-dimensional feature extraction module extracts image omics features, sub-region topological features, and sub-region spatial invasiveness features from habitat sub-regions. The feature selection and training module initially screens the features extracted by the multi-dimensional feature extraction module using univariate analysis, then uses the mRMR algorithm to reduce the dimensionality and obtain a joint feature set. Subsequently, a genetic algorithm is used to further screen and optimize the joint feature set and clinical features. First, a chromosome population containing fixed clinical features is initialized. The XGBoost classifier is trained on the training set and its performance is evaluated on the validation set. The fitness is calculated using the area under the receiver operating characteristic curve as the main indicator. Then, genetic operations are performed to update the population through tournament selection, two-point crossover, and position flip mutation, while protecting the fixed clinical features from being changed. Then, the convergence is monitored, and optimization is terminated when the maximum number of iterations or the fitness plateau is reached. Finally, the feature subset with the best overall performance on the validation set is output, and XGBoost is retrained using the best subset to obtain the final model. Step 3. Prediction and Output: Input the joint feature set of the patient to be predicted into the trained prediction model to obtain the predicted probability or classification label of the patient reaching TDS, and output the prediction result. The trained model is then transformed into a stable and interpretable clinical decision support tool to achieve automated prediction and result output for new patients.

2. The method according to step 1, characterized in that, Step 1 includes: Step 101: Data Acquisition Collect pelvic multiparameter magnetic resonance images of patients with pathologically confirmed locally advanced rectal cancer (LARC) after neoadjuvant chemoradiotherapy (nCRT) and before surgery. Step 102: Bias field correction, commonly using the N4 algorithm, which iteratively solves for the bias field through maximum likelihood estimation. B(x)B(x) ; Step 103: Image registration. Align images from different sequences or time points of the same patient spatially to ensure consistency of the analysis area in subsequent steps. Use rigid or affine transformations to register DWI images with high-resolution T2WI images as a fixed reference. Step 104: Resampling. Standardize the image voxel size to eliminate anisotropic resolution differences caused by different scanning protocols. Use linear interpolation or B-spline interpolation to resample all images to isotropic voxels. The interpolation formula is as follows: (1,1) in The grayscale value of the target image. The grayscale value of the original image; Step 105: Intensity standardization, to make the grayscale values ​​of images from different scanning devices and different patients comparable. Z-score standardization is calculated using the following formula: (1,2) in μ and σ It represents the mean and standard deviation of an image or a specific region.

3. The method according to step 2 of claim, characterized in that, Step 2 specifically includes: Step 201: First, a simple linear iterative clustering algorithm is used to segment the tumor region into multiple compact, uniform three-dimensional hypervoxels; then, a K-means clustering algorithm is used to cluster the hypervoxel feature vector sets to discover different image phenotype "habitats". N Each hypervoxon is divided into K Within each cluster, the intra-class squared error is minimized; the number of habitat subregions K is optimized using at least one evaluation index among the silhouette coefficient, Calinski-Harabasz index, and Davies-Bouldin index, and the optimal value is determined by combining the elbow rule with the silhouette coefficient. K value; Step 202: Extract radiomics features from each of the identified habitat subregions; Step 203: To further quantify the spatial configuration and invasive potential of intratumoral heterogeneity, extract the spatial topological features of subregions; Step 204: To focus on the morphological complexity of the tumor-normal tissue boundary or subregion margins to assess their invasiveness, extract the spatial invasiveness features of the subregions; Step 205: Extract clinicopathological features, including at least one of the following: modified tumor regression grade mTRG, postoperative N stage, carcinoembryonic antigen (CEA) level, and carbohydrate antigen (CA19-9) level; Step 206: Standardize and select features from the extracted high-dimensional features; first, use univariate analysis to screen out features that are significantly different between the TDS and non-TDS groups; then use the mRMR algorithm to further reduce the dimensionality, and finally retain the most predictive features to form a joint feature set; Step 207: Subsequently, a genetic algorithm is used to screen and optimize the joint feature set and clinical features. A binary chromosome population containing fixed clinical features is initialized. For the feature subset represented by each individual, the XGBoost classifier is trained on the training set and its performance is evaluated on the validation set. The fitness value is calculated using the area under the receiver operating characteristic curve as the core indicator, and a feature quantity penalty term is introduced to control the model complexity. Then, genetic operations are performed. A tournament selection strategy is used to select superior individuals. A new generation of population is generated using two-point crossover and position flip mutation operators. During the crossover and mutation process, it is ensured that the clinical feature gene positions are not changed. An elite retention mechanism is used to maintain the optimal solution of the population. The convergence status of the evolutionary process is continuously monitored. The optimization is terminated when the preset maximum number of iterations is reached or the fitness has not significantly improved for several consecutive generations. Finally, the feature subset with the best overall performance on the validation set is output, providing optimized feature input for subsequent model construction.

4. The method as described in claim 3, characterized in that, The radiomics features mentioned in step 202 include: Shape characteristics: Describe the three-dimensional geometric morphology of the tumor or its subregions; First-order statistical characteristics: describe the distribution of voxel intensity within a region; Texture features: describe the spatial arrangement and interrelationships of voxel intensities, and are calculated using methods such as gray-level co-occurrence matrix, gray-level run-length matrix, gray-level size region matrix, and neighborhood gray-level difference matrix.

5. The method as described in claim 3, characterized in that, The subregional spatial topological features mentioned in step 203 include: Dominant subregion analysis features: Calculate the volume percentage of each subregion relative to the total tumor volume to identify dominant subregions; Subregion spatial distribution characteristics: By calculating the spatial contact surface area and voxel adjacency relationship between different subregions, or by using spatial correlation analysis, it can be determined whether a specific subregion tends to be adjacent or separated in space; The formula for calculating the spatial relationship of subregions is as follows: (2,1) in and The grayscale value of a pixel.

6. The method as described in claim 3, characterized in that, The subregional spatial invasiveness features described in step 204 include: The degree of edge wrinkling is calculated using the following formula: (2,2) It is the size of pixels. N Minimum number required for coverage The spatial regularity characteristic is calculated using the following formula: (2,3) V For volume, A This represents the surface area; the closer the value is to 1, the more regular the surface; the smaller the value, the more irregular the surface. Interface sharpness features are used to assess the steepness of the gradient change in image signal intensity between a specific subregion and its adjacent normal tissue or other subregions. Blurred interfaces suggest invasive growth. At the boundary between the subregion and adjacent normal tissue, a narrow band is extracted along the normal direction, and the average gradient amplitude of the DWI signal within this narrow band is calculated using the following formula: (2,4) Band is a set of pixels, and p is a single pixel. N The total number of pixels, the lower G The value indicates blurred boundaries, suggesting invasive growth under a microscope.

7. The method according to step 3, characterized in that, Step 3 specifically includes: Step 301: Use the area under the receiver operating characteristic curve (AUC), calibration curve, and decision curve analysis to comprehensively evaluate the model's discrimination, calibration, and clinical usefulness using DCA. Step 302: Save the final trained optimal model in a standardized format to ensure that it can be stably loaded and called in the production environment.

8. A system for implementing the method according to any one of claims 1-7, characterized in that, include: The data acquisition and preprocessing module is used to acquire and preprocess DWI-MRI and clinicopathological data; The tumor segmentation and habitat subregion division module is used to divide and perform tumor segmentation and habitat subregion division. A multi-dimensional feature extraction module is used to extract multi-dimensional radiomics features; The feature selection and training module is used to build, train, and store the TDS prediction model; The predictive inference module is used to load the predictive model, perform TDS prediction on new patient data, and output the results. The user interface is used to input data, trigger analysis processes, and visualize prediction results and habitat subregion maps.

Images