Multiplexed imaging analysis of the tumor-immune microenvironment reveals predictors of outcome in triple-negative breast cancer

Triple-negative breast cancer, the poorest-prognosis breast cancer subtype, lacks clinically approved biomarkers for patient risk stratification and treatment management. Prior literature has shown that interrogation of the tumor-immune microenvironment may be a promising approach to fill these gaps. Recently developed high-dimensional tissue imaging technology, such as multiplexed ion beam imaging, provide spatial context to protein expression in the microenvironment, allowing in-depth characterization of cellular processes. We demonstrate that profiling the functional proteins involved in cell-to-cell interactions in the microenvironment can predict recurrence and overall survival. We highlight the immunological relevance of the immunoregulatory proteins PD-1, PD-L1, IDO, and Lag3 by tying interactions involving them to recurrence and survival. Multivariate analysis reveals that our methods provide additional prognostic information compared to clinical variables. In this work, we present a computational pipeline for the examination of the tumor-immune microenvironment using multiplexed ion beam imaging that produces interpretable results, and is generalizable to other cancer types. Patwa, Yamashita et al. utilize multiplexed imaging to demonstrate that profiling cell-to-cell interactions in the microenvironment can reveal predictors of recurrence and overall survival in triple-negative breast cancer, especially highlighting the relevance of immunoregulatory proteins. The authors also use multivariate analysis to provide additional prognostic information compared to clinical variables.

T riple-negative breast cancer (TNBC) is a subtype of breast cancer that is negative for estrogen receptor, progesterone receptor, and human epidermal growth factor receptor 2. Representing an estimated 10-20% of breast cancers, it is characterized by aggressive behavior, including earlier onset, larger tumor size, and a more advanced grade 1,2 . TNBC is the subtype of breast cancer with the poorest prognosis 3 , having a lower chance of survival 4,5 and higher risk of recurrence, especially within a short timeframe 6,7 . The absence of common breast cancer hormonal targets and high heterogeneity among TNBC tumors makes treatment management difficult, creating a need for more advanced interrogation of cellular processes within TNBC tumors 8 . Currently administered treatments, such as checkpoint inhibitors, only provide benefit to a small proportion of treated patients and are associated with high cost and toxicity 9 . Their effectiveness is limited, necessitating further interrogation of cancer-cell cues, factors in the tumor microenvironment, and host-related influences 10 . Currently, physicians are unable to separate patients with a low risk of recurrence from patients with a high risk of recurrence, making it difficult to deescalate treatments for those who may not need it and pursue more aggressive treatments for those who do 11,12 . To risk-stratify for overall survival, the American Joint Committee on Cancer staging system is the most commonly used technique in clinical practice; it is based on variables such as tumor size, nodal status, and the presence of distant metastasis. However, its survival estimates vary considerably because other prognostically relevant factors are excluded 13 . There is a need to identify additional biomarkers of TNBC to aid prognosis [14][15][16] . Identifying predictors of recurrence and survival in TNBC patients would allow improved patient stratification and targeted treatment plans, which would lead to better outcomes and spare patients from unnecessary aggressive therapies 17 .
The tumor-immune microenvironment (TIME) is a dynamic system comprising cancer cells, immune cells, and the surrounding extracellular matrix and vasculature 18 . The TIME is modulated by the expression and secretion of proteins that contribute to angiogenesis, immune suppression, and the coordination of the immune response 19 . Previous research has sought to discover the features of the TIME that are tumorpromoting or tumor-rejecting using transcriptomic and proteomic data [20][21][22][23] .
However, until recently, conventional histological techniques lacked the ability to measure the expression of a multitude of proteins at subcellular resolution while preserving spatial information 24,25 . Advancements in high-dimensional multiplexed imaging, such as multiplexed ion beam imaging (MIBI), have allowed for more direct interrogation of the TIME 26 while boosting standardization and reproducibility of results 27 . MIBI uses secondary ion mass spectrometry to image antibodies tagged with isotopically pure elemental reporters 28 . It is compatible with formalin-fixed paraffin-embedded (FFPE) tissue samples, the foremost preservation method of solid tissue in routine clinical pathology. MIBI enables in-depth analysis of the TIME, measuring the expression of more than 40 proteins simultaneously while preserving spatial information 29 and avoiding spectral overlap 30 and autofluorescence 31 .
This study builds on the work of Keren et al. 25 , who found structure in the composition and spatial organization of the TIME. TIME architecture was broadly classified as immune cold, mixed, or compartmentalized, based on the amount of immune infiltration into the tumor. The immune architecture was associated with patient survival.
However, previous research did not test the association between single-cell features of the TIME and clinical outcomes such as recurrence/survival. Although previous work has identified macro-level features associated with survival, there is still a need to study more granular features of the TIME at subcellular resolution, such as the expression patterns of individual proteins, which can add prognostically relevant information 32 and the characteristics of cell-to-cell interactions 33 .
In this work, we aim to uncover features of the TIME that are associated with recurrence and overall survival by analyzing MIBI scans of TNBC tissue 25,28 . The primary focus is to profile the proteins involved in cell-to-cell interactions and establish a link between the spatial organization of cells with varying expression patterns and clinical outcomes. We examine interactions involving functional proteins and immunoregulatory proteins in particular. As corollary aims, we demonstrate an association between protein co-expression patterns and recurrence/survival, examine proteins whose overall expression is associated with recurrence/ survival, and test associations between immune composition and recurrence/survival.

