Bidirectional crosstalk between epithelial–mesenchymal plasticity and IFNγ-induced PD-L1 expression promotes tumour progression

Epithelial–mesenchymal transition (EMT) and immunoevasion through upregulation of programmed death-ligand 1 (PD-L1) are important drivers of cancer progression. While EMT has been proposed to facilitate PD-L1-mediated immunosuppression, molecular mechanisms of their interaction remain obscure. Here, we provide insight into these mechanisms by proposing a mathematical model that describes the crosstalk between EMT and interferon gamma (IFNγ)-induced PD-L1 expression. Our model shows that via interaction with microRNA-200 (miR-200), the multi-stability of the EMT regulatory circuit is mirrored in PD-L1 levels, which are further amplified by IFNγ stimulation. This IFNγ-mediated effect is most prominent for cells in a fully mesenchymal state and less strong for those in an epithelial or partially mesenchymal state. In addition, bidirectional crosstalk between miR-200 and PD-L1 implies that IFNγ stimulation allows cells to undergo EMT for lower amounts of inducing signal, and the presence of IFNγ accelerates EMT and decelerates mesenchymal–epithelial transition (MET). Overall, our model agrees with published findings and provides insight into possible mechanisms behind EMT-mediated immune evasion, and primary, adaptive, or acquired resistance to immunotherapy. Our model can be used as a starting point to explore additional crosstalk mechanisms, as an improved understanding of these mechanisms is indispensable for developing better diagnostic and therapeutic options for cancer patients.


Introduction
Activating invasion and metastasis and evading immune destruction are well-established hallmarks of cancer, complementary capabilities that enable tumour growth and metastatic dissemination [1]. Because metastasis is the main cause of cancer mortality, a thorough understanding of the interaction between these hallmarks is essential for developing therapeutic approaches, yet they are often studied separately [2]. Consequently, the interplay between metastatic dissemination and immune evasion remains poorly understood.
In recent years, epithelial-mesenchymal plasticity (EMP), the ability of cells to interconvert between intermediate E/M phenotypes along the epithelial-mesenchymal transition (EMT) spectrum [3], has been extensively studied because of its crucial role in invasion and metastasis (reviewed in [4][5][6]). Moreover, EMP is associated with therapeutic resistance (reviewed in [5,7]), including resistance to immunotherapy (reviewed in [8]). One of the mechanisms through which cancer cells acquire immune resistance is by activation of immune checkpoint pathways, which under physiological conditions are indispensable for self-tolerance and modulation of the immune response [9].
A frequently upregulated checkpoint protein in tumours is programmed death-ligand 1 (PD-L1) [10]. Its corresponding receptor programmed death 1 (PD-1) is expressed on the cell membrane of T cells, and PD-1-PD-L1 interaction induces a variety of immunosuppressive effects, such as inhibition of T cell proliferation, survival and effector functions [11]. Tumour cells can express PD-L1 by two general mechanisms: adaptive immune resistance, where PD-L1 expression is upregulated in response to inflammatory factors (such as interferon gamma (IFNγ)) produced by an anti-tumour immune response, and innate immune resistance, where PD-L1 expression is upregulated in response to constitutive oncogenic signalling [9]. Such constitutive signalling could, for example, be caused by EMT, and indeed several links between EMT and PD-L1-mediated immunoevasion have been reported in the literature (reviewed by Jiang and Zhan [12], who conclude that additional mechanistic studies are urgently needed).
Here, we study the post-transcriptional regulation of PD-L1 by the microRNA-200 (miR-200)-zinc finger E-box-binding homeobox 1 (ZEB1) axis, which is one of the proposed mechanisms underlying the interplay between EMT and immune resistance [13][14][15]. The binding of miR-200 to the mRNA of PD-L1 inhibits mRNA translation and in general such binding can stimulate miRNA decay [16,17]. The miR-200-ZEB1 axis is 'a motor of cellular plasticity' [18] and is considered to be part of the EMT core regulatory network [19]. Various mathematical models of EMT that include the miR-200-ZEB1 axis have been developed, which have contributed to a better mechanistic understanding of EMT (reviewed in [3,20,21]). Recently, Sahoo et al. [22] presented a mathematical model interconnecting a minimal EMT network and PD-L1 using shifted Hill functions to investigate immune-evasive strategies of hybrid E/M states, assuming an indirect effect of PD-L1 on EMT. Model analysis showed that both the stable hybrid and full EMT phenotypes resulting from the model were associated with high PD-L1 levels. However, although EMT scores determined from gene expression levels across a large panel of cell lines indeed exhibited a pattern of gradually increasing PD-L1 expression with the increasing EMT score, there was no clear dichotomy as expected from the mathematical model. This difference may be explained by the indirect nature of the feedback of PD-L1 to EMT as implemented in the model of Sahoo et al. [22]. Moreover, the impact of adaptive, IFNγ-driven PD-L1 expression, i.e. the primary reason for PD-L1 upregulation [23,24], was not taken into account in their model.
To describe the full interactions expected between adaptive and innate immune resistance and to study the impact of direct, mutual feedback between EMT and PD-L1, we connected a 'core' EMT model to a model for IFNγ-induced PD-L1 expression, which we developed based on an extension of a published JAK-STAT model [25]. For EMT, we used the ternary chimera switch (TCS) model by Jolly et al. [26], which is a simplification compared with the prior work [27]. Combining these two models was achieved by adding a mutually inhibitory feedback loop between miR-200 and PD-L1, which we described using appropriate miRNA-mRNA dynamics [28]. The analysis of our model shows that IFNγ-induced PD-L1 expression is expected to greatly accelerate EMT and decelerate the reverse mesenchymal-epithelial transition (MET) process. Moreover, IFNγ-induced PD-L1 lowers the required level of EMT-inducing signal, leading to an overall larger probability of EMT in tumours with a high IFNγ expression compared with tumours with a low expression. Vice versa, and consistent with the study by Sahoo et al. [22], a full EMT induced via other signals greatly upregulates PD-L1 expression, which IFNγ further amplifies. However, in our model, the hybrid EMT phenotype only moderately affects PD-L1 expression, which is due to the mutual feedback between EMT and PD-L1. Finally, we show that our model findings are broadly consistent with published experimental results in an extensive comparison with experimental data. Overall, our analysis illustrates how crosstalk between EMP and IFNγ-induced PD-L1 production can result in immune evasion and contribute to invasion and metastasis.

