Tumor-mediated immunosuppression and cytokine spreading affects the relation between EMT and PD-L1 status

Epithelial-mesenchymal transition (EMT) and immune resistance mediated by Programmed Death-Ligand 1 (PD-L1) upregulation are established drivers of tumor progression. Their bi-directional crosstalk has been proposed to facilitate tumor immunoevasion, yet the impact of immunosuppression and spatial heterogeneity on the interplay between these processes remains to be characterized. Here we study the role of these factors using mathematical and spatial models. We first designed models incorporating immunosuppressive effects on T cells mediated via PD-L1 and the EMT-inducing cytokine Transforming Growth Factor beta (TGFβ). Our models predict that PD-L1-mediated immunosuppression merely reduces the difference in PD-L1 levels between EMT states, while TGFβ-mediated suppression also causes PD-L1 expression to correlate negatively with TGFβ within each EMT phenotype. We subsequently embedded the models in multi-scale spatial simulations to explicitly describe heterogeneity in cytokine levels and intratumoral heterogeneity. Our multi-scale models show that Interferon gamma (IFNγ)-induced partial EMT of a tumor cell subpopulation can provide some, albeit limited protection to bystander tumor cells. Moreover, our simulations show that the true relationship between EMT status and PD-L1 expression may be hidden at the population level, highlighting the importance of studying EMT and PD-L1 status at the single-cell level. Our findings deepen the understanding of the interactions between EMT and the immune response, which is crucial for developing novel diagnostics and therapeutics for cancer patients.

where µ is the amount of miRNA, n represents the number of binding sites, l i is the translation rate of mRNA bound to i miRNA molecules, and C i n and M i n (µ) are calculated with In the latter equation, µ 0 represents the equilibrium constant for the binding and unbinding of miRNA to mRNA. µ 0 is calculated with where r µ+ and r µ− are respectively the binding and unbinding rate of miRNA to mRNA. The total translation of mRNA is calculated with where m 0 represents the total mRNA concentration. The miRNA-facilitated mRNA degradation rate is calculated with where γ mi is the individual degradation rate for an mRNA bound to i miRNA. Similarly, the miRNA-facilitated miRNA degradation rate is calculated with where γ µi represents the individual degradation rate for a miRNA [see Supplementary Information of 1, for derivation and more details].

IFNγ-PD-L1-EMT model
The IFNγ-PD-L1-EMT model [2] combines the simplified TCS model [3] with a model for IFNγ-induced PD-L1 expression, which is based on an extension of a published JAK-STAT model [4]. The combined model consists of the following ODEs: PD-L1 mRNA : This model uses SNAIL1 and IFNγ as inputs and includes appropriate TF-TF dynamics (H S ) and miRNA-mRNA dynamics (L and Y functions) from the theoretical framework by Lu et al. [1]. The parameters used in the L(µ), Y m (µ), and Y µ (µ) functions are shown in Table S1 [following 2]. Table S2 lists the variables and remaining parameters used in the model. For (IFNγ)-JAK-STAT signaling, the IFNγ-PD-L1-EMT model uses a steady-state approximation of the JAK-STAT model by Quaiser et al. [4]. The relationship between IFNγ and [STAT1p_2] is described with a Gompertz function: with parameter values d = 414.56394, b = −51.52794, and e = 0.02833.
[STAT1p_2] in nm is converted to STAT1p_2 in number of molecules by multiplying with 6020 following Jolly et al. [3]. Table S3 shows all components in the IFNγ-PD-L1-EMT model and their units.

Cascading Bistable Switches model
The revised CBS model ( [5]; [originally presented by 6]) is built on a theoretical framework for miRNA-mediated regulation of mRNA different from the one defined in Lu et al.    (2) and (3)). General miRNA, mRNA, and protein dynamics are modeled with Protein : where k and kd denote respectively production and degradation rates, [miRNA] and [mRNA] are the amount of free miRNA and mRNA, respectively, λ is the recycle ratio for an miRNA following degradation of the miRNA-mRNA complex R1, and k s0 and k s1 are the translation rates of free mRNA and the R1 complex, respectively. R1 is calculated with In this equation, the equilibrium constant K is calculated with where k on and k off are the binding and unbinding rates of miRNA to mRNA, respectively. The amount of miRNA-mRNA complex with i miRNAs bound is calculated with Conservation equations for miRNA and mRNA are where n is the number of binding sites, and C i n is calculated with Eq. (S2) [see Supplementary Information of 5, for derivation and more details].
The TGFβ and SNAIL1-miR-34 modules of the revised CBS model consist of the following equations: Free TGFβ mRNA : Total TGFβ protein : In the TGFβ module, [T GF 0] represents the input concentration of exogenous TGFβ in µm. Table S4 lists the variables and parameters used in the TGFβ and SNAIL1-miR-34 modules.