Results
Patient population. Our study examines 38 TNBC patients with no neoadjuvant treatment, a subset of the 41 patients examined by Keren et al. 25 . FFPE slides of breast tissue were taken from patients, scanned using MIBI, and subsequently segmented to demarcate cell boundaries 25 . Patient data regarding age, tumor grade, stage, cancer site, and clinical outcome-recurrence and overall survival (OS)-were also gathered (Table 1). We additionally gathered MIBI images of breast tissue of 8 healthy patients, a subset of the patients examined by Risom et al. 34 .
Dataset. MIBI scans produce images of protein expression from FFPE tissue, where each image has 44 channels; each channel conveys the expression of a certain marker on the tissue sample (Fig. 1a). Cellular segmentations for both TNBC and healthy patients' images were provided by Keren 25,34,35 . Cell type assignment for TNBC patients' images was also performed by Keren et al. 25 through a hierarchical methodology (Fig. 1b) The immune composition of the microenvironment is not associated with recurrence or survival. We examined whether the prevalence of certain cell populations in the TIME was associated with recurrence and survival. We measured the number of cells of each cell type in each patient and represented that number as a proportion of the total number of cells in that patient's sample. We then performed univariate Cox regression and performed a two-sided t-test of the variable coefficient to determine whether each cell type's prevalence was related to recurrence and overall survival. After performing Benjamini-Hochberg correction to account for multiple comparisons 36 , there were no cell types whose coefficients were significant for either recurrence (Table 2) or overall survival (Table 3).
Single-cell expression levels of functional proteins are not associated with recurrence or survival. We examined whether the expression of functional proteins in the cells of the tissue samples was associated with recurrence and survival (Fig. 2a). We calculated the per-pixel expression levels of each protein in each patient. The histograms of expression for several proteins are shown in Fig. 2b, and the histograms for all proteins are shown in Supplementary Fig. 1. For this analysis, we included only functional proteins, which stand in contrast to proteins used solely for lineage assignment; their expression is modulated according to the functional state of the cell. Proteins whose expression had been implicated by previous literature as having an important role Fig. 1 Overview of the computational pipeline. a Drawing of the layered structure of MIBI scans. Each MIBI image has dimensions of 2048 × 2048 pixels with 44 channels, where each channel represents expression for each protein; i.e., each pixel in the image at each channel conveys the concentration of that protein at that location. b Color-mapped image of cell segmentation performed on a MIBI image. The cell segmentation map has one channel with dimensions of 2048 × 2048. Each cell has its own cell type represented in colors referenced in the color bar. From these cell segmentation maps and the original MIBI images, we extract cell counts, measure protein expression, and quantify co-expression. c Voronoi tessellation diagram of the cell segmentation map. Each polygon corresponds to a cell in the original segmentation, such that each point in the area of the polygon is closer to the centroid of the corresponding cell than any other cell. Each polygon borders a finite number of other polygons, simulating adjacencies between cells. d Using Voronoi diagrams, we analyze interactions between neighboring cells. e An interaction matrix is computed for each patient, with the entry at row A and column B representing the number of times a cell positive for protein A was adjacent to a cell positive for protein B (top). The top half triangle of the matrix, split across the diagonal, is selected, as shown with the purple rectangles. These rectangles are then flattened to form one feature vector, i.e., interaction features, for each patient. f Interaction features are used to cluster patients, and the two patient clusters are compared with regard to recurrence/survival using Kaplan-Meier curves and the log-rank test.   Table 3). Keratin6 (coefficient = 0.025, HR = 1.025, p = 0.034) and HLA-DR (coefficient = −0.018, HR = 0.982, p = 0.045) were significantly associated with survival before correction. We placed Keratin6 and HLA-DR in a multivariate model to assess their relative prognostic relevance; Keratin6 remained significant (p = 0.04), whereas HLA-DR did not (p = 0.06). Expression of Keratin6 has been associated with poor survival outcomes in previous work 37 , which our finding loosely corroborates. CD45RO (coefficient = −0.019, HR = 0.981, p = 0.051) was nearly significantly associated with recurrence before correction. CD45RO has previously been discussed in the literature for its role in anti-tumor immunity, especially with regards to its expression in memory T cells 38,39 . Our findings loosely corroborate this, as CD45RO expression was associated with favorable recurrence outcomes.
Within this cohort, the expression levels of functional proteins did not hold reliable prognostic relevance. As such, we decided to move away from macro-level interrogation of the TIME, opting to add spatial context to our analysis by quantifying protein coexpression and cell-to-cell interactions.
Co-expression of functional proteins in patients' cells is associated with recurrence and survival. We sought to develop a computational pipeline to test the association between localized coordination of immune activity and recurrence/survival. We calculated the number of times that pairs of functional proteins were co-expressed across all cells of a patient, summarizing this information in a "co-expression matrix." (Fig. 2c).
The co-expression matrices provide information regarding the phenotypes of the cells present in each patient, placing the expression of proteins in a single-cell context. We used the coexpression information as features to describe each patient. Patients were grouped by hierarchical clustering, and the tree was cut to form two patient clusters (Fig. 2d). Our choice to select two clusters in this analysis, as well as all hierarchical clustering analyses, was motivated by silhouette score analysis 40 , which showed that division into two clusters would maximize intercluster dissimilarity (Supplementary Table 4). The recurrence/ survival outcomes of the two patient clusters were compared using two-sided log-rank tests. They diverged according to recurrence (χ 2 (df = 1, n = 38) = 3.75, p = 0.053), and survival (χ 2 (df = 1, n = 38) = 2.80, p = 0.094) (Fig. 2e). We also tested patient stratification when three clusters were chosen (Supplementary Fig. 2). The log-rank test (df = 2) p-value for recurrence was 0.093 and 0.222 for survival.
We assessed the relative importance of individual co-expression features using random forest variable importance. The four most important co-expression features were CD45RO + H3K27me3 (score = 0.822), CD45RO + H3K9ac (score = 0.767), CD45RO + HLA Class 1 (score = 0.646), and HLA-DR + IDO (score = 0.604). These results show that calculating the co-expression of proteins, namely the combinations listed above, can aid patient stratification. CD45RO's co-expression with HLA Class 1, an antigen used to promote cytotoxic T cell activation, is aligned with existing literature on melanoma 41