Modelling IFNγ-induced PD-L1 expression and EMT
To study the interplay between IFNγ-induced PD-L1 expression and EMT, we first need a quantitative description of these major processes. Although interferon gamma (IFNγ)-induced PD-L1 expression via the JAK-STAT pathway (figure 1a) has been extensively studied [29,30], to the best of our knowledge, no mathematical model of this entire process exists. Existing models focus on IFNγinduced JAK-STAT signalling (reviewed in [31]) and usually do not take into account direct targets of STAT such as interferon regulatory factors (IRFs) and interferon-stimulated genes. To model IFNγinduced PD-L1 expression, we used the simplified JAK-STAT model by Quaiser and Mönnigmann [32] (figure 1b, left box), which is a truncated version of a more detailed model by Yamada et al. [33]. In this model, phosphorylated signal transducer and activator of transcription 1 (STAT1) homodimer (STAT1p_2) reaches its maximum at about 30 min after IFNγ exposure (figure 1c), in agreement with  Quaiser et al. [25] simplified TCS. Jolly et al. [26] promotor GAS P inducing signal (e.g. TGFβ) The eventual phenotype also depends on the initial state of the system, illustrated by temporal dynamics after a sudden increase in SNAIL1 from 2.0 × 10 5 molecules to 2.15 × 10 5 molecules given three different initial states (coloured lines).
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 220186 Brysha et al. [34]. This time scale is much shorter than the time scale of hours at which PD-L1 dynamics occur (figure 1d), in agreement with Kleijn et al. [35] and Ghosh et al. [36]. Therefore, we simplified this model further by describing the steady state of its output, i.e. STAT1p_2, with a Gompertz function (electronic supplementary material, Methods). We then extended this model with the production of IRF1 mediated by STAT1, and production of PD-L1 mediated by IRF1 (figure 1b, middle box, and electronic supplementary material, figure S1A). We modelled both of these processes with shifted Hill functions, which led to realistic IRF1 [37] and PD-L1 [36] dynamics (electronic supplementary material, figure S1B).
For EMT, various mathematical models have been developed (reviewed in [20,21]). Two well-known models of the miR34-SNAIL1 and miR200-ZEB1 axes, i.e. the core regulatory machinery of EMT [19], are the cascading bistable switches model (Tian et al. [38]; later revised by Zhang et al. [39]) and the TCS model [27]. We elected to use the simplified TCS model developed by Jolly et al.