IFNγ
Because PD-L1 expression on the tumor cell membrane should depend on the local IFNγ level, we simulated IFNγ diffusion in our multi-scale models using a PDE layer describing the IFNγ concentration. Within the Morpheus framework, secretion and degradation rates are expressed per lattice site [7]. Hence, we divided the T cell production and tumor cell uptake rates of IFNγ (explained below) by the actual area of T cells and the target area of tumor cells, respectively. IFNγ is primarily secreted by activated lymphocytes such as CD8 + and CD4 + T cells, natural killer (NK) cells, and NK T cells [8]. For simplicity, we considered tumor-infiltrating CD8 + T cells to be the only source of IFNγ in our model. T cells secrete IFNγ at an average rate of 1200 molecules min −1 for several hours after activation [9]. Therefore, we set the baseline production of IFNγ by T cells at this rate. In simulations with PD-L1-mediated inhibition of IFNγ, the total IFNγ production rate of a T cell depended on the average PD-L1 level of its neighboring lattice sites, following Eq. (4) in the main text. In simulations with TGFβ-mediated inhibition of IFNγ, the total IFNγ production rate of a T cell depended on its local TGFβ concentration, calculated as described in the next section, following Eq. (8) in the main text. Produced IFNγ was distributed over all lattice sites covered by the generating T cell. Upon secretion, IFNγ diffused into the TME with a diffusion coefficient of 5.43 × 10 3 µm 2 min −1 , which was previously measured in vivo in murine lymph node tissue [10]. As it remains unclear to what extent T cell-derived IFNγ can spread within the TME, and this likely depends on TME composition [11], we simulated two different extents of IFNγ spreading by modifying the rate of cellular uptake of IFNγ by tumor cells. This rate has previously been estimated at 2.1 min −1 based on an in vitro experiment with cultured fibroblasts [12,13]. We opted to consider differential uptake rates to describe the two spreading scenarios because the uptake rate may well be different in in vivo tumors compared to in vitro fibroblasts. We modeled the long-and short-range spreading of IFNγ with uptake rates of 2.1 × 10 −2 min −1 and 2.1 × 10 3 min −1 , respectively. Using these rates, the IFNγ concentration in molecules cell −1 decreases by a factor of 2.7 within 6 and 1 cell layer(s), respectively, when simulating a single IFNγ-producing T cell surrounded by IFNγ-consuming tumor cells. Following Beck et al. [13], we set the rate of disappearance of IFNγ in the medium at 20% of the cellular uptake rate of tumor cells. We considered T cells not to consume IFNγ. The Morpheus framework employs a Forward-Euler diffusion solver that has a time-step limit depending on the diffusion rate [7]. To make our simulations feasible with regard to the computation time involved, we decreased the diffusion coefficient and cellular uptake rates by a factor of 100 (i.e., the earlier-mentioned values for these parameters were divided by 100). Importantly, this only influences the absolute IFNγ concentration and the time until a steady state is reached, and not the extent of IFNγ spreading.
We approximated the local IFNγ concentrations in nm cell −1 by dividing the total number of molecules located on the lattice sites belonging to a cell by the volume that the cell would have in three dimensions based on a diameter of 24 µm [14] and considering tumor cells to be perfect spheres. This is the diameter of a circular tumor cell with an area that equals the tumor cell target area. The resulting IFNγ concentration was used as an input signal for the ODE model. For short-range IFNγ spreading simulations, the calculated IFNγ concentrations were lowered (by multiplying them with 0.6) to create a concentration range that led to a limited amount of phenotype-switching events roughly consistent with the low numbers of E/M and M tumor cells in in vivo tumors. We employed the Neumann boundary condition representing zero flux at the boundary in all simulations. For long-range IFNγ spreading simulations, we additionally set the IFNγ disappearance rate in lattice sites along the simulation space boundary to 1 min −1 to prevent IFNγ accumulation.

TGFβ
Simulations either had a uniform TGFβ field or a TGFβ gradient. We modeled the former by setting the TGFβ concentration in nm cell −1 as a fixed parameter instead of explicitly modeling a field. The TGFβ concentration was 0.07 nm cell −1 and 0.1 nm cell −1 for long-and short-range IFNγ spreading simulations, respectively. These parameters were taken differently in order to tune the phenotype switching rates between these two scenarios. We modeled the TGFβ gradient with a TGFβ field varying between approximately 13.7 and 2.1 molecules per lattice site from left to right in the simulation space. The local TGFβ concentration in nm cell −1 was calculated by dividing the sum of all TGFβ molecules located on the sites occupied by the cell by the estimated volume the cell would have had in three dimensions. This volume was calculated based on a diameter determined from either the cell's target area (for tumor cells) or the cell's actual area (for T cells), considering cells to be perfect spheres.

Intratumoral heterogeneity
We implemented intratumoral heterogeneity by randomly drawing all ODE model parameter values initially assigned to tumor cells from a Gaussian distribution with a mean of the default parameter value and a standard deviation (SD) of 1% of the default parameter value. We selected such a low SD because a greater SD triggered spontaneous EMT in a subset of tumor cells even in the absence of IFNγ. For parameters with either a positive or negative default value, drawing of parameter values from a Gaussian distribution for individual tumor cells could hypothetically lead to a negative or positive parameter value. Due to the relatively small SD, the probability of such events was only approximately 1.34 × e −2174 for any given parameter. In these rare cases, a value of 0 was used for the tumor cell. We additionally mimicked heterogeneity in IFNγ and TGFβ cell surface receptor expression by multiplying the sensed IFNγ and TGFβ concentrations by two additional parameters with a default value of 1 and heterogeneity as described above.  Table 1) were varied from the default model values (black) by +10/20% (red shades) and -10/20% (blue shades), while n BA parameters for these functions were varied from the default model values (black) by +1 (red) and -1 (blue).