A method for predicting the efficacy of neoadjuvant chemotherapy for breast cancer by DCE-MRI

Abstract

This invention discloses a method for predicting the efficacy of neoadjuvant chemotherapy for breast cancer using DCE-MRI, which relates to the fields of medical image analysis and artificial intelligence. The method includes: cropping and intensity normalizing the region of interest (ROI) of four pre-treatment dynamic contrast-enhanced magnetic resonance imaging (MRI) images; calculating peak intensity, clearance rate, and early enhancement rate from the RIO to form a TIC feature vector, generating modulation weights via a multilayer perceptron, and re-weighting the backbone features to obtain modulation features; dividing the tumor voxels into four habitat regions based on the sample adaptive quantile threshold, extracting the features of each region as graph nodes, and then aggregating and propagating these nodes through a graph attention layer and a graph convolutional layer before mean pooling to obtain spatial habitat features; concatenating the modulation features and spatial habitat features along the feature dimension and outputting the pCR prediction probability through a two-stage classification head. This invention achieves high-precision efficacy prediction in a single pre-treatment scan and exhibits strong cross-institutional generalization ability.

Description

Technical Field

[0001] This invention relates to the fields of medical image analysis and artificial intelligence technology, and in particular to a method for predicting the efficacy of neoadjuvant chemotherapy for breast cancer using DCE-MRI. Background Technology

[0002] Breast cancer is the leading cause of cancer death worldwide. Neoadjuvant chemotherapy (NAC) is a crucial treatment for locally advanced breast cancer, capable of reducing tumor stage, increasing breast-conserving surgery rates, and improving patient prognosis. Pathological complete response (pCR) is the gold standard for evaluating NAC efficacy, defined as the absence of residual invasive lesions in the surgical specimen, and is a strong predictor of long-term survival. However, pCR can only be determined through postoperative pathological examination and cannot be known before treatment. Therefore, accurate pre-treatment pCR prediction has significant clinical value, helping oncologists develop individualized treatment plans and avoid the toxic side effects of ineffective chemotherapy.

[0003] Dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI), through multi-phase sequence acquisition, simultaneously captures tumor morphological and hemodynamic features, and has become the main imaging method for evaluating the efficacy of NAC (negative arterial cancer). The core diagnostic value of DCE-MRI lies in the time-intensity curve (TIC), whose dynamic characteristics include peak enhancement intensity (PI), washout rate (WR), and early enhancement rate (EER), reflecting the tumor's neovascularization density, permeability, and perfusion status.

[0004] The main drawbacks of the existing technology are as follows: Limitations of single-phase DCE-MRI methods: Most existing deep learning methods utilize only single-phase DCE-MRI phase images, ignoring the rich pharmacokinetic information encoded by multi-phase sequences. Although some methods attempt to stitch multi-phase images along the channel dimension (early fusion strategy), this approach is merely a simple juxtaposition along the channel direction and cannot explicitly model the temporal kinetics of contrast agent uptake, resulting in low efficiency in utilizing TIC dynamic information and a significant increase in the number of parameters.

[0005] Insufficient temporal feature mining: Some methods use TIC parameters as independent classification features, failing to integrate pharmacokinetic information into deep feature representation. How to effectively integrate TIC dynamic information into deep feature learning while preserving the spatial heterogeneity of voxel-level TIC parameters remains a key challenge.

[0006] Insufficient modeling of spatial heterogeneity within tumors: Most existing methods treat the tumor as a whole, neglecting the spatial enhancement heterogeneity of different functional regions within the tumor. Although radiomics methods based on habitat imaging attempt to segment tumor subregions, these methods rely on manually extracted radiomics features and do not model the spatial interactions and biological interdependencies between adjacent habitat regions.

[0007] Insufficient utilization of complementary spatiotemporal information: Tumor response to NAC is closely related to its spatiotemporal heterogeneity, reflected in the temporal changes of TIC enhancement kinetics in multi-phase DCE-MRI sequences and the spatial distribution of different enhancement behaviors within the tumor. Existing methods have not fully utilized information from these two complementary dimensions simultaneously, resulting in insufficient mining of key pharmacokinetic information. Summary of the Invention

[0008] The purpose of this invention is to provide a method for predicting the efficacy of neoadjuvant chemotherapy for breast cancer using DCE-MRI, thereby solving the problem that the existing technology cannot effectively integrate the temporal pharmacokinetic information of multi-phase DCE-MRI with the spatial heterogeneity information of the enhancement pattern inside the tumor, thus failing to achieve accurate pre-treatment prediction of complete pathological remission of breast cancer with neoadjuvant chemotherapy.