IFNγ amplifies the increase in PD-L1 caused by EMT primarily for mesenchymal cells
Because we aimed to study the impact of IFNγ-mediated PD-L1 on EMT and vice versa, we connected the separate model parts (figure 1b, boxes) by adding mutual inhibition between PD-L1 and miR-200 [13,14,40] using the theoretical framework for miRNA-transcription factor (TF) dynamics by Lu et al. [28]. This results in direct inhibition of PD-L1 by miR-200, as miR-200 binds to the mRNA of PD-L1 [13]. However, there is automatically also an inhibitory effect of PD-L1 on miR-200 because the modelled miRNA-mRNA interaction leads to additional miR-200 degradation (Lu et al. [28], see electronic supplementary material, Methods for details). This PD-L1 inhibitory effect on miR-200 also affects the miR-200 level in the absence of IFNγ because in the coupled model PD-L1 has a non-zero expression level (driven by EMT status). The basal miR-200 degradation rate in the TCS model should contain the influence of degradation by any molecule not described explicitly in that model, including PD-L1. Thus, to aim for a fair comparison between these structurally different models (i.e. TCS and our combined model (figure 1b)), we decreased the basal miR-200 degradation rate in the combined model until it approximately matched the bifurcation diagram for the TCS model (electronic supplementary material, figure S2).
Next, we studied how the system responds to different levels of SNAIL1 (activated via, for example, transforming growth factor beta (TGFβ)) and IFNγ (figure 2). In the absence of IFNγ, the previously reported tristable ZEB1 levels [27] translate to similar tristability in PD-L1 expression at the cell membrane (blue line in figure 2a). Here, the highest PD-L1 expression occurs for mesenchymal cells and the lowest for epithelial cells. Note that the PD-L1 expression of hybrid cells is only slightly increased compared with epithelial cells, a prediction that is different from the recent study by Sahoo et al. [22], who predict an almost equal level of PD-L1 for the hybrid and mesenchymal phenotype. Notably, the mutual inhibition of miR-200 and PD-L1 did not result in additional bifurcation points. However, the combined model predicts that exposure of cells to IFNγ greatly amplifies PD-L1 expression for all phenotypes, thereby also amplifying the differences in PD-L1 expression between those EMT phenotypes (figure 2a, orange line; figure 2b). Nevertheless, in the presence of IFNγ, the hybrid phenotype has a PD-L1 expression level similar to that of epithelial cells, and only mesenchymal cells are expected to have substantially higher PD-L1 levels.