and may evidence coordination between memory T cells and cytotoxic T cells in cancer.
Cell-to-cell interactions contain prognostically relevant information. We examined cell-to-cell interactions by creating Voronoi tessellation diagrams out of the segmented MIBI images (Fig. 1c). Voronoi diagrams have been used previously to define spatial organization and cellular morphology 31,42 . Each cell's Voronoi polygon is created from the location of its centroid; its polygon will border some number of polygons from other cells 43 . These borders can be used to model cell-to-cell interactions (Fig. 1d); cells whose polygons share a border can be considered adjacent (Fig. 3a). Due to the geometry of the Voronoi tessellation algorithm, polygons will only border their immediate neighbors, which restricts the area of influence of a certain cell to the cells that are closest nearby.
We created an interaction matrix for each patient to describe the characteristics of the patient's cell-to-cell interactions by counting the number of times that specific pairs of proteins were involved in interactions. The entry in the matrix at row A and column B represents the number of times a cell positive for protein A was adjacent to a cell positive for protein B (Fig. 3b).
Our method of quantifying cell-to-cell interactions reveals that the spatial proximity of functional proteins contains valuable prognostic information; the proteins involved in interactions can be used as features to cluster patients into groups with significantly different outcomes.
By contrast, quantifying interactions involving lineage proteins does not hold prognostic relevance. Hierarchical clustering on features of lineage protein interactions did not result in clusters that differed in recurrence and survival outcome significantly ( Supplementary Fig. 3).
A drawing comparing the clusters formed from clustering on functional protein interaction features to the morphology distinction described by Keren et al. 25 is shown in Supplementary  Fig. 4.
Interactions involving immunoregulatory proteins predict recurrence and survival. We further examined a subset of functional proteins, the immunoregulatory proteins PD-1, PD-L1, IDO, and Lag3, which are in consideration as immunotherapy targets 9,44-48 . Prior research did not answer whether interactions involving these four proteins are associated with recurrence and survival, information that would be valuable in understanding their roles in TNBC progression.
To answer this question, we quantified spatial interactions between cells expressing immunoregulatory proteins, excluding all other proteins from the analysis (Fig. 4a). We reasoned that if interactions between cells positive for these proteins were associated with recurrence or survival, the result would point to the prognostic relevance of these proteins. Similar to the previous analysis, the counts of interactions were used as features to cluster patients (Fig. 4b). The Kaplan-Meier curves of the clusters formed from this analysis diverged significantly according to recurrence (χ 2 (df = 1, n = 38) = 7.60, p = 0.0058) (Fig. 4c). We also tested patient stratification when three clusters were chosen ( Supplementary Fig. 5). The three clusters diverged significantly with regards to recurrence (χ 2 (df = 1, n = 38) = 5.40, p = 0.020), demonstrating that the efficacy of risk-stratification was robust to the number of clusters chosen.
Ablation analyses reveal prognostically relevant groups of features. We further examined the cell-to-cell interaction data by performing multiple ablation analyses.
We also examined "homotypic" interactions-interactions involving the same protein. Homotypic interactions are found in the diagonal of the interaction matrix-they represent the number of times in a patient that a cell positive for protein A is adjacent to a cell positive for protein A (Fig. 4e). This information communicates the spatial proximity of cells with similar expression patterns. We used all of the homotypic interactions of functional proteins (the entire diagonal) as features for each patient and repeated the clustering analysis. The Kaplan-Meier curves diverged according to recurrence, χ 2 (df = 1, n = 38) = 3.43, p = 0.064, and diverged significantly according to survival, χ 2 (df = 1, n = 38) = 4.90, p = 0.027, indicating that the frequency of homotypic interactions is relevant information for survival prognosis.
We calculated the importance of interaction features by fitting a random forest model with interactions as predictors and cluster assignments as the response variable. Feature importance was scored using the mean decrease in Gini Index. The highestimportance feature was the Beta Catenin + CD45RO interaction feature (score = 0.794), followed by CD45RO + HLA-DR (score = 0.738), PD-1 + CD45RO (score = 0.716), PD-1 + H3k27me3 (score = 0.709), Lag3 + CD45RO (score = 0.706), IDO + PD-1 (score = 0.694), and Lag3 + PD-1 (score = 0.647). CD45RO was present in 4 of the 7 most important interactions, PD-1 was present in 4, and Lag3 was present in 2. These results point to interactions involving these proteins as being particularly useful for patient stratification; they contributed the most to clustering, and the resulting clusters differed significantly in terms of recurrence and survival.
Extracted features differ between healthy samples and TNBC samples. To confirm the validity of the features we extracted, we tested whether they differed between healthy tissue samples and TNBC tissue samples. The healthy tissue used in our analysis came from a different study, which profiled a different set of markers. There were 6 proteins common between the healthy images and TNBC images: FoxP3, IDO, Ki67, PD-1, PD-L1, and phospho-S6.
We also validated our method of profiling cell-to-cell interactions on the healthy tissue by subjecting it to our computational pipeline and testing whether the cell-to-cell interaction features of healthy tissue would be different from TNBC tissue. We reduced the interaction features for healthy and TNBC patients to two dimensions using uniform manifold approximation and projection (UMAP) 49 and plotted the reduced features for visualization. The resulting scatterplot showed separation between healthy and TNBC tissue ( Supplementary  Fig. 6b).
These results suggest that our computational pipeline succeeded in extracting tumor-specific single-cell spatial features that are prognostic for recurrence and overall survival. In addition, they demonstrate that our computational pipeline is applicable to a variety of MIBI datasets, as we applied the same methods to the two distinct datasets.
Multivariate analysis reveals features with independent prognostic relevance for recurrence and survival. To assess the prognostic importance of the features we identified, we fitted three multivariate Cox regression models, each of which included one of the cluster variables, two clinical variables (grade and age), and the immune architecture distinction described by Keren et al. 25 . We obtained coefficients and hazard ratios to determine whether the cluster variables added prognostic information.
Both of the clusters formed from cell-to-cell interaction features contained additional prognostic information for at least one clinical outcome. The immunoregulatory proteins interaction cluster contained independent prognostic information for recurrence (coefficient = −1.32, HR = 0.27, p = 0.02). The functional proteins interaction cluster contained independent prognostic information for survival (coefficient = −1.24, HR = 0.29, p = 0.04). These results suggest that our computational pipeline was able to extract additional prognostically relevant features and use them to risk-stratify patients.
Next, we assessed the relative prognostic relevance between each of the cluster variables. To do this, we fit random forests with six predictors: the three cluster variables, two clinical variables (tumor grade and age), and the immune architecture distinction. We then measured variable importance by calculating SHAP (Shapley additive explanations) values 50 and overall goodness-of-fit using Harrell's c-index 51 .
The random forest analysis corroborated our results from the multivariate Cox regression analysis. The immunoregulatory protein interactions cluster was the most relevant feature for recurrence (Fig. 5a), and the functional protein interactions cluster was the most relevant feature for survival (Fig. 5b). These features were more important than tumor grade, age, and tumor architecture. The c-index for the recurrence model was 0.718, and the c-index for the survival model was 0.731, indicating a good fit.