[0009] To achieve the above objectives, this invention provides a DCE-MRI method for predicting the efficacy of neoadjuvant chemotherapy in breast cancer, comprising the following steps: Step 1: Obtain the four-phase dynamic contrast-enhanced magnetic resonance imaging image sequence of breast cancer samples before treatment, and record them as the baseline image, the first phase after enhancement, the second phase after enhancement, and the third phase after enhancement. The four-phase dynamic contrast-enhanced magnetic resonance imaging images are cropped and their intensity is normalized to obtain the normalized four-phase voxel intensity field. Step 2: Using a 3D residual convolutional network as the backbone network, extract backbone feature vectors from the normalized second-stage enhanced image; calculate the overall peak intensity, overall clearance rate, and overall early enhancement rate based on the phase-by-phase average voxel intensity of the region of interest within the cropped four-stage dynamic contrast-enhanced magnetic resonance imaging image, and stack them into TIC feature vectors; feed the TIC feature vectors into a multilayer perceptron and generate modulation weights through the Sigmoid activation function; multiply the modulation weights by the backbone feature vectors to obtain the modulation features; Step 3: Calculate the voxel-by-voxel early enhancement rate and voxel-by-voxel clearing rate for each voxel in the region of interest; calculate the adaptive quantile threshold for the sample based on the distribution of the voxel-by-voxel early enhancement rate and voxel-by-voxel clearing rate of all voxels in the current sample; divide the voxels in the region of interest into four habitat regions according to the threshold, and generate a corresponding mask volume for each habitat region; extract the region feature vector from each mask volume through four independent 3D residual convolution branches, and stack the four region feature vectors to form a graph node feature matrix; calculate the adaptive adjacency matrix from the graph node feature matrix through a multi-head graph attention layer, then propagate and aggregate the graph node feature matrix through two graph convolution layers in sequence, and finally perform mean pooling on the output of the last graph convolution layer along the node dimension to obtain the spatial habitat features; Step 4: Concatenate the modulation features with the spatial habitat features to obtain a fused feature vector; input the fused feature vector into the two-stage classification head and output the pCR prediction probability.

[0010] Therefore, the present invention employs the above-mentioned DCE-MRI method for predicting the efficacy of neoadjuvant chemotherapy in breast cancer, which has the following beneficial effects: 1. By using a pharmacokinetic feature modulation module to generate modulation weights with a lightweight multilayer perceptron to reweight the backbone features through channels, the expansion of backbone network parameters caused by direct splicing of multi-phase images along channels is avoided. Global time-intensity curve dynamics information is explicitly injected with only a small number of parameters, which can stably improve the accuracy of efficacy prediction compared with single-phase baseline methods. 2. By using a habitat graph convolutional network based on time-intensity curves to divide tumor voxels into four habitat regions with sample adaptive quantile thresholds, and using graph attention mechanism to model the interaction between different regions, while preserving the spatial heterogeneity of voxel-level time-intensity curve parameters, we can provide an interpretable characterization of tumor pharmacokinetic heterogeneity and locate the habitat region that contributes the most to the therapeutic effect. 3. By splicing modulation features and spatial habitat features along the feature dimension to form a unified fusion representation, and by utilizing information from two complementary dimensions—global time-intensity curve temporal dynamics and tumor intratumoral spatial enhancement heterogeneity—the prediction accuracy can be further improved compared to methods that rely on only a single dimension.

[0011] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0012] Figure 1 This is a flowchart of a DCE-MRI method for predicting the efficacy of neoadjuvant chemotherapy in breast cancer according to the present invention; Figure 2 This is a TIC curve diagram of the first sample of the present invention; Figure 3This is a TIC curve diagram of the second sample of the present invention; Figure 4 This is the TIC curve of the third sample of the present invention; Figure 5 This is the TIC curve of the fourth sample of the present invention; Figure 6 This is a scatter plot showing the Spearman correlation between the modulation weight α and PI in this invention. Figure 7 This is a Spearman correlation scatter plot of the modulation weight α and WR in this invention; Figure 8 This is a Spearman correlation scatter plot of the modulation weight α and EER in this invention; Figure 9 This is a t-SNE visualization of the delineated region of the present invention; Figure 10 This is a visualization of the t-SNE of the continuously enhanced region in this invention; Figure 11 This is a t-SNE visualization of the platform area of ​​this invention; Figure 12 This is a visualization of the t-SNE in the low-intensity region of the present invention; Figure 13 This is a habitat area importance analysis diagram for the present invention. Detailed Implementation

[0013] The following detailed description of embodiments of the invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely illustrates selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the invention without inventive effort are within the scope of protection of the invention.

[0014] Please see Figure 1 A DCE-MRI method for predicting the efficacy of neoadjuvant chemotherapy in breast cancer includes the following steps: Step 1: Acquisition and preprocessing of multi-phase DCE-MRI data.

[0015] Specifically, step 1 includes the following sub-steps: Step 101: Obtain DCE-MRI image sequences of the breast cancer sample in the four pre-treatment phases, and record them as baseline images. Image after first-stage enhancement Second-stage enhanced image The third phase enhanced image ,in Images were collected before intravenous injection of contrast agent. , , These are three phases of enhanced images acquired at set time intervals after contrast agent injection. Each image is cropped to a region of interest centered on the tumor, and the voxel size is uniform. Voxel.