IFNγ promotes the occurrence of EMT
We next asked what the impact of IFNγ is on ZEB1, i.e. a central EMT-TF. In the absence of IFNγ, cells in our model undergo EMT for similar levels of SNAIL1 compared with the original TCS model (blue curve in figure 2c compared with the grey curve). Small differences between our model and the TCS model are royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 220186 caused by the bidirectional influence between PD-L1 and miR-200. In the presence of IFNγ, this mutual influence causes the bifurcation diagram to shift to the left (figure 2c, orange curve), because the expression of PD-L1 leads to a reduced amount of miR-200, in turn affecting EMT. This implies that according to our model IFNγ promotes EMT because both partial and full EMT more readily occur than in the absence of IFNγ.
We created a phase diagram to provide an overview of how the stability of EMT phenotypes depend on IFNγ and SNAIL1 levels (figure 2d ). In this phase diagram, the IFNγ-induced leftward shift is clearly visible. In addition, it shows that the total range of SNAIL1 for which the hybrid E/M phenotype, a particularly aggressive EMT phenotype (Lüönd et al. [41], and reviewed in Jolly et al. [42]), can (co-)exist is unaffected. As the different EMT phenotypes express different levels of PD-L1 we additionally quantified these levels for each phenotype and for the average PD-L1 expression per stability region (electronic supplementary material, figure S3).
To investigate the role of the included JAK-STAT regulation, we created a simplified crosstalk model where we replaced JAK-STAT signalling with a generic inducing signal I (electronic supplementary material, figure S4A). The influence of this inducing signal is more gradual than the influence of IFNγ in the full model, yet the same leftward shift is present (electronic supplementary material, figure S4B, compared with figure 2b), suggesting that IFNγ-induced JAK-STAT signalling makes the response highly 'switch-like'. In addition, we used a sensitivity analysis on parameters used in the negativefeedback loop between miR-200 and PD-L1, which showed that both the predicted leftward shift of SNAIL1 levels and the amplification of PD-L1 by IFNγ are robust to parameter change (electronic supplementary material, figure S5). In conclusion, our model predicts that local presence of IFNγ royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 220186 promotes partial and full EMT by lowering the threshold for additional EMT-inducing signals to induce (partial) EMT.

IFNγ-driven PD-L1 expression accelerates EMT
Besides the steady states for EMT status and PD-L1 expression that the combined regulatory network can reach in the long run, in practice, it will also matter how long it takes to reach these states. Therefore, we also investigated our combined regulatory network's temporal dynamics and compared them with the temporal dynamics of the separate models.
To study the impact of PD-L1 expression on the temporal dynamics of an EMT transition and vice versa, we started with a stable epithelial phenotype (SNAIL1 = 1.7 × 10 5 molecules) and simulated transition to a fully mesenchymal phenotype by a sudden increase of SNAIL1 (to 2.3 × 10 5 molecules). In the absence of IFNγ, there is only a minor change in the temporal dynamics of EMT in our combined model relative to the simplified TCS model (blue versus grey lines in figure 3). However, as a consequence of the double-negative feedback loop between miR-200 and PD-L1, the decrease in miR-200 during EMT affects the PD-L1 mRNA and membrane PD-L1 (which are not present in the TCS model) at a similar time scale as miR-200. Thus, the effect of EMT on PD-L1 expression takes place on a much slower time scale than for PD-L1 expression driven by IFNγ (for comparison, see figure 1d).
In the presence IFNγ, the increased level of PD-L1 over time speeds up the miR-200 decrease (compared with the TCS model that does not consider an interaction with PD-L1). This is due to fast PD-L1-driven miR-200 degradation during the initial phase of EMT (electronic supplementary material, figure S6). Because of miR-200's role as a critical suppressor of several EMT-TFs [43,44], this  royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 220186 accelerates the EMT process by about 10 days, mainly by reducing the time cells remain in a state close to the hybrid (E/M) phenotype before converting to a fully mesenchymal phenotype (orange line in figure 3, left). Again, the PD-L1 expression level evolves slowly compared with the rapid increase observed on a time scale of hours for the model describing just IFNγ signalling (figure 1d). This difference in time scales of initial IFNγ-driven PD-L1 increase and of further PD-L1 increase due to EMT is especially apparent when considering a scenario where cells are simultaneously exposed to IFNγ and a TGFβ (electronic supplementary material, figure S7, left column).
In summary, the double-negative feedback loop between miR-200 and PD-L1 accelerates EMT in the presence of IFNγ.

