Phenotypic heterogeneity of melanoma cells contributes to drug tolerance, increased metastasis, and immune evasion in patients with progressive disease. Diverse mechanisms have been individually reported to shape extensive intra- and inter- tumoral phenotypic heterogeneity, such as IFN-gamma signaling and proliferative-invasive transition, but how their crosstalk impacts tumor progression remains largely elusive. Here, using dynamical systems modeling and analysis of publicly available transcriptomic data at both bulk and single-cell levels, we demonstrate that the emergent dynamics of a regulatory network comprising MITF, SOX10, SOX9, JUN and ZEB1 can recapitulate experimental observations about the co-existence of diverse phenotypes (proliferative, neural crest-like, invasive) and reversible cell-state transitions among them. These phenotypes have varied levels of PD-L1, thus driving variability in immunosuppression. We elucidate how this heterogeneity in PD-L1 levels can be aggravated by combinatorial dynamics of these regulators with IFN-gamma signaling. Our model predictions are corroborated by analysis of PD-L1 levels in pre- and post-treatment scenarios both in vitro and in vivo. Our calibrated dynamical model offers a platform to test combinatorial therapies and provide rational avenues for clinical management of metastatic melanoma.