[0016] Step 102 involves intensity normalization for each region of interest, including the following sub-steps: (1) Calculate the 1st percentile of the intensity of all voxels in the current period region of interest. With the 99th percentile Cut off voxel intensity to the interval Internally, to suppress extreme noise; (2) The mean of the truncated image with standard deviation conduct Normalization yields the normalized four-phase voxel intensity field, calculated using the following formula: ; in, Indicates the stage number of the DCE-MRI; Represents the three-dimensional coordinates of a voxel; Indicates the first The coordinates of the region of interest during the period are The strength after voxel truncation; , They represent the first Mean and standard deviation of the truncated image; Indicates the normalized i-th Period in coordinates The voxel intensity at each coordinate. The voxel intensity of all coordinates in each period constitutes the voxel intensity field for that period.

[0017] This normalization strategy can eliminate signal intensity differences between different acquisition devices, different scanning protocols, and different samples, thereby improving the robustness of subsequent feature extraction.

[0018] Step 2: The Pharmacokinetic Feature Modulation Module (PFMM) extracts modulation features from the region of interest.

[0019] The Pharmacokinetic Feature Modulation (PFMM) module, with backbone feature extraction and TIC global modulation as its core, includes the following sub-steps in step 2: Step 201, Backbone Feature Extraction: A pre-trained R3D-18 3D residual convolutional network is used as the backbone network for spatial feature extraction. The first 3D convolutional layer of the backbone network is replaced with a single-channel equivalent layer. The convolutional kernel weights of the single-channel equivalent layer are initialized by averaging the weights of the pre-trained three-channel RGB convolutional kernels along the channel dimension to preserve the existing spatial filter patterns of the pre-trained model. The normalized second-stage enhanced image is input into the R3D-18 backbone network, and a 512-dimensional backbone feature vector is output after global average pooling. : Step 202, TIC Feature Vector Extraction: Based on the cropped four-phase dynamic contrast-enhanced magnetic resonance imaging images, calculate the average voxel intensity within the region of interest in each phase. Based on this, three semi-quantitative TIC scalar parameters are then calculated to form a 3D TIC feature vector. The specific calculations are as follows: (1) The average voxel intensity over time is calculated using the following formula: ; in, Represents the set of coordinates of all voxels within the region of interest; Represents a set The total number of median voxels; Indicates the first The average voxel intensity within the region of interest during the period.

[0020] (2) Overall peak intensity: Reflects the maximum contrast agent uptake capacity of tumor tissue, and the calculation formula is as follows: ; in, Indicates the overall peak intensity. This indicates taking the maximum value.

[0021] (3) Overall clearance rate This reflects the clearance rate of contrast agent within tumor tissue. A positive value indicates significant contrast agent clearance. It serves as a characteristic marker of TIC in malignant and chemotherapy-sensitive lesions. The calculation formula is as follows: ; in, This means truncating the values ​​within the parentheses to a specific range. ; This indicates the overall clearance rate.

[0022] (4) Overall early enhancement rate The initial contrast agent uptake rate from baseline to stage one is quantified, reflecting the early response of tumor tissue to the contrast agent. The calculation formula is as follows: ; in, This means truncating the values ​​within the parentheses to a specific range. ; This indicates the overall early enhancement rate.

[0023] Stacking the three TIC scalar parameters above yields a 3D TIC (Time-Intensity Curve) eigenvector, calculated as follows: ; in, Represents the TIC feature vector. This indicates transpose.

[0024] Step 203, Feature Modulation: Input the TIC feature vector into a multilayer perceptron consisting of three fully connected layers, a batch normalization layer, and a Dropout layer cascaded together. The network structure of the multilayer perceptron is as follows: The output value of the last layer after the Sigmoid activation function ranges from 1 to 10. The modulation weight is calculated using the following formula: ; in, This represents a nonlinear mapping for a three-layer multilayer perceptron; This represents the Sigmoid activation function; The modulation weights are data-driven estimates of the overall TIC dynamics of the tumor.

[0025] Modulation weights With backbone feature vector Channel reweighting is performed to obtain modulation features, calculated using the following formula: ; in, This represents the 512-dimensional modulation feature output by the pharmacokinetic feature modulation module.

[0026] In this way, global TIC enhancement dynamics are injected into spatial feature representation in a channel reweighting manner without modifying the backbone network architecture.

[0027] Step 3: Extract spatial habitat features from the voxel intensity field using the Habitat-GCN habitat graph convolutional network based on TIC.

