Rahman, N.; Jackson, T. L.

HER2 heterogeneity and reversible phenotypic plasticity play a central role in breast cancer progression and therapeutic resistance, yet how their interaction shapes treatment response remains poorly understood. Experimental and clinical evidence indicates that HER2-positive and HER2-negative tumor cell states can dynamically interconvert, enabling compensatory population shifts that undermine monotherapies targeting a single phenotype. Because stochastic lineage effects, phenotypic witching, and local cell interactions are averaged out in mean-field population-level ODE models, we develop a spatially resolved agent-based model (ABM) that explicitly represents individual-cell dynamics and heterogeneous tumor growth. The model incorporates phenotype-specific proliferation, migration, and death, division-coupled HER2 state transitions, and therapy-induced selective pressures. We consider two therapeutic interventions with complementary mechanisms of action: paclitaxel, which preferentially suppresses HER2-positive proliferation, and Notch inhibition, which targets HER2-negative populations and alters phenotypic composition. Starting from single-cell lineages, we validate the ABM against theoretical predictions from a population-level switching model and against single-cell-derived experimental measurements, demonstrating quantitative agreement with early lineage dynamics and long-term phenotypic equilibria. Simulation results show that monotherapies induce compensatory phenotypic shifts and spatial reorganization that permit tumor persistence. In contrast, combination therapy simultaneously targeting HER2-positive and HER2-negative populations disrupts phenotypic replenishment, fragments spatial structure, and can achieve sustained tumor control across a broad range of treatment strengths. To quantify robustness across heterogeneous tumor parameter regimes, we pair the ABM with an interpretable Random Forest surrogate trained on ensemble simulation data. Using only pre-treatment and early-trajectory features, the surrogate predicts long-term response, identifies growth-rate asymmetries as dominant drivers of resistance, and interpolates across previously unseen parameter combinations within the sampled domain. Together, this integrated mechanistic and data-driven framework clarifies how HER2-mediated plasticity, spatial organization, and competitive growth dynamics shape therapy resistance and provides a scalable approach for predicting and optimizing treatment strategies in HER2-heterogeneous breast cancer. Sample MATLAB code for the agent-based model (ABM) used in this study is available on GitHub.