Discussion
TNBC is the most aggressive breast cancer subtype, with a higher risk of recurrence and lower probability of survival. It lacks clinically approved biomarkers for patient risk stratification, making treatment planning and management difficult. Previous research involving TNBC and multiplexed imaging did not analyze the prognostic relevance of protein co-expression patterns and cell-to-cell interactions. In this study, we aimed to examine the association between these features and recurrence/survival in TNBC patients by constructing a computational pipeline for the analysis of MIBI.
Our contributions are threefold. First, we identify possible predictors of recurrence and overall survival in TNBC, demonstrating that the information contained within cell-to-cell interactions and protein co-expression patterns can aid patient stratification and therapeutic design, as proven through evaluation of patient groups, statistical tests, and predictive modeling. Second, we demonstrate that the immune composition of the TIME does not always hold prognostic relevance, and should therefore be examined with caution. Third, we present a computational pipeline for the interrogation of the TIME that produces interpretable and conclusive results, making it potentially viable in a clinical setting.
Our primary focus was to examine cell-to-cell interactions in the TIME for the purpose of patient risk stratification and treatment management. Our findings showed that the type and number of cell-to-cell interactions involving functional proteins quantified by our pipeline were associated with both recurrence and survival in our cohort, and could possibly serve as a tool for prognosis. The two most important interaction pairs were CD45RO + Beta Catenin and CD45RO + HLA-DR, a finding that corroborates underlying biology. CD45RO marks memory T cells, which have been shown to mediate anti-tumor immunity 38,39 . Beta Catenin is expressed on tumor cells primarily 52,53 , so its interaction with CD45RO evidences the antitumor actions of memory T cells. HLA-DR is expressed on antigen-presenting cells 54 , so its interaction with CD45RO + cells evidence coordination between different immune cells to suppress tumor growth. The biology behind these interactions demonstrates that computational analysis of cell-to-cell interactions can elucidate immunological mechanisms playing a role in patients' tumors. The biological relevance of the proteins involved in these interactions could be further investigated through biological analysis of animals or clinical trials. Further, we found that "homotypic" interactions-interactions involving the same protein-hold predictive power. This finding indicates a coordinated immune response characterized by the localized enrichment of functional proteins.
Multivariate analysis revealed that the interactions of functional proteins contained independent prognostic information for survival, even when compared to clinical variables like tumor grade, age, and the immune architecture distinction determined by Keren et al. 25 , pointing to the potential efficacy of this technique for patient stratification and treatment management of TNBC. We also profiled the cell-to-cell interactions of healthy tissue. The interactive features of the healthy tissue were distinct from the interaction features of TNBC tissue in our cohort, indicating that profiling cell-to-cell interactions may also be able to convey information regarding the pathological state of the tissue. However, further research with larger sample size is necessary to confirm this.
We further analyzed four immunoregulatory proteins-IDO, Lag3, PD-1, and PD-L1-which are in consideration as immunotherapy targets 9,44-48 . We found that the expression profiles of these proteins were prognostically relevant, suggesting that these proteins can play a role in modulating tumor progression. A host of literature has described the importance of these proteins in TIME processes 9,32,44,55 , but only a small subset of such literature examines them in the context of paired cellular interactions. Interestingly, the individual expression levels of these proteins were not prognostically relevant; after Benjamini-Hochberg correction, none of the proteins had expression levels significantly associated with recurrence or survival. However, the cluster variable formed from their interactions was highly predictive of recurrence, as shown by multivariate analysis. Profiling the cellto-cell interactions involving immunoregulatory proteins revealed independent prognostic information when compared to tumor grade, age, and the immune architecture distinction.
Our methods differ from the previous analysis of these data in several ways. Keren et al. calculated interaction matrices by defining a distance of 39 micrometers to establish adjacent cells 25 ; however, the features within these interaction matrices did not result in patient clusters that differed significantly with respect to clinical outcome. This may suggest that using a set distance for adjacency is of insufficient spatial resolution to differentiate microenvironments. Our analysis also used a much lower threshold for cell protein positivity. This lower threshold may have improved the detection of important interactions. Voronoi diagrams and Delaunay triangulation have been used previously to define and examine cellular neighborhoods in colorectal cancer 31,56 . In contrast, we use Voronoi diagrams to examine protein expression in pairwise cellular interactions, rather than larger neighborhoods. We then use these pairwise cellular interactions to explain higher levels of abstraction, building interaction matrices to summarize patients' TIME overall.
We found that the co-expression profile of functional proteins in patients' cells was associated with recurrence and survival. The four most important co-expression pairs were CD45RO + H3K27me3, CD45RO + H3K9ac, CD45RO + HLA Class 1, and HLA-DR + IDO. These results point to highly specific cellular phenotypes, a trademark of a complex TIME 25,57,58 . Our computational pipeline presents an efficient, interpretable way to identify co-expression patterns and use them to risk-stratify patients.
Our methodology allowed for analysis of the cell types present in the TIME as a whole, providing a macro-level view of immune coordination. Our findings indicate that caution should be exercised when using the immune composition as a biomarker in clinical settings. After Benjamini-Hochberg adjustment, there were no cell types with significant prognostic value. This does not corroborate existing literature regarding the prognostic relevance of certain cell types, including the monocyte/neutrophil cell type 59 , the dendritic cell/monocyte cell type 60 , natural killer cells 61,62 , CD8+ T cells 27,63 , macrophages 64 , B cells 65 , CD4+ T cells 66,67 , and CD3+ T cells 68 .
The subcellular resolution achieved by MIBI allowed us to quantify the expression of individual molecules on a single-cell basis. After Benjamini-Hochberg correction, there were no proteins significantly associated with either recurrence or survival. Keratin6 and HLA-DR were associated with survival before adjustment, and Keratin6 remained significantly associated with survival when placed in a multivariate model with HLA-DR. This aligns with some existing literature. CD45RO was almost significantly associated with recurrence before adjustment, but its prognostic relevance was more clearly highlighted through its cell-to-cell interactions and co-expression patterns. This suggests that adding spatial context to the TIME can reveal otherwise hidden prognostically relevant information, a potential benefit of our developed computational pipeline for MIBI analysis.
A limitation of our work is that our results are derived from a sample of 38 TNBC patients that were treated at Stanford hospital from 2002 to 2015-further work is needed to validate these results on a larger cohort of patients. Although it was known that the patients had not undergone neoadjuvant treatment, further data regarding treatments pursued was not available; future research is necessary to examine associations across treatment types. In addition, this study was retrospective and performed with patients at a single institution. Our cell type classifications were found computationally, derived only from the expression of molecules that were a part of our chosen assay-future work should repeat this analysis using other biologically relevant molecules.
Nonetheless, this study presents a computational pipeline for the robust interrogation of multiple features of the TIME. We demonstrate the potential for cell-to-cell interactions and protein co-expression to improve prognosis and patient stratification. We found several statistically significant results within a limited cohort, suggesting that they may have large effect sizes and merit further exploration. Our methods produce interpretable results, which may make them beneficial in therapeutic design 69 . Furthermore, they can be applied to other cancer types, as they are generalizable to any MIBI scan.