[0028] The Habitat-GCN habitat graph convolutional network based on TIC has voxel-level TIC dynamic partitioning and graph convolution modeling as its core. Step 3 specifically includes the following sub-steps: Step 301, Voxel-by-Voxel TIC Dynamics Calculation: For each voxel in the region of interest... The early enhancement rate and clearance rate of each voxel were calculated separately, using the following formulas: ; ; in, This represents a numerical stability constant set to prevent the denominator from being zero. Representing coordinates Early enhancement rate of the body at the site, Representing coordinates The clearance rate of each element at the location.

[0029] Step 302, Calculation of adaptive quantile threshold for samples: Based on the voxel-by-voxel TIC feature distribution of all voxels in the region of interest of the current sample, calculate the adaptive quantile threshold for the three samples. The calculation formula is as follows: ; ; ; in, This indicates that the specified feature of all voxels within the region of interest of the current sample is selected. Quantiles , These represent the 75th and 25th quantiles of the early enhancement rate of each voxel, respectively. This represents the 50th quantile of voxel clearance. and These represent the voxel-by-voxel early enhancement rate and voxel-by-voxel clearance rate distributions for all voxels within the region of interest of the current sample, respectively. The sample adaptive quantile threshold can be divided based on the enhancement dynamics distribution of each sample, avoiding the use of a fixed population-level threshold.

[0030] Step 303, Habitat Region Delineation: Based on the voxel-by-voxel TIC features and the sample adaptive quantile threshold, each voxel within the region of interest is assigned to one of the four habitat regions, and a corresponding binary mask and mask volume are generated. The assignment rules are as follows: (1) Clearing area High early enhancement rate and high clearance rate correspond to an enhancement pattern of rapid contrast agent clearance after rapid inflow. ; (2) Continuous reinforcement zone High early enhancement rate and low clearance rate correspond to an enhancement pattern of continuous contrast agent retention. ; (3) Platform area Medium-intensity early enhancement rate, corresponding to moderately enhanced tissue: ; (4) Low-intensity zone Low early enhancement rate corresponds to poor perfusion or necrotic areas. ; in, Indicates the first The set of voxel coordinates contained in a habitat region. , Represents the logical AND operation.

[0031] For each habitat area Generate a binary mask of the same size as the region of interest. (that is, in) The value inside is 1. (The value of the outer region is 0), and then the region of interest of the normalized second-phase voxel intensity field is multiplied element-wise with the binary mask of each habitat region to obtain the first... The mask volume of a habitat region is calculated using the following formula: ; in, Indicates the first Binary mask for each habitat region This indicates element-wise multiplication. Indicates the first Habitat area coordinates The mask volume at that location, finally obtained Indicates the first The mask volume of a habitat region.

[0032] If a certain habitat area If the current sample is an empty set (i.e., the sample has no voxels corresponding to the enhancement mode), then the mask volume of the region is uniformly filled with the global mean of the region of interest of the second-phase voxel intensity field after normalization, as a backup mask volume to ensure the consistency of the input of the subsequent four branches and the stability of training.

[0033] Step 304, Region Feature Extraction: Construct four independent R3D-18 3D residual convolution branches. These four branches have the same network architecture as the backbone network in Step 201, but their parameters are trained independently. The mask volumes of the four habitat regions are then extracted. , , , The four branches are input one-to-one, and the regional feature vectors of each habitat region are output through global average pooling. The four regional feature vectors are stacked row by row to form a graph node feature matrix. The calculation formula is as follows: ; in, Indicates the first Regional feature vectors of each habitat region ; This represents a graph node feature matrix consisting of four habitat region features, the first of which is... OK That is, the first Features of each graph node.

[0034] Step 305, Graph Attention Calculation of Adaptive Adjacency Matrix: The four habitat regions are treated as four nodes in a fully connected graph, and the multi-head graph attention layer calculates the adjacency matrix from the graph node feature matrix. The adaptive adjacency matrix is ​​calculated using the following formula: ; ; in, , Let represent the learnable weight matrices for query projection and key projection, respectively. , The query matrix and key matrix represent the node features, respectively. The scaling factor represents the feature dimension of the query / key projection. Indicates matrix transpose. express function, This represents an adaptive adjacency matrix, whose elements... Quantifying habitat areas and The strength of feature interactions between them.

[0035] Step 306, Graph Convolutional Layer Information Propagation: The graph node feature matrix is ​​propagated and aggregated sequentially using two graph convolutional layers. The calculation formula is as follows: ; in, , The input feature matrix of the graph convolutional layer is... and They represent the first and The output feature matrix of the layer graph convolutional layer, Indicates the first Learnable weight matrix of layered graph convolutional layers; This indicates a modified linear unit activation function; Indicates the first The output feature matrix of the layer graph convolutional layer.

[0036] Step 307, Node Feature Mean Pooling: This involves applying mean pooling to the output of the second graph convolutional layer. Mean pooling is performed along the node dimension to obtain spatial habitat features, calculated using the following formula: ; in, The second layer of the graph convolutional layer represents the first feature matrix output by the second layer. OK, This represents the spatial habitat characteristics output by Habitat-GCN.