IFNγ-driven PD-L1 expression decelerates MET
Finally, we investigated whether similar changes occur in the temporal dynamics of the reverse process in which mesenchymal cells transition to an epithelial phenotype, i.e. MET. For this purpose, we started with a stable mesenchymal phenotype, where SNAIL1 = 2.3 × 10 5 molecules, which we instantaneously decreased to 1.7 × 10 5 molecules. In the simplified TCS model, this SNAIL1 decrease led to a direct transition from the mesenchymal into the epithelial phenotype (grey dashed line in figure 3, right), i.e. without a substantial slowing down of the dynamics around a hybrid phenotype.
As anticipated from the observed EMT dynamics showing that the double-negative feedback loop between PD-L1 and miR-200 favours the mesenchymal phenotype, this loop also affects MET, albeit in the opposite direction, i.e. it decelerates MET. This is primarily due to the long-lasting suppression of miR-200 (orange lines in figure 3, right). However, this effect is considerably smaller than the effect in the forward transition, which is probably the result of MET not transitioning through the hybrid state at a slow pace since the PD-L1-miR-200 interaction primarily accelerated EMT by shortening the time spent there ( figure 3, left). Note that simultaneous lowering of SNAIL1 and removal of IFNγ do lead to an MET occurring at a similar speed as in the original TCS model (electronic supplementary material, figure S7, right column), demonstrating the important role of IFNγ in EMP.
In conclusion, the PD-L1-miR-200 double-negative feedback loop is expected to slow down MET considerably in the presence of IFNγ.

Discussion and conclusion
Here, we have used a mathematical model to describe the crosstalk between IFNγ-induced PD-L1 expression and EMT. We showed that merely adding the reported interaction between PD-L1 and miR-200 gives rise to tristability in the PD-L1 levels, where a mesenchymal state corresponds with high PD-L1 expression and an epithelial state with low PD-L1 expression. The difference in PD-L1 levels between the stable EMT states is amplified by adding IFNγ stimulation. In addition, we showed that adding this crosstalk reduces the amount of SNAIL1 required to undergo EMT in the presence of IFNγ. Finally, we showed that in the presence of IFNγ, this crosstalk accelerates the forward EMT process and decelerates the reverse MET process.
To assess the extent to which the model results agree with recent studies, we have compiled a summary of published experimental reports on the link between EMT and PD-L1 and indicated which findings are qualitatively consistent with our results (table 1). Note that we omitted several papers that propose different mechanisms underlying EMT/PD-L1 crosstalk and lacked experimental findings that could be compared with our simulations [54][55][56][57]. Overall, the simulated results of our model are in good agreement with the reported experimental findings.
Although findings of our model and experimental data are in good agreement, some results are partially conflicting: firstly, Noman et al. [14] reported that while PD-L1 expression in EMT-activated breast cancer cell lines is regulated via the SNAIL1/ZEB1/miR-200 axis, treatment of MCF7-2101 cells with TGFβ or inhibition of TGFβ in MCF sh-WISP2 cells did not modulate PD-L1 expression. However, the MCF7-2101 cell line already has a mesenchymal phenotype, evident from their vimentin expression and lack of E-cadherin; hence, additional TGFβ is not expected to further increase PD-L1 expression by triggering EMT. Moreover, although the TGFβ signalling inhibitors repressed SMAD2 activation in MCF7 sh-WISP2 cells, phosphorylated SMAD2/3 is only essential for SNAIL1 activation and not for the persistence of SNAIL1 expression [58]. Thus, the performed experiments cannot rule out that TGFβ signalling modulates PD-L1 expression in these cell lines. Secondly, PD-L1 expression royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 220186 Table 1. Summary of experimental reports regarding the crosstalk between EMT and PD-L1 and their consistency with our ODE model.

