<p>Intratumor phenotypic heterogeneity plays a pivotal role in shaping the evolutionary dynamics of tumor growth and drug responses. Nonlocal reaction-diffusion equations have been developed to model the spatio-temporal dynamics of tumor cells with heterogeneous phenotypes. However, rigorous mathematical analysis of such nonlocal reaction-diffusion equations with high spatial and phenotypic dimensions remains challenging. Additionally, previous studies have seldom explored key determinants of phenotypic heterogeneity by integrating models with realistic spatial data. To address these challenges and gaps, we present a rigorous mathematical framework for dissecting phenotypic heterogeneity and evolution within vascularized tumors, utilizing a nonlocal reaction-diffusion model coupled with spatial gene expression data. First, we establish the mathematical foundation of the model by proving the well-posedness, boundedness, and nonnegativity of solutions, as well as the existence of stationary positive solutions. Subsequently, we develop a stable and convergent alternating direction explicit-implicit scheme to numerically solve the nonlocal partial differential equations. Furthermore, we investigate the impact of key parameters on the evolution of tumor phenotypic heterogeneity, revealing that the epimutation rate governs the concentrated distribution of tumor cell phenotypes. Additionally, we leverage spatially resolved transcriptomics data to quantify tumor cell spatial density, phenotypic states, and vessel distribution within the model. Our integrative analysis uncovers significant associations between spatial transcriptomics features and the long-term progression of intratumor heterogeneity. This study provides a rigorous mathematical analysis of intratumor phenotypic heterogeneity evolution using a nonlocal reaction-diffusion model and identifies pivotal determinants of phenotypic heterogeneity through the integration of novel spatial data.</p>