Accurately predicting pathological complete response (pCR) to neoadjuvant treatment (NAT) in breast cancer remains challenging due to tumor heterogeneity. This study enrolled 2279 patients across 12 centers and develops a novel multi-modality model integrating longitudinal magnetic resonance imaging (MRI) spatial habitat radiomics, transcriptomics, and single-cell RNA sequencing for predicting pCR. By analyzing tumor subregions on multi-timepoint MRI, the model captures dynamic intra-tumoral heterogeneity during NAT. It shows superior performance over traditional radiomics, with areas under the curve of 0.863, 0.813, and 0.888 in the external validation, immunotherapy, and multi-omics cohorts, respectively. Subgroup analysis shows its robustness across varying molecular subtypes and clinical stages. Transcriptomic and single-cell RNA sequencing analysis reveals that high model scores correlate with increased immune activity, notably elevated B cell infiltration, indicating the biological basis of the imaging model. The integration of imaging and molecular data demonstrates promise in spatial habitat radiomics to monitor dynamic changes in tumor heterogeneity during NAT. In clinical practice, this study provides a noninvasive tool to accurately predict pCR, with the potential to guide treatment planning and improve breast-conserving surgery rates. Despite promising results, the model requires prospective validation to confirm its utility across diverse patient populations and clinical settings.