cancer type description match
association between EMT markers and PD-L1 expression in different cancer types [12] yes breast cancer Induction of EMT in human breast cells using TGFβ upregulated PD-L1 expression in a manner that was dependent on the activation of the PI3K/AKT and MEK/ERK pathways. In addition, manipulation of PD-L1 modulated the EMT status of breast cancer cells. Downregulation of PD-L1 using specific shRNA in mesenchymal-like cells led to partial EMT reversal. Conversely, overexpression of PD-L1 led to EMT [40]. Mechanistically, PD-L1 activated the EMT process via Ras binding and activating the downstream MEK/ERK signalling [48]. in certain cell lines is reportedly affected by silencing of ZEB1 but not of SNAIL1 [14,45]. Although SNAIL1 drives ZEB1 expression in the EMT core regulatory network model we used, this apparently contradictory finding can be reconciled by appreciating that the EMT transcriptional response is highly context specific [59,60]. Indeed, sustained expression of ZEB1 can be achieved by various EMT regulators and even by ZEB1 itself, resulting in an irreversible switch to the mesenchymal phenotype [39,61,62]. Thus, SNAIL1 silencing may be insufficient to affect PD-L1 in settings where ZEB1 is maintained via other means. Even though partially conflicting findings between our model and published data can be reconciled, it is also likely that additional mechanisms by which PD-L1 and EMT mutually influence each other play a role depending on the studied cell line or cancer. For instance, PD-L1 expression can induce EMT by activating the TF SREBP-1c in hepatocellular and renal cell carcinoma [49,53], or by preventing SNAIL1 ubiquitination in triple-negative breast cancer [63]. In addition, ZEB1 can bind and silence IRF1, a TF of PD-L1 [64], and the activity of the CMTM family, which stabilizes PD-L1, correlates with EMT-TFs such as SLUG [65,66].
In the related mathematical model presented by Sahoo et al. [22], an indirect feedback mechanism of PD-L1-mediated E/cadherin inhibition was included. Results of their model agree qualitatively with those of our model in terms of the prediction that hybrid E/M cells are expected to be PD-L1 positive and hence immune-evasive. However, the predicted PD-L1 expression levels for the different EMT phenotypes differ substantially: Sahoo et al. [22] predict that the hybrid and mesenchymal phenotypes have almost equally high PD-L1 expression, whereas we predict that hybrid cells are closer to epithelial cells with respect to PD-L1 expression. They substantiated that PD-L1 expression increases with EMT scores by an analysis of gene expression in pan-cancer datasets and analysis of responses in several cell lines. For some of these data, intermediate EMT scores (interpreted as hybrid cells) have the tendency to only moderately increase the PD-L1 expression (e.g. figs 1E and F and 5E in [22]),  royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 220186 9 which is more in line with our model predictions. A key difference in our respective approaches is that Sahoo et al. [22] employed RACIPE (random circuit perturbation) [67], which yields semiquantitative predictions based on network topology, whereas we used a mechanistic approach based on appropriate miRNA-mRNA dynamics. As a consequence, only the inhibition of miR-200 on PD-L1 is present in the model by Sahoo et al. [22], because the inhibition of PD-L1 by miR-200, which arises as a consequence of miRNA-mRNA dynamics, cannot be captured by using the Hill-type dynamics in RACIPE. Nevertheless, the experimental work is required to investigate whether miR-200 degradation is enhanced upon binding to PD-L1. Moreover, the involved mechanisms and factors affecting PD-L1 expression and EMT are probably context and cell-line specific, and further research will be required to unravel this complicated interplay between factors. Our model forms a basis to include such factors.
A limitation of our model is its relatively crude description of IFNγ-induced PD-L1 expression. To keep the model simple, we used a steady-state approximation of the simplified model by Quaiser et al. [25], which only describes the first 15 min of JAK-STAT signalling before transcriptional feedback occurs. Nevertheless, our model's predicted IRF1 temporal dynamics roughly agree with the JAK-STAT model by Rateitschak et al. (cf. fig. 4A in [37]). Comparison of these results with those of a highly simplified IFNγ-driven PD-L1 expression model showed that the role of JAK-STAT signalling primarily creates a 'switchlike' response to IFNγ (cf. figure 2b and electronic supplementary material, figure S4b).
The predicted decrease of SNAIL1 levels required for EMT in the presence of IFNγ (as indicated by the leftward shift in the bifurcation diagrams in figure 2c,d) is noteworthy in itself; although several studies report on the modulation of PD-L1 expression as a result of EMP, the mechanism by which PD-L1 could modulate EMT remains unclear [12]. Here, we show that the additional degradation of the PD-L1-miR-200 complex in the presence of IFNγ could be affecting miR-200 levels sufficiently to affect EMT, without the need for additional mechanisms.
Another interesting result is the acceleration of EMT by IFNγ as a result of the increased interaction between miR-200 and PD-L1. This is reminiscent of recent theoretical work that showed that intermediate EMT states or unobservable microstates might accelerate the EMT process [68,69], i.e. an increased PD-L1 level could be viewed as such a microstate. (Single-cell) transcriptomics experiments with sufficient time resolution would help decipher the different involved mechanisms and how they influence the EMT temporal dynamics.
The results reported here may have diagnostic and therapeutic implications. PD-L1 has been proposed as a biomarker to predict the efficacy of PD-1/PD-L1 blockade therapy [70], yet the multifactorial mechanisms underlying PD-L1 membrane expression complicate its use as an exclusive biomarker [71,72]. For example, oncogene-driven PD-L1 expression caused by EMT is constitutive, diffuse and distinct from inflammation-driven PD-L1 expression, which may occur more focally and during a limited time window. The latter is often associated with the presence of an immune infiltrate, while the former is associated with a lack thereof [71], and it has been proposed that this combination of PD-L1 expression and the presence of immune infiltrate affects the response to PD-1/ PD-L1 blockade therapy [73]. Our findings highlight how oncogenic and inflammation-driven PD-L1 expression might interact and give rise to strongly increased PD-L1 levels. PD-L1 on tumour cells is an important contributor to immune evasion through inhibition of CD8+ T cell cytotoxicity [74]. Consistently, in metastatic urothelial cancer, lack of response to anti-PD-L1 treatment occurred particularly in patients with tumours that showed exclusion of CD8+ T cells [75]. Combined treatment with TGFβ-blocking and anti-PD-L1 antibodies invoked anti-tumour immunity and tumour regression by facilitating T cell infiltration, which could be due to EMT-independent TGFβ signalling [76], but possibly also via the EMT-dependent path we have presented here. Another opportunity could be the use of therapeutic siRNA [77]: in contrast to anti-PD-L1 antibodies, siRNA affects PD-L1 mRNA, which, according to our model, should free up miR-200 for its anti-EMT effect, a prediction that is supported by ( partial) EMT reversal following downregulation of PD-L1 using siRNA/shRNA [40,52].
In conclusion, we developed a mathematical model that describes the crosstalk between IFNγinduced PD-L1 and EMT, which is in good agreement with experimental findings. In addition, our model sheds light on potential mechanisms behind EMT-mediated immune evasion, and primary, adaptive, or acquired resistance to immunotherapy. Our (simplified) model can serve as a starting point to explore additional EMP and immune crosstalk mechanisms. In particular, we propose embedding the presented model in a multi-scale model to explicitly describe the local effects of the adaptive immune response and the effects of TGFβ on the tumour microenvironment. Improved royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 220186 understanding of the interaction between the immune response and EMP is indispensable for developing better diagnostic and therapeutic options for cancer patients.