[0037] Step 4, Feature fusion and pCR classification.

[0038] Specifically, step 4 includes the following sub-steps: Step 401, Feature Concatenation and Fusion: The modulation features and spatial habitat features are concatenated along the feature dimension to obtain the fused feature vector. The calculation formula is as follows: ; in, Represents the fused feature vector. This represents a concatenation operation along the feature dimension.

[0039] The splicing-based fusion approach preserves the independent contribution of each feature stream, enabling subsequent classification heads to simultaneously leverage global TIC to enhance dynamics and tumor-internal spatial heterogeneity.

[0040] Step 402, Two-stage classification head mapping: fusing feature vectors Input a two-stage classification head to obtain the predicted probability of complete pathological remission. The calculation formula is as follows: ; in, This represents the learnable weight matrix of the first-stage fully connected layer; This represents the learnable weight matrix of the second-stage fully connected layer; This indicates the first Dropout layer. This represents the second Dropout layer (in this embodiment) Pre-dropout rate (It operates before the first-stage fully connected layer). This indicates a batch normalization layer, and a dropout rate is set between the first-stage fully connected layer and the second-stage fully connected layer. Dropout layer (in this embodiment) Pre-dropout rate ) express function; This represents a binary classification prediction probability vector, where the first and second dimensions correspond to the prediction probabilities of non-pCR and pCR, respectively.

[0041] It should be noted that the model training process is as follows: We employ a weighted random sampler and Focal Loss to address the imbalance between positive and negative samples. We use the AdamW optimizer (learning rate 1e-4, weight decay 1e-4) and a cosine annealing learning rate scheduling strategy. During training, we utilize data augmentation methods such as random flipping, random rotation, brightness / contrast jitter, and Gaussian noise. On the internal dataset, we use sample-level hierarchical five-fold cross-validation to determine the optimal hyperparameter configuration, with the random seed fixed at 42 to ensure repeatability.

[0042] To further illustrate the implementation process and beneficial effects of the present invention, a multicenter DCE-MRI breast cancer dataset will be used as an example.

[0043] First, data preparation: Taking the multicenter DCE-MRI breast cancer dataset as an example, the internal dataset comprises three cohorts: the I-SPY2 public dataset (389 cases, pCR rate 34.96%), the Duke University breast cancer MRI dataset (281 cases, pCR rate 22.06%), and the Zunyi First People's Hospital (FPHZ, 72 cases, pCR rate 18.06%), totaling 742 cases (pCR: 211 cases, 28.44%; non-pCR: 531 cases, 71.56%). Stratified randomization at the sample level (80 / 20) resulted in a training set of 594 cases and an internal test set of 148 cases. The external test cohorts included Guizhou Provincial People's Hospital (GPPH, 264 cases) and Sun Yat-sen Memorial Hospital of Sun Yat-sen University (SYMH, 256 cases). Informed consent was obtained from all participants before data collection.

[0044] All samples underwent DCE-MRI scanning, and four-phase sequentially enhanced sequences were extracted. Tumor ROIs were manually segmented by radiologists, and all volumes were cropped and resampled. Voxel.

[0045] Second, data preprocessing: Intensity normalization is performed on each image: the 1st and 99th percentiles are calculated, the signal intensity is truncated to this range, and then the mean is subtracted and divided by the standard deviation. This normalization strategy effectively eliminates signal intensity differences across institutions and devices, improving the model's generalization ability.

[0046] Third, the PFMM module implementation: Backbone feature extraction: Using R3D-18 pre-trained on Kinetics-400, the first convolutional layer is replaced with a single-channel equivalent layer (the weights are initialized by averaging the original RGB weights along the channel dimension). The normalized second-phase voxel intensity field is input, and a 512-dimensional backbone feature vector is output after global average pooling.

