In the tumor microenvironment (TME), various types of immune cells interact with each other and with cancer cells, playing critical roles in cancer progression and treatment [1]. Numerous studies have reported that the infiltration levels of specific immune cells are associated with patient prognosis and response to immunotherapies [2, 3]. For instance, the density of pre-existing tumor infiltrating lymphocytes in the TME has been found to positively correlate with patient responses to anti-PD-1 treatment in triple-negative breast cancer (TNBC) [3]. However, the relationship between immune cell-cell interactions (CCIs) in the TME and patient clinical outcomes remains unclear due to the limited availability of large-scale datasets for systematic CCI investigation. Recently, imaging mass cytometry (IMC) has been utilized to characterize the immune landscape within the TME of tumor samples [4]. IMC can detect 30 to 40 protein markers on a single tissue slide, enabling the visualization of spatial distributions of various cell types at single-cell resolution. Analyzing IMC data enables the quantification of interactions between all TME cell types by examining their spatial distributions. In this study, we conducted a systematic analysis to investigate the association between CCIs and treatment responses of cancer patients using a large IMC dataset. The dataset comprises the immune landscape of 660 tumor samples from 279 TNBC patients enrolled in a randomized clinical trial [4]. These patients were treated with either neoadjuvant chemotherapy (n = 141) or immunochemotherapy therapy (chemotherapy combined with anti-PD-L1 immunotherapy, (n = 138), with tumor samples collected at three time points for IMC analysis: baseline, early on-treatment, and post-treatment. We applied a modified method introduced by Windhager et al. [5] to quantify interactions for all pairs of cell types captured by IMC and examined their associations with patient outcomes (Supplementary Material and Methods). Our results indicated that compared to the infiltration levels of immune cells, CCIs between specific immune cell types were more strongly correlated with patient responses. Notably, we found that the interaction between regulatory T cells (Treg) and GZMB+ cytotoxic CD8+ T (Tc) cells in pre-treatment samples was predictive of patient response to immunochemotherapy but not to chemotherapy alone in TNBC. The processed IMC data provides the coordinates of all single cells along with their cell type annotations. To quantify the interaction from cell type X to Y (X→Y), we calculated the average number of X cells among the 10 nearest neighbors of each Y cell and standardized this as a Z-score by comparing it with a null distribution generated through permutations. In each permutation, we shuffled the labels of all cell types except epithelial cells (Figure 1A). We applied this method to the TNBC IMC dataset, which included 20 non-epithelial cell types: endothelial cells, fibroblasts, myofibroblasts, PDPN+ stromal cells, CA9+ cells, Treg cells, CD4+PD1+ T cells, CD4+TCF1+ T cells, CD8+TCF1+ T cells, CD8+ T cells, CD8+PD1+ exhausted T cells (CD8+PD1+T_Ex), CD8+GZMB+T cells, CD79a+ plasma cells, CD20+ B cells, CD56+ NK cells, PD-L1+ antigen-presenting cells (APCs), PD-L1+IDO+APCs, dendritic cells (DCs), M2 macrophages (M2Mac), and neutrophils. In total, for each IMC image, we calculated Z-scores for 380 between-cell interactions and 20 self-interactions (Figure 1B). In the TME all cell types tend to cluster with their own type, with fibroblasts, stromal cells and endothelial cells displaying stronger clustering tendencies compared to immune cell types. Additionally, several between-cell interactions are prevalent in baseline samples, such as the interactions of CD56+NK cells and CA9+ cells. The baseline interaction between Treg and Tc (CD8+GZMB+T) is associated with patient response to immunochemotherapy in triple-negative breast cancer. (A) Schematic diagram illustrating how the interaction between two cell types (X and Y) were quantified by using a CCI Z-score. For X→Y, the average number of X neighbors of all Y cells in an image was calculated and transformed into a Z-score based on a null distribution estimated from permutations. In each permutation, the positions of all epithelial cells were persevered, while the positions of other cells were shuffled. Created with BioRender.com. (B) Interactions between 20 non-epithelial cell types across all baseline tumor samples. Dot size represents the fraction of samples with significant Z-scores for a cell pair (P < 0.001). Cell pairs with significant interactions in over 20% of samples are highlighted with a circle. Color represents the average Z-score for a cell pair across all samples. (C) Volcano plot illustrating CCIs with significantly different Z-scores between responders and non-responders in the immunochemotherapy arm. (D) Comparison of Treg and Tc (CD8+GZMB+T) cell interactions between responders and non-responders at baseline, prior to immunochemotherapy treatment. (E) Volcano plot showing that immune cell fractions are less informative in distinguishing responders from non-responders. (F) ROC curves comparing the classification accuracy of responder status between CCI Z-scores (CD8+GZMB+T→Treg and Treg→CD8+GZMB+T) and the fractions of CD8+PD1+T_Ex and PD-L1+APCs cells. The values in the parentheses are AUC scores. (G) Temporal trends in Treg→CD8+GZMB+T Z-scores from baseline to on-treatment and post-treatment in responders, with lines connecting samples from the same patient collected at the three time points. (H) Changes in the fractions of Treg and CD8+GZMB+T in responders and non-responders during immunochemotherapy. (I) Volcano plot illustrating CCIs with significantly different Z-scores between responders and non-responders in the chemotherapy arm. (J) Baseline comparison of CD79a+ Plasma and CA9+ cell interactions in responders and non-responders, prior to chemotherapy. (K) Schematic diagram summarizing the link between baseline Treg-Tc interaction and patient response to immunochemotherapy. Created with BioRender.com. Note: Panels (C-H) pertain to the immunochemotherapy arm, and Panels (I-J) pertain to the chemotherapy arm. Abbreviations: NR, non-responder (RD, Residual Disease); On, On-treatment; Post, Post-treatment; R, Responder (pCR, pathological Complete Response). Next, we compared the CCI Z-scores between responders and non-responders to identify baseline CCIs associated with patient response to immunochemotherapy. At a significance level of P < 0.01, we identified 13 significant CCIs (Figure 1C and Supplementary Table S1). For example, the Z-scores for both CD8+GZMB+T→Treg (P < 0.001) and Treg→CD8+GZMB+T (P < 0.001) were significantly higher in responders than in non-responders at baseline (Figure 1D). This finding suggests that patients sensitive to neoadjuvant immunochemotherapy tend to exhibit stronger interactions (i.e., spatially closer) between these two cell types in the pre-treatment samples. CD8+GZMB+T cells are activated Tc cells, which serve as primary killers of tumor cells and are correlated with improved treatment outcomes [6-8]. Conversely, Tregs are suppressive immune cells that inhibit the activation and function of effector T cells, including Tc cells [9, 10]. Another example is the baseline interaction between PD-L1+IDO+APCs and CD56+NK cells, which is negatively correlated with patient response (Supplementary Figure S1). Among the 13 significant CCIs, CD8+PD1+T_Ex→CD8+PD1+T_Ex is the only self-interaction, showing significantly higher Z-scores in responders compared to non-responders (Supplementary Figure S2). In contrast, when the fractions of immune cells among all TME cells were analyzed, only PD-L1+IDO+APCs (P = 0.002) were significantly associated with patient response to immunochemotherapy, followed by CD8+GZMB+ T cells (P = 0.031) (Figure 1E and Supplementary Table S2). These results indicate that the spatial interactions between specific cell types offer a more effective set of features for predicting patient response to therapeutic treatment than the infiltration levels of immune cells. More specifically, when the Z-scores for CD8+GZMB+T→Treg and Treg→CD8+GZMB+T were used to classify responders and non-responders, they achieved area under the curves (AUC) of 0.717 and 0.706, respectively (Figure 1F). These values significantly outperformed the immune fractions of the two most predictive immune cell types, PD-L1+IDO+APCs (AUC = 0.651) and CD8+GZMB+T cells (AUC = 0.613) (Figure 1F). Wang et al. [4] reported that the fraction of proliferating CD8+TCF+ T cells (high Ki67 expression) was highly correlated with patient response, yielding an AUC of 0.622 according to our computation (Supplementary Figure S3). When patients were stratified into two groups based on the median CCI Z-score for CD8+GZMB+T→Treg, the group with strong interactions demonstrated a response rate of 66.1%, in contrast to only 25.8% in the group with weak interactions (P < 0.001, Supplementary Figure S4). As the interaction between CD8+GZMB+T cells and Tregs most effectively classified responders versus non-responders among TNBC patients, we examined how this interaction changes during immunochemotherapy. In the TNBC immunochemotherapy arm, 75 patients had samples collected at all three timepoints. Using IMC data for these patients, we compared the Z-scores for CD8+GZMB+T→Treg and Treg→CD8+GZMB+T interactions across baseline, on-treatment, and post-treatment samples. Interestingly, the Z-scores showed a significant decrease from baseline to on-treatment (P = 0.002, paired Wilcoxon test) and further to post-treatment (P < 0.001, paired Wilcoxon test) in responders (Figure 1G). In contrast, no significant changes were observed in non-responders (Figure 1G). During immunochemotherapy, both responders and non-responders showed a significant increase in Tregs infiltration levels during treatment, followed by a post-treatment decrease, as measured by their fractions among all TME cells (Figure 1H). However, CD8+GZMB+T cells exhibited a significant increase during treatment and a decrease in post-treatment only in responders, with no significant changes in non-responders (Figure 1H). To examine the influence of neighboring Tregs, we compared baseline marker protein levels, measured by IMC, in CD8+GZMB+T cells with at least one Treg among their top 10 nearest neighbors (Treg-proximal) to those without Treg neighbors (Treg-distal). We observed significantly higher PD-L1 and PD-1, but lower expression of GZMB, in Treg-proximal CD8+GZMB+T cells compared to Treg-distal cells at baseline (Supplementary Figure S5). These findings suggest that Tc cell function is suppressed when they are in close proximity to Tregs. Using the same method, we next investigated the chemotherapy arm and identified baseline CCIs associated with the response of patients treated with neoadjuvant chemotherapy alone. In total, we identified 6 response-associated CCIs at a significance level of P < 0.01 (Figure 1I and Supplementary Table S3). For example, patients with weaker interactions between CD79a+ plasma cells and CA9+ cells were more likely to respond to chemotherapy (Figure 1J). None of these six response-associated CCIs identified in the chemotherapy arm were significant in the immunochemotherapy arm. However, when comparing the t-statistics (responder vs. non-responder) for all 400 CCIs between the two arms, we observed a weak but significant correlation (R = 0.27, P < 0.001), suggesting that some CCIs may have similar effects on patient response to both treatments (Supplementary Figure S6). Taken together, our analyses of the TNBC data indicate that the baseline interaction between Treg and Tc cells (CD8+GZMB+T) plays a critical role in patient response to immunochemotherapy but not to chemotherapy alone (Figure 1K). In this study, we modified a commonly used method to quantify the strength of spatial CCIs for all possible combinations of TME cell types captured in each IMC image. The resultant CCI Z-scores were used as features to correlate with the immunochemotherapy responses of TNBC patients. Compared to conventional features based on immune cell fractions, these spatial interaction features showed a stronger correlation with patient treatment responses. Specifically, TNBC patients with stronger interactions between Tregs and Tc cells in their pretreatment samples were more likely to respond to neoadjuvant immunochemotherapy (anti-PD-L1 combined with chemotherapy) but not to neoadjuvant chemotherapy alone. Consistently, these interactions were significantly reduced during treatment, from baseline to early on-treatment, and further reduced post-treatment in responders, but not in non-responders. In the presence of neighboring Tregs, Tc cells exhibited significantly higher levels of PD-L1 and PD-1, but lower levels of GZMB, suggesting a molecular impact of their spatial interactions. These results highlight the critical roles of CCIs in cancer progression and treatment. Interestingly, a high ratio of Treg to CD8+ T cells (Treg/CD8+ T) has been reported to be associated with poor patient outcomes in gastric cancer when treated with immune checkpoint blockade therapy [10]. Although here we identified the Treg-Tc interaction from IMC data, it can also be measured by conventional imaging techniques such as immunohistochemistry, making it readily translatable into clinical applications. It should be noted that a limitation of this study is the lack of independent IMC datasets for further validation. Additionally, due to the limited availability of molecular features, it is challenging to determine the molecular mechanisms underlying the link between Treg-Tc interactions and immune responses. Further validation and mechanistic investigation in TNBC and other cancer types will be critical when more data become available. With the increasing application of IMC and other single-cell imaging and transcriptomic technologies, we anticipate incorporating CCI metrics into prognostic and predictive models for more accurate clinical outcome predictions in cancers. Chao Cheng designed the concept and the method. Chao Cheng, Xiang Wang, Jing Dong, and Jian-Rong Li wrote the main manuscript text and prepared the figures. Yupei Lin and Bikram Sahoo contributed to the preparation of figures and the manuscript. Yanhong Liu, Yong Li, Robert Taylor Ripley, Jia Wu, Jianjun Zhang, and Christopher I Amos revised the manuscript and figures. All authors reviewed the manuscript. Not applicable. The authors declare no competing financial interests. This study is supported by the Cancer Prevention Research Institute of Texas (CPRIT) (RR180061) and the National Cancer Institute of the National Institute of Health (1R01CA269764). Not applicable. The function for calculating the Z-scores of CCIs based on IMC data and example code for its implementation is included in GitHub https://github.com/wangsky137/Cancer_IMC_CCI. The original TNBC dataset utilized in this study are available for download as referenced in Wang et al. [4]. Please note: The publisher is not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing content) should be directed to the corresponding author for the article.