Simplified ternary chimera switch model
The simplified TCS model [26] is built on the theoretical framework for microRNA-TF chimera toggleswitches defined in the study by Lu et al. [28]. See electronic supplementary material, Methodsfor the model definition, theoretical framework, and used parameters.

IFNγ-PD-L1 model
The IFNγ-PD-L1 model consists of two parts: (IFNγ-)JAK-STAT signalling, for which we use a steadystate approximation of an existing model, and STAT-PD-L1 regulation which we developed following the theoretical framework underpinning the simplified TCS model. We discuss these parts in more detail below.

(IFNγ-)JAK-STAT
For the IFNγ-PD-L1 model, we use the JAK-STAT model presented in [25] (see electronic supplementary material, Methods). The dynamics of IFNγ-STAT and STAT-PD-L1 are well separated in time (figure 1c, d), with JAK-STAT dynamics happening on a time scale of minutes [34] and STAT-PD-L1 dynamics happening on a time scale of hours [35,36]. Therefore, we simplified our IFNγ-PD-L1 model by assuming the JAK-STAT model is in aquasi-steady state; we fitted its steady state by a Gompertz function (see electronic supplementary material, figure S8) and inserted this approximate relationship between IFNγ and STAT in our STAT-PD-L1 model.

STAT-PD-L1
The STAT-PD-L1 submodel (see also electronic supplementary material, figure S1A) consists of the following ordinary differential equations (ODEs): This model uses the abundance of STAT (in molecules) as input and includes IRF1 and PD-L1 using appropriate TF-TF dynamics (H S ) from the theoretical framework by Lu et al. [28] (see electronic supplementary material, Methods). Transport of PD-L1 through various cellular compartments is modelled with a constant rate from compartment to compartment, with rates assumed to be similar to Lippincott-Schwartz et al. [78]. Model parameters are provided in table 2; some parameters were considered to be similar to those of the simplified TCS model.