[0047] TIC Feature Calculation: Based on the Mean Intensity of ROI in Four Periods ( The overall peak intensity, overall clearance rate, and overall early enhancement rate are calculated and stacked to form a 3D TIC feature vector.

[0048] Feature modulation: Three-layer MLP ( With BN and Dropout Mapping TIC feature vectors to modulation weights The modulation features are obtained by multiplying the modulation weights by the backbone feature vector.

[0049] Fourth, the Habitat-GCN module implementation: Voxel-by-voxel TIC feature calculation: Calculate EER and WR for each voxel, using the sample adaptive quantile threshold (EER). / , WR's The voxels are divided into four functional areas: the clearance zone, the sustained reinforcement zone, the plateau zone, and the low reinforcement zone.

[0050] Region feature extraction: Four independent R3D-18 branches (with the same architecture as the backbone network but with independent parameters) process the masked volume of each habitat region. Each branch outputs a 512-dimensional region feature vector. The four region feature vectors are stacked to form a graph node feature matrix. .

[0051] Graph Attention and Graph Convolution: Multi-head Graph Attention Layer (4 attention heads, 128 projection dimensions) computes the adaptive adjacency matrix. Two graph convolutional layers propagate information sequentially, and the spatial habitat features are obtained by mean pooling of the features from the four nodes. Fifth, feature fusion and training: The modulation features and spatial habitat features are concatenated to obtain the fused feature vector. After two-stage classification (Dropout(0.5)-Linear( )-BN-ReLU-Dropout(0.25)-Linear( The pCR prediction probability is output by Softmax.

[0052] Training configuration: AdamW optimizer (learning rate = Weight decay = Cosine annealing learning rate scheduling ( ), period length Let FocalLoss be the total number of training iterations. =reciprocal of the number of samples in class pCR WeightedRandomSampler handles class imbalance and gradient clipping. The gradient accumulation steps are 4. The batch size is 4, and the training lasts for 75 epochs. The random seed is fixed at 42.

[0053] Five-fold cross-validation determines the optimal hyperparameters on the internal dataset, and finally retrains on the entire training set to obtain the final model.

[0054] Sixth, performance evaluation: The pCR prediction performance of the method of this invention was verified on different test sets. Performance was evaluated using four metrics: AUC (Area Under the Operating Characteristic), ACC (Accuracy), PPV (Positive Predictive Value), and NPV (Negative Predictive Value). On the internal test set (148 cases), the model (PHGN) constructed using the method of this invention achieved a prediction performance of AUC=0.890, ACC=0.884, PPV=0.844, and NPV=0.933. Without any fine-tuning, it achieved an AUC=0.856 (ACC=0.856, PPV=0.810, NPV=0.860) on the GPPH external cohort and an AUC=0.850 (ACC=0.881, PPV=0.846, NPV=0.888) on the SYMH external cohort, demonstrating robust cross-institutional generalization ability. (See Table 1). Table 1. pCR prediction performance of the method of the present invention on different test sets.

[0055] Seventh, comparative verification: Under the same backbone network and training configuration, compared with the existing single-phase DCE-MRI baseline method, the introduction of PFMM improved the predicted AUC by approximately 0.116; compared with the early fusion strategy of directly stitching four-phase images along the channel, the introduction of PFMM improved the AUC by approximately 0.101 with only a very small increase in parameters. Habitat importance ablation analysis further showed that the clearance zone contributed the most to the predictive performance; ablation of this region alone reduced the AUC by approximately 0.023, consistent with its biological significance of high vascular permeability and high perfusion. These results validate the advantages of the method of this invention in terms of predictive accuracy, parameter efficiency, and biological interpretability.

[0056] Figure 2-5The TIC curves and their corresponding α values ​​for four samples are shown. The four samples are represented as the first sample (non-pCR, α=0.4735), the second sample (pCR, α=0.4732), the third sample (non-pCR, α=0.1509), and the fourth sample (pCR, α=0.2029). This shows that α differs significantly among different samples, and both pCR and non-pCR cases are distributed across the entire range of α, confirming that α acts as an adaptive scaling factor for samples rather than a population-level discriminative feature.

[0057] Figure 6-8 This is a Spearman scatter plot showing the correlation between the modulation weight α and the three TIC parameters (PI, WR, EER). The results show that α is strongly negatively correlated with PI (rs=-0.782, p<0.001), moderately negatively correlated with EER (rs=-0.316, p<0.001), and not significantly correlated with WR (rs=-0.015, p=0.796). This indicates that PI has the most dominant influence on α, and the sample adaptive recalibration effectively reduces the signal intensity variability between samples. The lack of correlation between WR and α indicates that the mean rinsing dynamics at the global ROI level has lost spatial distribution information, providing a mechanistic basis for the complementary design of PFMM and Habitat-GCN.

[0058] Figure 9-12 The t-SNE visualization of the four habitat regions in this invention displays the two-dimensional projection of the 512-dimensional feature vectors of each region. The clearance region shows the clearest separation between the pCR and non-pCR clusters, the separation between groups in the sustained enhancement region and the low enhancement region is comparable, and the plateau region has the greatest inter-group overlap, indicating that there are significant differences in the discriminative contribution of different habitat regions to pCR prediction, and the vascular microenvironment features of the clearance region have the greatest diagnostic value.

[0059] Figure 13 This is a habitat region importance analysis chart for this invention. Each bar represents the AUC decrease after ablation of the corresponding habitat region. Ablation of the clearance zone resulted in the largest AUC decrease (ΔAUC=0.023), the sustained enhancement zone and the low enhancement zone contributed equally (ΔAUC=0.012), and the plateau zone had the least impact (ΔAUC=0.004). The quantitative ranking is completely consistent with the t-SNE visualization results, mutually confirming the differentiated contribution of each habitat region to pCR prediction, indicating that the spatial distribution of clearance-type enhancement dynamics is the most discriminative imaging signal for predicting complete pathological remission.

[0060] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit them. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the technical solutions of the present invention, and these modifications or equivalent substitutions cannot cause the modified technical solutions to deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1. A DCE-MRI method for predicting the efficacy of neoadjuvant chemotherapy in breast cancer, characterized in that, Includes the following steps: Step 1: Obtain the four-phase dynamic contrast-enhanced magnetic resonance imaging image sequence of breast cancer samples before treatment, and record them as the baseline image, the first phase after enhancement, the second phase after enhancement, and the third phase after enhancement. The four-phase dynamic contrast-enhanced magnetic resonance imaging images are cropped and their intensity is normalized to obtain the normalized four-phase voxel intensity field. Step 2: Using a 3D residual convolutional network as the backbone network, extract backbone feature vectors from the normalized second-stage enhanced image; calculate the overall peak intensity, overall clearance rate, and overall early enhancement rate based on the phase-by-phase average voxel intensity of the region of interest in the cropped four-stage dynamic contrast-enhanced magnetic resonance imaging image, and stack them into TIC feature vectors; feed the TIC feature vectors into a multilayer perceptron and generate modulation weights through the Sigmoid activation function; The modulation weights are multiplied by the backbone feature vector to obtain the modulation features; Step 3: Calculate the voxel-by-voxel early enhancement rate and voxel-by-voxel clearing rate for each voxel in the region of interest; calculate the sample adaptive quantile threshold based on the distribution of the voxel-by-voxel early enhancement rate and voxel-by-voxel clearing rate of all voxels in the current sample; divide the voxels in the region of interest into four habitat regions according to the threshold, and generate a corresponding mask volume for each habitat region; extract the region feature vector from each mask volume through four independent 3D residual convolution branches, and stack the four region feature vectors to form a graph node feature matrix; The adaptive adjacency matrix is ​​calculated from the graph node feature matrix through a multi-head graph attention layer. Then, the graph node feature matrix is ​​propagated and aggregated through two graph convolutional layers. Finally, the output of the last graph convolutional layer is averaged along the node dimension to obtain the spatial habitat features. Step 4: Concatenate the modulation features with the spatial habitat features to obtain a fused feature vector; input the fused feature vector into the two-stage classification head and output the pCR prediction probability.

2. The method for predicting the efficacy of neoadjuvant chemotherapy for breast cancer using DCE-MRI according to claim 1, characterized in that, In step 1, the four phases of dynamic contrast-enhanced magnetic resonance imaging images are cropped and their intensity normalized to obtain the normalized voxel intensity fields for the four phases, including: The four phases of dynamic contrast-enhanced magnetic resonance imaging images were cropped to regions of interest centered on the tumor and uniformly sized to obtain the regions of interest for each phase. For each region of interest, first calculate the 1st percentile of the intensity of all voxels within it. With the 99th percentile Cut off voxel intensity to the interval Then, Z-score normalization is performed using the mean and standard deviation of the truncated image to obtain the normalized voxel intensity field for each period. ,in Indicates the issue number. Represents the three-dimensional coordinates of a voxel.

3. The method for predicting the efficacy of neoadjuvant chemotherapy for breast cancer using DCE-MRI according to claim 2, characterized in that, In step 2, using a three-dimensional residual convolutional network as the backbone network, the extraction of backbone feature vectors from the normalized second-stage enhanced image includes: Using a pre-trained R3D-18 3D residual convolutional network as the backbone network, the first 3D convolutional layer of the backbone network is replaced with a single-channel equivalent layer. The convolutional kernel weights of the single-channel equivalent layer are initialized by averaging the pre-trained three-channel RGB convolutional kernel weights along the channel dimension. The normalized second-stage enhanced image is input into the backbone network, and the backbone feature vector is output after 3D global average pooling.

4. The method for predicting the efficacy of neoadjuvant chemotherapy for breast cancer using DCE-MRI according to claim 3, characterized in that, In step 2, the overall peak intensity, overall clearance rate, and overall early enhancement rate are calculated based on the phase-by-phase average voxel intensity of the region of interest within the cropped four-phase dynamic contrast-enhanced magnetic resonance imaging image, and stacked into a TIC feature vector, including: Based on the cropped four-phase dynamic contrast-enhanced magnetic resonance imaging images, the average voxel intensity of all voxels in the region of interest in each phase is calculated to obtain the phase-by-phase average voxel intensity. The overall peak intensity, overall clearance rate, and overall early enhancement rate are calculated based on the phase-by-phase average voxel intensity, using the following formulas: ； ； ； in, This represents the period-by-period average voxel intensity in period k. , Indicates the overall peak intensity. Indicates the overall clearance rate. Indicates the overall early enhancement rate. This means truncating the values ​​within the parentheses to a specific range. , This means truncating the values ​​within the parentheses to a specific range. , This indicates taking the maximum value; The overall peak intensity, overall clearance rate, and overall early enhancement rate are stacked into a 3D TIC feature vector. , , This indicates transpose.

5. The method for predicting the efficacy of neoadjuvant chemotherapy for breast cancer using DCE-MRI according to claim 4, characterized in that, In step 2, the network structure of the multilayer perceptron consists of three fully connected layers, a batch normalization layer, and a Dropout layer cascaded together.

6. The method for predicting the efficacy of neoadjuvant chemotherapy for breast cancer using DCE-MRI according to claim 5, characterized in that, In step 3, the voxel-by-voxel early enhancement rate and voxel-by-voxel clearance rate are calculated for each voxel in the region of interest. The adaptive quantile threshold for the sample is calculated based on the distribution of the voxel-by-voxel early enhancement rate and voxel-by-voxel clearance rate of all voxels in the current sample, including: For each voxel within the region of interest, calculate the voxel-by-voxel early enhancement rate and voxel-by-voxel clearance rate using the following formulas: ； ； in, Voxel representation The early enhancement rate of the body's components, Voxel representation The clearance rate of each individual element. Represents the numerical stability constant; Early enhancement rate per voxel based on all voxels within the region of interest of the current sample. Distribution of individual element clearance The adaptive quantile threshold for the sample is calculated using the following formula: ； ； ； in, This indicates that the specified feature of all voxels within the region of interest of the current sample is selected. percentile, , These represent the 75th and 25th percentiles of the early enhancement rate of each voxel, respectively. This represents the 50th percentile of the voxel clearance rate.

7. The method for predicting the efficacy of neoadjuvant chemotherapy for breast cancer using DCE-MRI according to claim 6, characterized in that, In step 3, the voxels within the region of interest are divided into four habitat regions according to a threshold, and a corresponding mask volume is generated for each habitat region, including: Will simultaneously satisfy and Voxels are allocated to the clearance zone; Will simultaneously satisfy and Voxels are allocated to the sustained reinforcement zone; Will satisfy The voxels are distributed to the plateau region; Will satisfy Voxels are allocated to low-enhancement areas; Generate a binary mask of the same size as the region of interest for each habitat region. Multiply the region of interest of the normalized second-phase voxel intensity field with the binary mask of each habitat region element by element to obtain the mask volume of each habitat region.

8. The method for predicting the efficacy of neoadjuvant chemotherapy for breast cancer using DCE-MRI according to claim 7, characterized in that, In step 3, if a habitat region is an empty set in the current sample, the mask volume of the habitat region is uniformly filled with the global mean of the region of interest of the normalized second-phase voxel intensity field, as a backup mask volume.

9. The method for predicting the efficacy of neoadjuvant chemotherapy for breast cancer using DCE-MRI according to claim 8, characterized in that, In step 3, an adaptive adjacency matrix is ​​calculated from the graph node feature matrix through a multi-head graph attention layer. Then, the graph node feature matrix is ​​propagated and aggregated sequentially through two graph convolutional layers. Finally, the output of the last graph convolutional layer is average-pooled along the node dimension to obtain the spatial habitat features, including: The four habitat regions are treated as four nodes in a fully connected graph. The feature matrices of the graph nodes are used to obtain the query matrix and key matrix through query projection weights and key projection weights, respectively. Then, the adaptive adjacency matrix is ​​calculated using the following formula: ； ； in, Represents the feature matrix of graph nodes. and These represent query projection weight and key projection weight, respectively. and These represent the query matrix and the key matrix, respectively. Indicates the scaling factor. express The function, where A represents the adaptive adjacency matrix; The graph node feature matrix is ​​propagated and aggregated sequentially using two graph convolutional layers: ； in , and They represent the first and The output feature matrix of the layer graph convolutional layer, For the first The learnable weight matrix of a layered graph convolutional layer. This indicates a modified linear unit activation function; The output of the second graph convolutional layer Mean pooling is performed along the node dimension to obtain spatial habitat features, calculated using the following formula: ； in, The second layer of the graph convolutional layer represents the first feature matrix output by the second layer. OK, It indicates the characteristics of the spatial habitat.

10. A method for predicting the efficacy of neoadjuvant chemotherapy for breast cancer using DCE-MRI according to claim 9, characterized in that, In step 4, the fused feature vector is input into the two-stage classification head, and the output pCR prediction probability includes: The two-stage classification head consists of a first Dropout layer, a first-stage fully connected layer, a modified linear unit activation function, a batch normalization layer, a second Dropout layer, a second-stage fully connected layer, and a Softmax function, all cascaded in sequence. The pCR prediction probability is output from the two-stage classification head using the following formula: ； in, Represents the fused feature vector. This indicates the first Dropout layer. This indicates the second Dropout layer. and Let represent the learnable weight matrices of the first-stage fully connected layer and the second-stage fully connected layer, respectively. For batch normalization layer, This represents the pCR prediction probability.

Images