Methods
Patient population and dataset. Our study examined 38 TNBC patients who were treated at Stanford Hospital from 2002 to 2015, a subset of the cohort examined by Keren et al. 25 . None of the 38 patients had undergone neoadjuvant treatment. Although the original cohort contained 41 TNBC patients, 3 of the patients were unusable for our analysis. Patients 22 and 38 lacked recurrence outcomes, and Patient 30's images were corrupted. These patients had no special type, with estrogen receptor and progesterone receptor positivity less than 1% and HER2 unamplified. 1 mm cores were taken from each patient and H&E stained. All samples were then stained with an antibody mix and scanned using MIBI-TOF. A computational pipeline converted the output of MIBI-TOF into images.
The dataset included two separate sets of 2048 × 2048 pixel images, representing a region of 800 2 square micrometers. The first set of images were 44-channel TIFFs that represent protein expression levels, where each patient had one TIFF. Each channel in the TIFF corresponded to one of the 44 molecules profiled in the study. Of the 44 molecules, 36 were biological macromolecules, such as double-stranded DNA or IDO, and 8 were elemental reporters. Each pixel in the image had a value representing the expression of the protein in that location. The second image set contained 38 grayscale segmentation of cells in the patient's sample. Patient data regarding age, tumor grade, stage, and recurrence, and survival outcomes were also gathered.
We additionally gathered 8 MIBI images of breast tissue of healthy patients. These MIBI images were originally collected for a different study, and they profiled a different set of markers. Six proteins overlapped between the TNBC samples and healthy samples: FoxP3, IDO, Ki67, PD-1, PD-L1, and phospho-S6.
The cellular segmentation of MIBI images of healthy patients was performed by Risom et al. 34 using DeepCell. Two distinct segmentations were performed. The first applied a three-pixel radial expansion and a stringent threshold for splitting cells. The second applied a one-pixel radial expansion and a lenient threshold. A post-processing step gave preference to the lenient threshold when the two segmentations were combined.
The authors complied with all ethical regulations involving human clinical data. Informed consent was obtained for all participants by previous studies. The study protocol was approved by the Stanford University Institutional Review Board.
Analysis of cell prevalence. We examined whether the cellular composition of the TIME was associated with recurrence and survival. We quantified the number of cells of each cell type in each patient. To isolate specific cell types at a time, we created binary masks of each grayscale value. Then, we found the number of connected components in each mask, which provided the number of cells of each cell type. After noticing a large variation in the total number of cells per patient, we divided each patient's cell type count by the total number of cells in their TIME to control for this lurking variable. Univariate Cox regression was then performed for each cell type to assess its association with recurrence and survival. Regression coefficients were examined using two-sided t-tests. P-values were adjusted using the Benjamini-Hochberg method 36 .
Single-cell protein expression. We examined whether the expression levels of functional proteins within patients' TIME were associated with recurrence and survival. For this analysis, we analyzed functional proteins, which modulate the activity of the cells in the TIME. These proteins stand in contrast to proteins used solely for lineage assignment; their expression is related to the functional state of the cell.
We labeled connected components and created a binary mask of each component to isolate the space taken up by each cell. We then applied this mask to the MIBI protein expression images, summing the value in each channel of the TIFF for each pixel in the mask. This created a 44-length vector of protein expression per cell. We realized that per-cell expression levels are dependent on cell size, so we divided the 44-length expression vector by the size (in pixels) of the cell. This left a 44-length vector representing the average per-pixel protein expression for a certain cell.
We calculated protein positivity thresholds from the expression levels of the image background, which lacks cells and therefore can act as a negative control. We calculated total protein expression in all background pixels in all patients and divided these values by the total number of background pixels across all patients (~67,000,000 pixels). We used each protein's threshold value to determine whether a cell was positive for a certain protein. Then, we counted the number of cells in each patient that were positive for each protein and divided the number by the total number of cells in the patient. The result was the proportion of cells in the patient that were positive for this protein. For example, 100% of a patient's cells would be positive for DNA, but only 30% might be positive for PD-1. Univariate Cox regression was used to determine the association between protein expression proportions and clinical outcomes. Regression coefficients were examined using two-sided t-tests. P-values were adjusted using the Benjamini-Hochberg method.
Functional protein co-expression. The co-expression of functional proteins in a single cell reveals functional status and immune coordination. We assessed the association between co-expression of functional proteins and recurrence and survival. We had previously measured single-cell protein expression and determined a threshold to designate cells as "positive" for each protein. We defined co-expression as an instance in which an individual cell is positive for a pair of proteins. For example, if a particular cell is positive for IDO, Lag3, and PD-1, it would have three instances of co-expression: IDO/Lag3, IDO/PD-1, and Lag3/PD-1. We constructed an 18 × 18 co-expression matrix for each patient to summarize the number of cells in the patient that co-expressed each pair of proteins. Because these matrices were symmetrical (a co-expression of IDO/Lag3 is the same as a co-expression of Lag3/ IDO), we divided the matrix in two across the diagonal and flattened the top half to create a feature vector for each patient. To control for lurking variables, we standard scaled features across patients. We then performed hierarchical clustering to segment patients according to these features.
Silhouette analysis 40 revealed that choosing two clusters would lead to the optimal segmentation, so we cut the dendrogram into two distinct clusters and compared the two groups using a two-sided log-rank test. Then, to assess the importance of individual co-expression features, we fit a random forest with all of the co-expression pairs as predictors and the cluster assignment as the response. We assessed variable importance using a mean decrease in the Gini index.
Voronoi tessellation. Analyzing cell-to-cell interactions requires a method of defining cell adjacencies. We used Voronoi tessellation diagrams to model cellular adjacencies within the TIME. Voronoi tessellation divides a planar space into a number of regions such that each point in the plane has its own region in the tessellation 43 . The sides of each Voronoi polygon are constructed to bisect two input points. Therefore, each line segment in the Voronoi tessellation represents the borders between two input points. Voronoi diagrams have been applied to single-cell imaging technology in the past, specifically for visualizing the spatial organization of colorectal cancer cells 31,42 . Due to the geometry of the Voronoi tessellation algorithm, polygons will only border their immediate neighbors.
We labeled connected components from the cell segmentation images to find each cell's centroid, which was then used to create Voronoi diagrams for each centroid. Therefore, every cell in the original cell segmentation images has a corresponding Voronoi diagram. We considered cells with bordering Voronoi regions to be adjacent and therefore interacting. This created a reliable foundation for upstream analysis.
Cell-to-cell interaction analysis. We used the borders created by Voronoi diagrams to iterate over all cell adjacencies in each MIBI image. Each adjacency represented an individual interaction between two cells. We constructed two lists: List 1 contained the names of the proteins that Cell 1 was positive for, and List 2 contained the names of the proteins that Cell 2 was positive for. We took the Cartesian product of the two lists to find all of the combinations of proteins present in this interaction. For example, if Cell 1 was positive for PD-L1 and Lag3, and Cell 2 was positive for PD-1 and IDO, then we would count the following: PD-L1 + PD-1, PD-L1 + IDO, Lag3 + PD-1, and Lag3 + IDO. These pairs would be tallied in the overall interaction matrix for each patient, in which the value at row A and column B represents the number of times a cell positive for protein A was adjacent to a cell positive for protein B. For this analysis, we only counted interactions between functional proteins, excluding proteins used for lineage assignment.
We selected the top half of the symmetric matrix and flattened it to create feature vectors for each patient. Hierarchical clustering was performed and the dendrogram was cut to produce two clusters based on silhouette score analysis. These two clusters were compared using Kaplan-Meier curves, a two-sided log-rank test, and Cox regression. To assess the importance of individual interactions, we fit a random forest with all interactions as predictors and the cluster assignment as the response. We measured variable importance using the mean decrease in the Gini index.
Healthy tissue analysis. We applied our computational pipeline to a set of 8 MIBI images of healthy tissue to validate our methods. We examined 6 proteins: FoxP3, IDO, Ki67, PD-1, PD-L1, and phospho-S6. We calculated the expression levels of the proteins in an identical fashion to the previous analysis of the TNBC images. We performed a two-sided Wilcoxon rank-sum test to compare expression levels between TNBC tissue and healthy tissue.
We also profiled cell-to-cell interactions in the healthy tissue in an identical manner to previous analysis. After obtaining interaction features, we performed dimensionality reduction using UMAP 49 . The resulting reduced features were compared to UMAP-reduced features of TNBC images.
Multivariate analysis. To assess whether the features identified by our computational pipeline contained independent prognostic information, we performed multivariate Cox regression. We fit three Cox Proportional Hazard models, each of which contained one of the three cluster variables identified by our study (protein co-expression, functional protein interactions, immunoregulatory protein interactions), two clinical variables (tumor grade and age), and the immune architecture distinction determined by Keren et al. 25 . We found the hazard ratio of each cluster variable and hypothesis tested the coefficient of each cluster variable to determine whether the variables contained additional prognostic information.
We also fit random forests to measure relative variable importance. We included all six variables in the random forests and calculated SHAP 50 values to get stable estimates of variable importance. SHAP values quantify the change in the model prediction that would result from conditioning on a certain feature. They have been shown to be more aligned with human intuition regarding feature importance and attribution. We measured overall goodness-of-fit using Harrell's c-index.
Pseudocode explaining the steps of statistical analyses is shown in Supplementary Fig. 7. A runnable Jupyter notebook with the specific code used for experiments can be found at github.com/aalokpatwa/rasp-mibi/blob/main/ rasp_mibi_pipeline.ipynb.