Combined model
To create our combined model, we connected the IFNγ-induced JAK-STAT signalling model and simplified TCS model to the central STAT-PD-L1 model (see figure 1b). In this combined model, PD-L1 mRNA : dm P dt ¼ g mP H S ðF, l F,mP Þ À m P Y m ðmÞ À k mP m P ; ð4:7Þ PD-L1 in ER : dP ER dt ¼ g PER m P LðmÞ À k ER,G P ER ð4:8Þ and miR-200 : dm dt ¼ g m H S ðZ, l Z,m ÞH S ðS, l S,m Þ À m Z Y m,mZ ðmÞ À m P Y m,mP ðmÞ À k m m: ð4:9Þ The connection between miR-200 and PD-L1 is done using appropriate miRNA-mRNA dynamics (L and Y functions) from the theoretical framework by Lu et al. [28] (see electronic supplementary material, Methods). To accommodate a low number of binding sites for miR-200 on the mRNA of PD-L1 (i.e. two binding sites according to Chen et al. [13]), we adapted the parameters for L(μ), Y m (μ) and Y m,mP ðmÞ from electronic supplementary material, table S1 (n = 6) to the values as shown in table 3 (n = 2). Note that the inhibition of miR-200 on PD-L1 automatically leads to inhibition of PD-L1 on miR-200 as the included miRNA-mRNA dynamics [28] result in additional miR-200 degradation (m P Y m,mP ðmÞ in equation (4.9); see electronic supplementary material, figure S6). We approximately matched the results of our combined model in the absence of IFNγ with the results of the TCS model by compensating for this additional degradation of miR-200. We achieved this by decreasing the basal degradation rate of miR-200 from k m ¼ 0:05 (electronic supplementary material, table S2) to k m ¼ 0:0475 (see electronic supplementary material, figure S2). We considered μ 0 = 10 4 molecules, following Jolly et al. [26]. In addition, since the JAK-STAT model has [STAT1p_2] as output in nM we need to convert to number of molecules to use STAT1p_2 as input in the STAT-PD-L1 model. As in the study by Jolly et al. [26], we use a cell volume of 10 000 μm 3 , such that 1 nM amounts to roughly 6020 molecules (6.02 × 10 23 · 10 −9 · 10 000 × (10 −5 ) 3 ). Note that this cell volume is on the high side, as typical animal cells are 10-20 μm in diameter (approx. 500-4000 μm 3 ) [79, p. 529]). However, because of our IFNγ-STAT steady-state approximation, a decrease in cell volume by a factor 10 can be compensated by multiplying the inducing Table 2. Variables and parameters used for the STAT-PD-L1 model. The top panel shows variable names and production and degradation rates; the bottom panel shows parameters for the shifted Hill functions of the interactions. Degradation rates k are in h −1 , production rates g in molecules h −1 and thresholds B 0 A in molecules.