Regulatory network and spatial modeling reveal cooperative mechanisms of resistance and immune escape in ER+ breast cancer.

Despite significant progress, the treatment of estrogen receptor-positive (ER+) breast cancer remains clinically challenging due to reversible drug resistance and immune evasion. Drug resistance often arises as cells undergo a dynamic epithelial-to-mesenchymal transition (EMT), while elevated PD-L1 levels contribute to immune escape. While these phenotypic features can variably co-occur, the impact of co-occurrence on the availability of synergistic treatment strategies remains unknown. To investigate their interplay, we constructed an ER-EMT-PD-L1 gene regulatory network and simulated these networks as coupled ordinary differential equations with biologically informed parameters, to generate steady-state expression profiles. Our study revealed that the relevant overarching network generated antagonistic epithelial and mesenchymal modules, capable of producing monostable, bistable, and tristable dynamics. We further examined the link between phenotypes and immune evasion by quantifying average PD-L1 expression, and found that epithelial-sensitive states consistently exhibited low PD-L1. In contrast, hybrid- and mesenchymal-resistant states were associated with a non-linear, stepwise increase in PD-L1, highlighting a strong coupling between EMT, resistance, and immune evasion. Extending on these network-level insights, we further used a spatially explicit agent-based model seeded with GRN-derived phenotypes to probe tumor behavior under therapeutic pressure. Simulations revealed that sustained tumor expansion occurred only when resistance, motility, and immune evasion traits co-existed, and this requirement remained robust across GRN landscapes with differing stability. Plasticity and multistability increased the accessible phenotypic state-space and accelerated shifts toward high-fitness resistant states. We further identified combination therapies that significantly reduced phenotypic diversification and improved immune infiltration in silico. Taken together, our modeling work links regulatory dynamics with tumor-level adaptation and highlights strategies to reprogram resistant cell states toward sensitivity, which are difficult to infer from bulk or cross-sectional data alone. In addition, it provides a controllable in silico testbed to systematically evaluate candidate treatment combinations and their effects on tumor phenotypic transitions and spatial T cell access, thereby helping to prioritize experimental regimens for follow-up.