Abstract
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 high 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 tumor escape required co-occurrence of therapy resistance, motility, and immune suppression, with plasticity and multistability further promoting adaptive persistence. Lastly, we identified combination therapies predicted to constrain malignant diversification and enhanced immune accessibility. Taken together, our modeling work links regulatory dynamics with tumor-level adaptation and underscores potential strategies to therapeutically reprogram cell states toward sensitivity.