Metabolic characterization of tumor-immune interactions by multiplexed immunofluorescence reveals spatial mechanisms of immunotherapy response in non-small cell lung carcinoma (NSCLC)

Domain: neuroscience
Published: 2026-02-03
Editorial outcome: major_revision

Abstract
Immune checkpoint inhibitors (ICI) have improved clinical outcomes for some patients with advanced NSCLC, however a substantial proportion of patients remain treatment resistant. Here we analyze the NSCLC tumor micro-environment (TME) using multiplexed immunofluorescence (mIF) of biopsies taken from patients prior to ICI treatment. We apply a deep-learning model to classify the cellular phenotypes and probe functional and metabolic states of both tumor and immune cells, aiming to reveal predictive features of response to ICI. Tissue neighborhoods are generated to allow geometric profiling of spatial densities and interactions at a range of scales. Multivariate modelling of ICI response yields a model that predicts progression-free survival (PFS) over 24 months (AUC = 0.8). The selected features in the model imply a role for cell-cell proximities within discrete metabolic contexts. These tissue insights may supplement our understanding of the current paradigms around classical immunology in the NSCLC TME and its influence on immunotherapy outcomes.

Introduction
Underlying mechanisms of resistance to immune checkpoint inhibitors (ICI) are poorly understood, despite representing most clinical outcomes for patients with non-small cell lung cancer (NSCLC) (70–80%)1. Treatment resistance, either primary or acquired, to standard of care therapies poses a major barrier, making lung cancer the predominant cause of cancer related deaths worldwide2. Primary diagnosis at advanced and stage IV disease occurs in ~40% of cases, with a 5-year survival rate of less than 9%3. While surgical resection, chemotherapy, radiotherapy, and ICI are associated with benefit at early disease stages, treatment at advanced stages is less effective, in part due to ICI therapy resistance. Current generations of ICI therapeutics, including anti-programmed cell death 1/ligand (PD-1/PD-L1) and anti-cytotoxic T-lymphocyte associated protein 4 (CTLA-4)4, aim to block T-cell inhibitory signaling and tumor directed self-tolerance. The mechanisms of refractory disease are varied and include: low neoantigen generation, inhibition of antigen presentation and suppressed recruitment of antigen presenting cells, disrupted T cell migration, T cell exhaustion, a hypoxic tumor microenvironment (TME), impairment of T cell activation, and PD-L1 splice variants bypassing PD-1 blockade5,6. While immunology informs on the roles that cells play within the TME, the spatial cellular arrangement of tumors remains poorly defined. Recent studies indicate the wealth of predictive power that the cellular composition and the derived cellular niches have with treatment outcomes7, expounding the utility of spatial analysis for predictive modeling of successful individualized treatments.

Spatially profiling tumors to identify functionally relevant neighborhoods opens the potential for prognostic biomarker guided treatment. The integrated human lung cell atlas (HLCA)8 has provided a detailed stratification of the human lung, albeit the single-cell components, not a spatial blueprint. Of particular interest is the categorization of distinct markers from primary NSCLC that are associated with disease progression and can be used to better inform treatment prospects. At steady state, lung tissue architecture is complex, at tumor state the disrupted milieu proves fatal, generating a highly inflammatory lung injury pattern9. The heterogeneity of tumor, stroma, and immune cells within the TME and between individuals makes a one size fits all approach to treatment unsuitable.

In this work, we sought to combine robust computational methods with meaningful statistical reasoning to understand the biologically relevant measurements from a custom-designed mIF panel. In doing so we develop an analytical pipeline for the evaluation of putative predictive biomarkers for ICI treatment that can be tailored and applied to diverse cohorts. We demonstrate the predictive power of a multivariable model dominated by metabolic and interaction features, which yields high accuracy in predicting progression free survival over 24 months.

Results
Patient characteristics and cohort selection
Two tissue microarrays (TMAs) were constructed from independent cohorts from NSCLC patient tumors presenting clinically between 2011–2017 (YTMA404) and 2017–2019 (YTMA471). TMA slides were run in parallel and subjected to multiplex (44-plex) tissue staining (Fig. 1a). Cohort metadata prior to sample refinement consisted of eighty-two patients with available clinical annotations shown (Fig. 1b). RECIST best overall response (BOR) information was available for seventy-seven patients, consisting of twenty-nine with progressive disease (PD), twenty-six with stable disease (SD), seventeen with partial response (PR), and five with complete response (CR). Samples from non-primary lung tissues were removed (n=10 lymph node, n=2 bowel, n=4 brain, n=1 spinal, n=1 skin), and only specimens taken pre-ICI treatment were considered for analysis (Fig. 1d). The final cohorts consisted of fifty-five patients with twenty-five females and thirty males. Smoking status varied with five being current smokers, forty-five former smokers, and five never smokers. The sampled histologies comprised 36 adenocarcinomas, 17 squamous cell carcinomas, and 2 adenosqu2 adenosquamous carcinomas, with varying disease stages at time of diagnosis between stage II (n=5), stage III (n=6), and stage IV (n=44). ICI treatment was given in the advanced stage setting. The mean patient age was 67.1 years (±10.3 years) (Fig. 1b). Either objective response (partial or complete) or progression free survival above 6 months was used as the primary endpoint to group patients who received clinical benefit beyond 6 months post treatment (CB6=yes, n=46) vs those who progressed within 6 months (CB6=no, n=34) (Fig. 1e). This definition allowed both the clinically meaningful grouping of advanced stage disease management as well as overcoming the ambiguity in assessing practical responses for those patients with RECIST stable disease.

Cell phenotyping
Representative images for markers indicative for immune lineage (CD45, CD3e, CD8, CD4, CD20, CD14, CD68, CD11b, Vimentin, CD31, CD34), immune status (CD44, CD45RO, FoxP3, Granzyme B, HLA-A, LAG3, ICOS, IDO1, PD-L1, PD-1), dividing cells (Ki67), epithelial cells (E-cadherin, PanCK) and metabolic status (ASCT2, ATPA5, citrate synthase, CPT1A, G6PD, GLUT1, hexokinase 1, IDH2, NA/K ATPase, pNRF2, SDHA) are shown (Fig. 1c). A deep learning cell phenotyping platform developed by Nucleai was used to classify cells (Fig. 1f) according to a rules table shown diagrammatically (Supplementary Fig. 1a) Final defined cell lineages included tumor cells (41.4%), macrophages (14%), myeloid NOS (not otherwise specified) (1.3%), immune NOS (3%), granulocytes (2.5%), B Cells (1.6%), plasma cells (1.7%), CD4+ Tregs (1.9%), CD4+ T cells (8.6%), and CD8+ T cells (4.5%). Cell lineages were then assessed for functional phenotypes (Supplementary Fig. 1b). Overall expression of functional markers are shown per cell type, where metabolic expression patterns predominated in both immune and tumor cells (Fig. 1g).

Cell type proportions
The gross cellular makeup of patient tissues was first inspected for cellular frequencies and functional expression patterns. Tumor cell composition ranged from 10% to 90% within the cohorts and isolation of non-tumor cells revealed macrophages and CD4+ T cells to be the predominant immune cells, with significant numbers of fibroblasts and myofibroblasts contributing to the tissue architecture (Fig. 2a). Overall abundance of the functional status of cell types was measured (Fig. 2b), indicating tumor cell expression patterns for functional markers (IDO1, PD-L1, HLA-A, Ki67 and vimentin). Of immune cells, CD4+ T cells indicated positivity for all immune functional markers, while CD4+ T regs expressed patterns of PD-1 and ICOS positivity. CD8+ T cells demonstrated expression of granzyme B and PD-1, while macrophages showed positivity for all functional markers. Metabolic marker positivity among both tumor and immune cells indicated expression of markers consistent with mitochondrial respiration (oxidative phosphorylation: OXPHOS). Interestingly, both macrophages and tumor cells had high positivity counts for ATPA5, citrate synthase, Glut1, IDH2, and SDHA (Fig. 2b) suggesting high OXPHOS activity and distinct activated metabolic phenotypes in both tumor and innate immune cells.

To statistically evaluate these cell type positivity states at the gross tissue level for clinical outcome, we applied Mann–Whitney U tests. Fold-change enrichments were assessed within patient outcome groups: clinical benefit (green) and no clinical benefit (purple). Elevated Na/K ATPase+/IDO1+ immune NOS cells, CPT1a+ endothelial and myofibroblast cells, G6PD HIGH B cells, Na/K ATPase+ macrophage cells, were enriched in patients who received clinical benefit, consistent with an immune environment with a high energy demand. Conversely, Granzyme B+ plasma cells, granzyme B+ macrophages, G6PD+ fibroblast cells, and hexokinase 1 HIGH tumor cells were elevated in non-clinical benefit patients. Tumor cells positive for amino acid metabolism (ASCT2) and TCA enzymes (SDHA), commonly overexpressed in NSCLC10 and associated with poor prognosis, were elevated in patients who relapsed within 6 months of ICI treatment (Fig. 2c, d).

Examination of all marker positivity features after correcting for multiple testing revealed that only granzyme B+ macrophages were significantly associated with poorer outcomes after ICI treatment (Fig. 2e, f). Assessment of marker positivity for progression-free survival supported the role or the presence of cytotoxic granzyme B+ macrophages in treatment resistance, where higher levels were concomitantly associated with poorer PFS outcomes (Fig. 2g).

Tissue region dissection
Cellular neighborhoods (CNs)11 were used to construct unbiased tissue regions. K-means clustering on the KNN histograms of binary cell types (tumor or non-tumor) allowed tumor/stroma CNs (K=2) or tumor/stroma/interface (K=3) CNs to be generated (Fig. 3a, b). The K=2 CNs were designed to capture tumor/stromal regions, whereas the K=3 CNs were designed to provide resolution of the interface region between tumor cells and their surrounds, where functional immune cells and metabolic gradients are expected to play a critical role in the clinical response. Cells were assigned within these regions for assessment of spatially resolved frequencies. Overall, within tumor regions, increased plasma cells exhibited a trend towards clinical benefit, while higher granulocyte and CD4+ Treg infiltration appeared to trend negatively with response (Fig. 3c).

Marker positivity for cell types was assessed within our assigned tissue regions (Fig. 3d), where metabolically active B cells in tumor areas associated with worse clinical outcomes. B cell expression of GLUT1 HIGH and G6PD/hexokinase 1 coupled with Ki67 positivity indicated a highly glycolytic phenotype with concomitant activation of the pentose phosphate pathway known to be active in proliferating cells12. Interestingly, additional levels of citrate synthase HIGH, IDH2, and SDHA positivity suggested elevated levels of OXPHOS, which coupled with a proliferative phenotype, point to a role for active, expanding B cell populations within tumor regions in ICI refractory disease. The role for activated B cells was further evidenced by elevated levels of IDO1 in tumor regions13. We observed granulocyte expression of G6PD and GLUT1 in tumor regions to be associated with worse clinical outcome and localized the negatively associated granzyme B+ macrophages within the tumor compartment. Conversely, plasma cell expression of OXPHOS enzymes ATPA5 and GLUT1 within tumor regions was beneficial, suggesting an activated state that may influence localized antibody secretion and recognition within the tumor, priming the tissue for ICI response14.

Within the interface region, fibroblast cells expressing PPP pathway marker G6PD HIGH were enriched in non-responding patients, while levels of plasma cells expressing granzyme B and hexokinase 1, and B cells expressing PD-1 and SDHA also associated with poorer outcomes. Interestingly, stromal expression of tryptophan pathway enzyme, IDO1, in granulocytes, immune NOS cells, and macrophages as well as HLA-A in tumor cells associated with better outcomes, suggesting a benefit for lower inflammatory stromal microenvironment (Fig. 3d). Assessment of marker positivity levels for PFS indicated tumor localized granzyme B+ macrophages to be enriched in patients with early relapse, passing FDR adjustment (Fig. 3e) in accordance with our previous analysis.

Immune cell ratios are widely discussed in the context of anti-tumor responses15,16. We explored this within both definitions of tumor regions (K=2, K=3 CNs) and observed that higher ratios of granulocytes to CD8+ T cells within tumor regions associated with relapse, as well as high stromal B cells to tumor region myofibroblast and immune cell populations. When incorporating the tumor interface region, we observed that elevated levels of granulocytes in the interface region relative to stromal B cells associated with poorer outcome (Supplementary Fig. 2).

Metabolic neighborhoods
To investigate the role of metabolic activity in each TMA core, a metabolic density clustering approach derived from CNs11 (Fig. 4a) was applied to identify metabolic neighborhoods (MBNs; see methods). Four distinct clusters were described: minimal, low, medium and regulatory activity, and high metabolic activity. These neighborhood descriptions specified tissue regions which followed discrete levels of metabolic pathway activity with high MBN displaying high SDHA expression and low MBN displayed reduced CPT1a expression, concurrent with high and low energy demands. Overall characterization of cell types for their predominant pathways indicated that granulocytes highly expressed several OXPHOS pathways, while tumor cells had high ATP synthesis (Fig. 4b, c).

Examination of the frequency of each discrete metabolic pathway within our previously defined tumor/stroma/interface tissue regions indicated that OXPHOS pathways clustered together and were enriched within the tumor and interface regions, consistent with a predominant metabolic contribution by tumor cells. Interestingly, RECIST BOR responder groups (PR, CR) demonstrated lower levels of cellular energy production pathways in tumor and interface regions when compared to other BOR groups, consistent with the role for the Warburg effect and higher glycolytic activity in more aggressive cancers (Fig. 4d).

We probed changes in cell frequency within each MBN for associations with clinical benefit. We noted an increase in tumor cells within higher activity MBNs was associated with poorer outcomes, as well as the presence of granulocytes in lower activity MBNs. Interestingly, enrichment of immune NOS cells in high MBN regions appeared to be associated with clinical benefit (Fig. 4e), indicating beneficial lymphocyte infiltration to those regions. High metabolic activity was dictated by ATP synthesis, glycolysis, TCA Cyle, regulatory, and amino acid uptake pathways. Medium and regulatory activity was associated with increased FA oxidation and the PPP (Fig. 4f).

Higher level feature engineering to better model tissue architecture
Our analysis thus far focused on cellular frequency within tissue regions defined by tumor or metabolic descriptive regions. We sought to extend the depth of analysis by designing a comprehensive feature generation pipeline that combined regional properties with cell-cell interaction and proximity features. Features were engineered to describe density, proximity, and relative clustering of base cell types, functionality types, and metabolic types at a variety of distance scales within the core. Cell proportions and proportion ratios were specified for all cell types, including immune and tumor functionality types, within K=2 and K=3 CNs, in addition to the four MBNs. Large-scale interaction and proximity were characterized by the JSD Score feature, which provides asymmetric and normalized measure of the overlap of the cell type probability distributions. A Python implementation of the distance as previously described17 was computed for the cell/cell plus functional/metabolic density within each neighborhood. More localized measures of the proximity of cells were captured by the G-Cross feature, which represents the area under the curve (AUC) of the nearest neighbor cumulative radial distribution of cell/cell or functional/metabolic types, computed up to a radius of 150 µm. Since this feature captures nearest-neighbor distances, this radius was chosen to be above the scale of expected local interactions, but less than the core radius to reduce biases by missing cells outside of the core boundary. Localized boundary packing of cell types has been found to correlate with patient response to therapies in spatial biology18,19, therefore we built an edge-cell definition by constructing the concave hull of each CN and computed the G-Cross for these edge cells across CNs. More localized interactions were captured by the SCIMAP spatial interaction metric20, computed within each CN for each phenotype over a radius of 100 µm (Fig. 5a).

Testing of all engineered features against patient response using a traditional Mann–Whitney U test revealed that the univariate P values differed between feature types and feature families (Fig. 5b). After Benjamini–Hochberg P value correction, no features approached the false-discovery rate threshold, motivating more advanced statistical methods for the analysis to search for biomarkers and to accurately model patient outcomes. Particularly, interaction-related features such as JSD and G-Cross had higher variability in P values than cellular/cellular plus functional/metabolic ratios in compartments alone. As an example, an extreme value of JSD between macrophages and endothelial cells within the cohort are shown, where interspersion of the two cell types resulted in maximum overlap probability (left), whereas the self-clustering of macrophages away from endothelial cells resulted in minimal overlap in the core (right) (Fig. 5c). A profile of JSD and G-Cross features for cell-cell interactions in our analysis presented no obvious feature trends across patient response, highlighting the need for more functional or metabolic stratification in our feature space than base cell-type interactions (Fig. 5d).

Clinical benefit spatial feature selection
A unique challenge when analyzing spatial metrics are the numerous feature spaces generated. Here, over 1049 K engineered features were generated. To address this, the features were run through the Stabl package21, a biomarker and feature selection software designed to select statistically robust features that distinguish clinical outcomes from high-dimensional datasets. The package allows for stable selection of relevant features while also estimating a bound on the false-discovery rate (FDR) of uninformative features by artificial feature injection (Supplementary Fig. 3). Prior to feature selection, features were normalized with standard-scale normalization, missing values were imputed as zeros (with the exception of the JSD metric, where null values were imputed as 1), and finally grouped into feature family sets (cell proportions in K=2 CNs, JSD in K=2 CNs, JSD in K=3 CNs, G-Cross in MBNs etc.). Stabl was used to select features within each feature family (Fig. 6a).

From our high-dimensional feature space, 87 features were selected, of which 7 features were cell proportions within low activity MBNs, while the remaining consisted of JSD and G-Cross interaction values in tissue or metabolic regions. Regularization paths and false-discovery estimates from the feature selection method are shown in Supplementary Fig. 4, and feature distributions are shown in Supplementary Fig. 5. We focused on several of these cell-cell interaction feature sets to evaluate the interpretability of the model. From the G-Cross selected features that associated with worse clinical benefit within the tumor and interface regions, the proximity of tumor cells to fibroblasts, proximity of PD-1+ CD8+ T cells to PD-1+ macrophages, as well as the proximity CD4+/CD8+ T cells to other immune cells were selected.

Conversely, features associated with clinical benefit were immune cell proximity to fibroblasts, indicative of immune accessibility, granzyme B+ immune NOS cells in proximity to fibroblasts, and CD4+ Tregs in proximity to immune cells, facilitating priming, in the tumor and interface regions. Granulocytes near active tumor cells in metabolically low regions, and self-proximity of IDO1+ macrophages as well as functional tumor cell proximity to fibroblast cells, potentially showing tumor containment within stromal regions, were also positively associated. The maximum false-discovery rates of uninformative features were ~35%, 30%, 25%, and 20% for G-Cross features in K=2 CNs, K=3 CNs, across CN edge cells, and in MBNs, respectively (Fig. 6b, c).

Selected JSD interaction features associated with relapse within low metabolic or stromal regions were composed of PD-1+ T cells near tumor cells both expressing and not expressing HLA-A, T cells in proximity to immune NOS or endothelial cells, and tumor cells interacting with functional CD4+ T cells or Tregs. Within these regions, selected features associated with clinical benefit indicated higher densities of plasma cell infiltration. The estimated maximum false-discovery rate of uninformative JSD features was ~15%. Seven proportional features within low metabolic regions were selected, with an estimated false-discovery rate of 50%, where highly metabolic tumor cells, tumor cell sexpressing HLA-A, and granzyme B+ macrophages had a negative clinical effect. Here, IDO+ endothelial cells were associated with a positive clinical benefit (Fig. 6d–f).

Stabl selected 31 metabolic G-Gross features across the tumor, interface, and stroma, with a maximum estimated 50% false-discovery rate of uninformative features. Inactive macrophages near tumor cells within the tumor compartment were associated with negative clinical benefit. Inactive CD4+/CD8+ T cells near active tumor cells within the tumor and interface regions also had a negative association. Metabolically active macrophages near inactive structural or immune cells were also negatively associated (logHR>0). Glycolytic tumor cells near CD4+ Tregs were also associated with better clinical benefit (logHR<0) (Fig. 6g, h). These interaction features point to a dynamic signaling environment that underpins processes involved in both resistance and response to ICI.

Survival predictions based on selected features
Feature selection using binary endpoints enabled a robust statistical framework to be applied, however further modeling was required to assess the power of these spatial features to predict PFS events. Cox PH regression was performed on all selected features, providing log-Hazard Ratios for PFS (top) and OS (bottom). Only significant features are shown (Fig. 7a) where features whose range intervals did not cross the log HR=0 threshold were significant in the context of hazards. Top features indicative of progression included G-Cross interactions between tumor cells and macrophages, and granzyme B+ macrophage frequency in low metabolic activity MBNs. AUCs from k-fold Cox PH fits on PFS and OS indicated that the predictive AUC for OS is around ~0.7, and that the predictive AUC for PFS is around ~0.8 for a remarkable period of 24 months. A full list of selected features can be found in Supplementary Fig. 3a. A supplementary model based directly on time to event (PFS) feature selection identified many similar features (37/88), exhibiting common themes of metabolic pathways and interactions (Supplementary Fig. 3b–d) with similar predictive power (AUC ~ 0.9). The binary model is shown here for consistency with previous binary tests.

Furthermore, we sought to define the cross-cohort predictive power by fitting each cohort to prediction models using both the CB6 and time to event PFS selection method (CB6 model; cohort YTMA471 AUC 0.73, cohort YTMA404 AUC 0.65, PFS model; cohort YTMA471 AUC 0.87, cohort YTMA404 AUC 0.93) (Supplementary Fig. 3e). We additionally implemented a grouped bootstrap sampling method where only one cohort was utilized during each subsampling fit performed during the feature selection. This method pulled more significant features by equally weighting the YTMA404 cohort, which contained fewer patients but had lower average PFS times.

Notably, seven of these features were common to all three feature selection models (combined, YTMA404/YTMA471 independent) and were each significant by Kaplan–Meier tests, supporting generalizability of features across two independent cohorts (Supplementary Fig. 3e). Features that were predictive of benefit to ICI included the interaction between ICOS+ CD4 Tregs and fibroblasts in the tumor interface region, suggesting Treg exclusion, as well as CD4 Tregs interacting with glycolysis+ tumor cells in stromal regions, implicating a role for TME immune suppression of isolated tumor cells. Negatively associated features implicated macrophages in several scenarios, where self-aggregation of IDO1+ macrophages in low metabolic regions, and granzyme B+ macrophages in low metabolic regions were associated with poorer ICI outcomes (Supplementary Fig. 3e).

To summarize our feature selection models, features from the CB6 selection were combined to form a single prognosis score, where the HR of each feature was multiplied by the normalized feature value for each patient. Prognosis scores, split by median values, correlated with OS and PFS (Fig. 7b) although validation on a larger cohort is needed for cross-validation to confirm that this method gives an appropriate indication if a patient will progress or not (PFS p<1e−7, OS p=1e−7).

We further expanded this approach to derive a discrete signature that might describe either resistance (mean log HR>0, 43 features), or response (mean log HR<0, 44 features). These individual signatures indicated successful prediction of response for PFS (p=0.016), as well as for prediction of resistance for both OS (p=7e−6) and PFS (p<1e−7). Our study thus culminated in the successful development of both expanded and discrete feature sets that describe patient response to ICI therapy in our discovery cohorts of fifty-five total patients (Fig. 7c).

Methods
Patient selection
The tissue samples were collected and utilized under the approval of the Yale Human Investigation Committee (HIC), protocol #9505008219, with assurances filed with and approved by the U.S. Department of Health and Human Services. Retrospective NSCLC patient cohorts were generated at Yale School of Medicine35,36. Two independently collected cohorts obtained between 2011 and 2019, YTMA404 and YTMA471, were made from FFPE resected tumor samples from advanced stage NSCLC patients treated with ICI, and consisted of single 660µm and 600µm cores respectively per patient biopsy, with representative regions selected by pathologists in the Rimm Lab (Yale). Inclusion criteria were applied as follows: confirmed NSCLC by a registered pathologist at the time of resection, treatment with immune checkpoint inhibitors (with or without concurrent radiotherapy or chemotherapy), accessible prognostic follow up data (Supplementary Data 2). Clinical endpoints included information on overall survival (OS), progression free survival (PFS), as defined by time from start of ICI treatment, and best overall response (BOR) according to RECIST 1.1; complete response (CR), partial response (PR), stable disease (SD) and progressive disease (PD). Progression free survival above 6 months or objective response (partial or complete) was used to define clinical benefit from ICI (CB6).

Multiplexed immunohistochemistry using the PhenoCycler Fusion
TMA FFPE blocks were sectioned at 5µm onto SuperFrost plus slides. Sections were stained as described previously37 by Akoya Biosciences’ (USA) STEP program on the Phenocycler platform (Akoya Biosciences, USA). Target proteins, conjugated barcodes and attached fluorophores are listed in Supplementary Data 3. Briefly, sections were baked at 60°C for 1h to optimize tissue adherence, followed by deparaffinization in Histochoice Clearing Agent (VWR, H2779). Sections were hydrated in decreasing concentrations of ethanol and subjected to antigen retrieval buffer AR9 under pressure for 20 min (Akoya, AR900). Final images were obtained as qptiff files and pre-processed using QuPath for manual TMA dearraying. Dearrayed cores were exported in the ome.tif format and quality control (QC) was performed to remove defects (out of focus areas, large staining artifacts, and folded tissue). WSI areas that passed QC were then analyzed downstream and phenotyped by Nucleai’s deep-learning cell classification algorithms (Tel Aviv, Israel)38.

Cell segmentation
Multi-instance cell segmentation was performed using a deep learning model39. Nuclear segmentation was based on the DAPI channel and whole cell segmentation was based on the sum of all available membrane and cytoplasmic channels. Further post-processing was performed to match nuclear and whole cell masks, allowing the removal of segmented cells that did not contain nuclei, with nuclear merging of multinucleated cells, or the splitting of any nuclei assigned to multiple cells. In addition, a membrane segmentation mask was extracted by subtracting each nuclear mask from its matched whole cell mask, then regularizing it by a ring around the nucleus. Cell segmentation provided cell centroids for downstream deep learning tile classification, however, was not used in cell classification. Subsequent spatial calculations were performed on cell centroids.

Deep learning cell typing and cell state phenotyping
Multi-channel 54×54-pixel tiles were cropped around the center point of all annotated cell segmentation instances, were decomposed to single channel tiles, and were normalized using per-channel per-slide normalization as previously detailed38, to optimize the dynamic ranges for each channel on each slide, amplifying signal, reducing noise, and standardizing images across the cohort.

To infer marker positivity for each cell-marker combination, the normalized single channel tiles were fed to Nucleai’s single-channel classifier (DL binary classifier), which outputs the probability of marker positivity for each marker in every cell. Cohort-level single channel thresholds for binary classification were calculated by identifying the minima of the density plot of the positive prediction probabilities distribution (between 0.02 and 0.75 probabilities). Thresholds larger than 0.5 were reassigned as 0.5.

For subsequent cell typing from the single-channel binary classification of lineage proteins, we defined 14 mutually exclusive cell types based on known protein expression, in addition to an artifact cell class which was excluded from downstream analysis and was defined based on recurrent artificial proteins expression pattern which was visually inspected as an artifact (Supplementary Data 1). Segmented cells were assigned a cell type label if their predicted marker positivity profile exactly matched one of the pre-defined expected protein expression profiles. For cells that did not meet these criteria, cell type label was established by the most frequent cell type among their five nearest neighbors in the marker probability space of exactly matched cells.

Cellular and functional classification
Canonical markers were used to phenotype cells according to the following expression profiles; Pan-CK+/E-Cadherin+ (cancer cells), CD45+CD14+CD68+ (macrophages), CD45+CD14+CD68- (myeloid not otherwise specified (NOS) cells), CD45+CD31+CD34+ (endothelial cells), CD45+vimentin+ (fibroblasts), CD45+SMA+ (myofibroblasts), CD45+CD3- (immune NOS cells), CD45+CD3-CDllb+ (granulocytes), CD45+CD3-CD20+ (B cells), CD45+CD3-CD20-CD31+ (plasma cells), CD45+CD3+CD4+FOXP3+ (CD4+ Tregs), CD45+CD3+CD4+FOXP3- (CD4+ T cells), CD45+CD3+CD8+ (CD8+ T cells). Non-tumor cell functionality was classified into PD-1, PD-L1, Granzyme B, ICOS and IDO1 expression. Tumor cell functionality was classified according to PD-L1, Vimentin, Ki67, IDO1 and HLA-A expression.

Tissue area and neighborhood classification
Cellular neighborhoods (CNs) were identified as previously described40. Briefly, a window of each cell’s K nearest neighbors (KNN) in physical space was used to construct a neighborhood histogram representing the frequency of observing a neighboring cell of a given phenotype. K-means clustering was then performed on this histogram to assign each cell a CN. The method was adapted to identify minimalistic CNs, where coarser tumor and non-tumor cell phenotype classifications were used in place of cell types. With this, two different minimalistic CN views were produced. The first was a tumor-stroma CN view that was computed by clustering the neighborhood histogram from a KNN window of N=50 into K=2 CNs. The second was a tumor-interface-stroma CN view that was computed by clustering the neighborhood histogram from a KNN window of N=30 into K=3 CNs.

Metabolic neighborhood definition
Derived from CNs16, metabolic neighborhoods (MBNs) were defined by grouping cells based on similar cellular metabolic densities using KNN windows and K-Means clustering. Each cell was annotated as positive or negative over 7 metabolic pathways based on the positivity of a proxy metabolic marker (Amino acid uptake: ASCT2, ATP Synthesis: ATPA5, TCA Cycle: Citrate synthase, Fatty acid oxidation: CPT1A, pentose phosphate pathway: G6PD, Glycolysis: GLUT1, OxRegulatory: pNRF2). Metabolic neighborhood histograms were constructed by summing the number of neighboring cells from a KNN window of K=30 that were positive for each metabolic pathway. K-means clustering was then performed on this histogram to cluster cells into K=4 MBNs. The choice of four neighborhoods was motivated to select a minimal number of MBNs that distinguish regions of minimal, low, medium, and high metabolic activity. Additionally, the cluster inertia was assessed at each K, and elbow fitted to determine optimal number of clusters (Supplementary Fig. 6a) where we additionally found no significant gain in the number of differentiable regions when using a larger number of clusters.

G-cross function
To quantitively determine cell-cell and cell-function interactions, the cumulative distribution of pairs of cell and cell plus functionality types were computed using the G-cross function. This measures the cumulative distribution of the nearest neighbors of each target cell or cell plus functional/metabolic type with respect to each reference type. The Area Under the Curve (AUC) was computed for this distribution up to 150 µm from each reference type. Internal testing indicated that 150 µm captured a majority of cell type proximity characteristics while representing a biologically meaningful scale (Supplementary Fig. 6b).

Jensen-Shannon Distance (JSD)
Cell-cell interactions and their correlation within individual cores were determined using a JSD distance metric, refined from the DIMPLE17 R package, to compute pairwise distance matrices. This method applies the 2-dimentional spatial intensities of two cell types from the same image, based on cell typing, where cell-cell Gaussian kernel density estimates are used to profile spatial density and measure the probability density interactions of cell types. JSD is bound by complete overlap (0) and complete separation (1) and is calculated between each pair of normalized intensity functions.

Cellular proportions
The proportions of cellular phenotypes normalized to each CN compartment were computed. This was repeated independently for every combination of three cellular phenotype views (cell types, cell types with singular functional marker positivity, and cell types with combined metabolic marker positivity) within the three CN views (the tumor-stroma, tumor-interface-stroma, and metabolic neighborhoods). The ratio of cell counts between pairs of base cell types were computed within the tumor/stroma/interface compartments per patient and compared between CB6 groups by Mann–Whitney U tests.

Multivariate analysis with Stabl
Multivariate predictive modeling improves feature selection by leveraging the combination of many pair-wise interacting features to model complex biological interactions. Stabl21 was implemented as a framework specific to clinical translation of complex omics data into future predictive clinical outcomes. The software performs feature selection over sparse and high-dimensional data by artificial feature injection during stability selection to select features at a minimal false discovery rate. The process avoids applying conservative false-discovery correction methods that are discovery-prohibitive in the high-dimensional regime of this analysis. Features were imputed and normalized, and the feature columns were filtered to those with at least twenty-five unique values, to reduce feature sparsity. Features were selected by Lasso fits to the CB6 variable or by Cox PH fits to the PFS. Feature selection was performed with either knockoff or random permutation artificial feature injection in the Stabl framework. The selected features from each feature family were combined and used to fit regularized Cox Proportional Hazards model to form predictive models of patient response and overall survival. For the grouped cohort fits, a bootstrapping function was used to randomly select a 50% fraction of samples from each cohort for each subsample fit. To estimate CB6 from the PFS-selected features, logistic regression fits to the selected features were performed over k-folds.

Clinical endpoint analysis
Univariate features. Potential differences between the two CB6 patient groups with respect to each proportion and cell count ratio feature were evaluated independently using multiple Mann–Whitney U-tests. False discovery rate (FDR) p values accounting for multiple testing was performed using the Benjamini–Hochberg correction. Correlation with progression-free survival was also performed independently for each feature using a Cox partial hazards model.

Survival analysis. Selected features were agglomerated for all feature families and used to fit a Cox Proportional Hazards model to the censored progression free survival (PFS) and overall survival (OS) data. 10-fold l1-regularized Cox PH fits (75%–25% split) to the features provided time-dependent AUCs of both PFS and OS from the test sets, providing an indication of the models’ prediction accuracies on held out data. Confidence intervals on the 95%–5% prediction ranges over these k-fold fits were additionally derived and plotted.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "1. The study explores metabolic pathways and their interactions within the TME, a unique and novel approach using multiplexed techniques.\n2. Non-tumor and tumor cell functional markers depicted in Figure 1 might not be entirely appropriate; for instance, PD-L1 can be a tumor cell function marker, and Ki67 can be a marker for both tumor and non-tumor cells.\n3. Including a cell lineage/state marker expression heatmap similar to Figure 2b would enhance the clarity of the analysis.\n4. Regarding lines 140-144, clarification is needed regarding whether the discussion about the cellular makeup of CB6 response groups refers only to the clinically benefit group.\n5. Lines 147-148 state \"with much higher levels of positivity for metabolic functional states,\" and line 151 mentions \"macrophages showed significant positivity for all functional markers.\" Are there any p-value significances for such cell type enrichments?\n6. In Figure 2c, the color legend for the p-value (*) is missing.\n7. In Figure 2g, was the survival analysis based on the same cohort? If so, this might simply reflect the association that the clinically non-benefit group had shorter progression-free survival (as expected), rather than demonstrating the predictive value of granzyme B+ macrophage abundance.\n8. In the selected PR and PD tissue samples shown in Figure 3b, the interface seems to be overestimated. It might be worth investigating whether this is commonly observed in tissues with small isolated tumor islands. Additionally, providing serial/adjacent H&E images for confirming the accuracy of tissue region categorization would be helpful.\n9. Lines 201-215 use \"patient relapse\" and \"worse clinical outcome\" interchangeably, which is confusing. Do they refer to the same group of patients?\n10. A systematic and comprehensive TME feature engineering approach represents a useful contribution to the community. However, automating the determination of the number of clusters (K=2, K=3, MBNs=4) could enhance the manuscript.\n11. Regarding lines 301-303, the G-cross function was computed for radii up to 300um. Is there a rationale for choosing this distance threshold? What would be a physiologically meaningful cell-cell distance?\n12. What do the maximum FDR values (lines 368-369, 378, 379) signify in the context of Stabl-derived features?\n13. In Figure 6h, within the tumor region (top left pie chart), the label for the '1' component is missing.\n14. Regarding lines 406-432, were the features selected by Stabl in a supervised manner for survival analysis? If so, this analysis might be biased.\n15. To fulfill the objectives of TME characterization and biomarker identification for ICI treatment in lung cancer, more biological interpretation of the identified TME features is required.", "reviewer_2": "In this manuscript, Monkman et al. present a significant and original study characterizing the tumor-immune microenvironment to predict immunotherapy response in NSCLC. The methodology is sound, combining high-plex mIF, deep learning, and a spatial feature engineering pipeline with robust statistical selection. The development of a multivariate model predicting progression-free survival over 24 months with high accuracy (AUC=0.8), based on spatial interaction and metabolic features, is a particularly noteworthy achievement. The work is of high interest; the following points are offered to further strengthen the manuscript for publication.\nMajor point\n1. This is a valuable study, and its impact could be maximized by making the acquired images and analysis code publicly available. This practice greatly benefits the field by fostering transparency and reproducibility. The URL provided in the reporting summary was found to be inaccessible at the time of review. Verification and updating of this link are requested during revision.\nMinor point\n2. The discussion on model robustness against tumor heterogeneity could be strengthened. A comment on the model's performance consistency across different samples from the same patient (if such data exist) would be a valuable addition to substantiate the findings.\n3. The Discussion would be enriched by a more detailed comparison with established predictive methods (e.g., PD-L1 score, TMB, or other spatial analysis or omics-based approaches). This would allow readers to better appreciate the advantages of the methodology presented in this work.\n4. The readability of several figure panels requires improvement, as the text is currently difficult to decipher. Specifically, the text within Figures 3e and 3f, Figure 5d, Figures 6b, 6d, and 6g, and Figure 7a is challenging to read. Enhancing the clarity of these figures would greatly benefit the reader's understanding of the results.", "reviewer_3": "This is an interesting work with specific target to reveal the spatial mechanisms of immunotherapy response in lung cancer based on MPX data. My detailed comments below:\n1) the overall writing is clear, but some parts are very wordy and convoluted. Please furthere proofread.\n2) The survival prediction approach first selects features based on binary labels, then filters them using continuous survival data via Cox regression. It would be interesting to compare this with a method where feature selection is performed directly using time-to-event data, rather than relying on manual dichotomization.\n3) The exploration of spatial mechanisms in this work heavily depends on cell segmentation and classification. However, technical details regarding the robustness and accuracy of these methods are limited. Given the substantial variation and dynamic range typical of MPX images, more comprehensive evaluation of image processing quality is needed.\n4) The derivation of cell metabolism-related markers associated with survival and treatment response is compelling. Nevertheless, these findings have not been validated on an independent cohort. Given the high-dimensional nature of MPX data, external validation is essential to confirm the reliability and generalizability of the results from the discovery cohort."}, "author_response": "Response to Reviewers comments:\nThe authors thank the 3 reviewers for their time in reviewing our manuscript and for their constructive\ncomments and suggestions. We have addressed all points below (our responses are highlighted in red) in\nresponse to each comment below.\nTo improve the discovery of the paper for our target audience, we have updated the title to include:\nLines 3-4:\n“in non-small cell lung carcinoma (NSCLC)”\nReviewer #1:\n1. The study explores metabolic pathways and their interactions within the TME, a unique and novel\napproach using multiplexed techniques.\nWe thank the reviewer for this summary of the manuscript and for highlighting the novelty of our\napproach.\n2. Non-tumor and tumor cell functional markers depicted in Figure 1 might not be entirely appropriate;\nfor instance, PD-L1 can be a tumor cell function marker, and Ki67 can be a marker for both tumor and\nnon-tumor cells.\nWe have now altered figure 1 to avoid confusion with functional marker classifications of tumor or non\ntumor cells both in figure 1c and with addition of 1g. Ki67 label has been amended from ‘Tumor’ to\n‘Division.’ This has also been amended in the main text at:\nLine 113:\n“dividing cells”. PD-L1 was included as both tumor and non-tumor functional marker in original figure 1h.\n3. Including a cell lineage/state marker expression heatmap similar to Figure 2b would enhance the\nclarity of the analysis.\nWe have generated a cell lineage/state marker expression heatmap for better clarity and replaced panels\nFig1 g-h with this heatmap. These panels are now Supplementary Figure 1 (a-b), indicated in line 122. We\nhave added reference to this in:\nLine 122-124: “Overall expression of functional markers are shown per cell type, where metabolic\nphenotypes predominated in both immune and tumor cells (Fig.1g)”.\nThis has affected subsequent supplementary numbering, hence supplementary figure 1 is now\nsupplementary figure 2 (line 241); supplementary figure 2 is now supplementary figure 3 (line 429 and\n432); Supplementary Figure 3 is now supplementary figure 4 *line 364); Supplementary figure 4 is now\nsupplementary figure 5 (line 364). The figure legend for Figure 1 has been updated to reflect the changes\n(lines 137-139).\n(Left is Figure 1 before revision, right is Figure 1 after revision)\n4. Regarding lines 140-144, clarification is needed regarding whether the discussion about the cellular\nmakeup of CB6 response groups refers only to the clinically benefit group.\nThese lines have been amended to avoid confusion regarding clinical comparisons, as this paragraph\nonly explores cell-state expression patterns, not clinical endpoints:\nLines 142-144:\n“The gross cellular makeup of patient tissues was first inspected for cellular frequencies and functional\nexpression patterns. Tumor cell composition ranged from 10-90% within the cohort and isolation of non-\ntumor cells revealed macrophages…”\n5. Lines 147-148 state \"with much higher levels of positivity for metabolic functional states,\" and line 151\nmentions \"macrophages showed significant positivity for all functional markers.\" Are there any p-value\nsignificances for such cell type enrichments?\nThe purpose of figure 2b is to introduce the patterns of cell type functional classification (on the log\nscale) within the entire cohort. To avoid confusion, this has been amended:\nLines 146-148:\n“Overall abundance of the functional status of cell types was measured (Fig. 2b), indicating tumor cell\nexpression patterns for functional markers (IDO1, PDL1, HLA-A, Ki67 and vimentin)”.\n6. In Figure 2c, the color legend for the p-value (*) is missing.\nWe have amended the figure legend to reference the color legend for Fig.2c on:\nLine 186-187: “Significant values (*) are shown, those in yellow are further evaluated in d (p<0.05)”.\n7. In Figure 2g, was the survival analysis based on the same cohort? If so, this might simply reflect the\nassociation that the clinically non-benefit group had shorter progression-free survival (as expected),\nrather than demonstrating the predictive value of granzyme B+ macrophage abundance.\nOur approach aimed to identify overall functional cell frequencies that associated with patient\noutcomes. The progression free time characteristics of our clinical benefit at 6 months is shown in Fig 1e.\nFigure 2 demonstrates the association between cell frequency and both response at 6 months (CB6) as\nwell as time to event (PFS) analysis. Increased granzyme B+ macrophages appeared to have an\nassociation with both CB6 and PFS time, and Fig.2g illustrates the association by splitting the cohort at\nmedian frequency levels. We agree with the reviewer that this may not be predictive but it is rather a\nreflection of this association within this cohort. Further cross validation predictive modelling of granzyme\nB+ macrophages is provided later in the paper, where this feature within tumor regions is selected and\npassed STABL feature selection criteria in both cohorts (Fig.3d-e).\n8. In the selected PR and PD tissue samples shown in Figure 3b, the interface seems to be\noverestimated. It might be worth investigating whether this is commonly observed in tissues with small\nisolated tumor islands. Additionally, providing serial/adjacent H&E images for confirming the accuracy\nof tissue region categorization would be helpful.\nIn spatial analyses using cellular neighborhoods (CNs) (Schurch, Cell 2020), the neighborhood definition\nis formed from each cell’s KNN neighborhood phenotype histogram. In this case, the “interface” is not\nnecessarily the bordering cells that separate tumor cells from everything else. The tumor-stroma\n“interface” will be composed of a harmonic balance of tumor and non-tumor cells encompassing the\nregions surrounding the edge of the tumor and encroaching into the stromal region. The border region\n(i.e. direct interface cells) is more likely to be captured by the k=2 CN edge cells defined in this analysis\nvia the concave hull of the CN boundaries.\nIn the case of isolated tumor nests, the nest would not be classified within the Tumor CN if the relative\nproportion of neighboring cells did not match those of a larger primary tumor. For example, if the CN\ndefinition (assuming k=2) was for the distribution of the 50 closest cells, and the composition in the nest\nwas 20 tumor cells and 30 non-tumor cells, this would fall below the harmonic balance of tumor/non-\ntumor cells and would likely fall under the Stroma CN. This is why the k=3 tumor-stroma interface CN is\ncritical, as it allows the method to be more sensitive to dynamics along the interface where more\nflexibility in the cell-type distribution is allowed, instead of being dominated by interactions close to the\nbulk tumor mass. We have clarified this aspect of the analysis in:\nLines 200-203 stating “The K=2 CNs were designed to capture tumor and stromal/immune regions,\nwhereas the K=3 CNs were designed to provide resolution to the “interface” region between structural\ncells and the tumor, where functional immune cells and metabolic gradients are expected to play a\ncritical role in the clinical response”.\nThis highlights that this interface region captures tumor cells on the edge of the tumor bulk along with\nimmune cells and structural cells – we posit that these play an outsized contribution to the dynamics of\nthe tumor microenvironment and compute spatial metrics in this constrained region as a result. Capture\nof this critical region is also important in the context of tumor nests, where immune action against\nstromal-invasive tumor regions may play a role in preventing tumor outgrowth.\nWe evaluated H&Es available for these cohorts, but the serial sections were too far apart to make a\nmeaningful comparison between annotated tumor regions and the CNs. However, when examining the\nraw immunofluorescence data of PanCK positive regions, we see that bulk tumor regions are captured in\nthe Tumor CNs, whereas invasive margins, bulk edges, and neighboring immune cells are captured\nwithin the Tumor/Stroma Interface CNs, highlighting how metrics within this region will be enriched in\ninformation containing immunological-tumor interactions that contribute to the downstream clinical\nmodeling. We have included the below images from cores in Fig. 3b and updated the figure legend\naccordingly:\nLines 247-247: “Immunofluorescent images of core positivity for PanCK (cyan) and DAPI (blue) are shown\nto illustrate neighborhood concordance”.\nAbove, PanCK (cyan) and DAPI markers (blue) on the left subpanels. Tumor CN (orange), Tumor/Stroma\nInterface CN (green), and Stroma CN (blue) annotated cell locations on the right subpanels, showing that\nstromal regions are captured by the Stroma CN and tumor-interacting regions are captured by the\ninterface.\nWe have modified Figure 3 and the figure legend to include fluorescence images of the PanCK-stained\ntumor regions. Before (left) and after (right) changes to this figure are shown above.\n9. Lines 201-215 use \"patient relapse\" and \"worse clinical outcome\" interchangeably, which is confusing.\nDo they refer to the same group of patients?\nWe thank the reviewer for picking up this discrepancy. We are indeed referring to the same group of\npatients. To clarify this, we have replaced the word “relapse” with “worse clinical outcome” (line 210-\n211).\n10. A systematic and comprehensive TME feature engineering approach represents a useful contribution\nto the community. However, automating the determination of the number of clusters (K=2, K=3, MBNs=4)\ncould enhance the manuscript.\nThe choice for k=2 and k=3 CNs is first principles based, and we have clarified in Lines 203-206 on the\nreasons for this, as in comment #8 above. While the choice in number of CNs for k > 3 can be automated\nby examination of the Jaccard Index or K-means inertia as a function of the “k” and picking the optimal\nnumber of clusters at the elbow-point, CNs for k > 3 are not-meaningful for this analysis as the\nproportions and distance metrics become extremely sparse for tissue cores of this size. This greatly\nhampers both the feature selection method and interpretability of the selected features. Cellular\nneighborhoods of these types are also typically derived using cell phenotypes rather than tumor/non-\ntumor phenotypes and are typically used to find unique cell-type distributions which can be quite\nspecific to the dataset and subject to interpretation. In contrast, our approach of finds unique cell-cell\ninteractions within minimalistic neighborhoods that are simpler yet familiar and interpretable to the field\n(tumor, stroma, interface, etc.) and enriched in specific immunological and metabolic processes in a\ncore-standardized manner.\nThe MBNs defined in this analysis are not defined by exclusive labels, unlike cell-types in the CNs, so the\noptimal point was chosen by examination of heatmaps of the total metabolic pathway expression.\nCluster inertia at each K was also assessed to aid the choice of the optimal number of clusters (elbow at\nk=4). We found that choosing a number of MBNs greater than 4 yielded no further gain in differentiable\nclusters. We also sought for a small number of differentiable clusters that did not cause a similar sparsity\nissue as the one seen when one chooses a large number of CNs.\nWe have clarified this in text:\nMethods Line 702:\n“The choice of 4 neighborhoods was motivated to select a minimal number of MBNs that distinguish\nregions of minimal, low, medium, and high metabolic activity. Additionally, the cluster inertia was\nassessed at each K, and elbow fitted to determine optimal number of clusters (Supplementary Fig.6)\nwhere we additionally found no significant gain in the number of differentiable regions when using a\nlarger number of clusters”.\nSupplementary Figure 6a: MBN cluster inertia with elbow fit to determine K-means fit\n11. Regarding lines 301-303, the G-cross function was computed for radii up to 300um. Is there a\nrationale for choosing this distance threshold? What would be a physiologically meaningful cell-cell\ndistance?\nAs a feature, G-Cross AUC is designed to capture interactions at a variety of scales, making the choice of\nradius somewhat arbitrary – assuming that this radius will cover most meaningful cell-cell pairwise\ninteractions. Internal tests show that if the G-Cross function plateaus at a certain radius, for example 200\nµm, then the G-Cross AUC is typically the computed value up to 200um plus that additional “100 um”.\nThis makes the radius of this feature somewhat invariant in the subsequent feature selection for plentiful\ncell types due to the standard scaling that’s performed in the pipeline. This “100 um” extra will be\nremoved in the mean subtraction calculation of the feature normalization. The rationale then for not\nchoosing an arbitrarily large radius is that the calculation may be biased by the missing data at the edge\nof the tissue/core, so a radius is chosen to be a fraction of the core radius.\nAn optimized physiologically meaningful distance would need to be rigorously defined and is outside the\nscope of this paper – an optimal distance would be cell-pair-type specific, dependent on the core radius,\nand dependent on both expected biological interactions and on the sparsity of the cell-type pairs for each\ncomparison. For example, cell pair interactions such as receptor-ligand interactions on cell surface have\ndifferent interaction ranges than the suppressant effects of T-Regs in the stroma, and metabolic\ngradients exist in the tissue at a range of scales, making a cell-pair specific radius difficult to generalize\nto other tissue/disease types. The biases existing in calculations of ecological metrics of this type for rare\ncell types are known but is an unsolved problem in the field. Biases in the G-Cross calculation exist in\nedge cases, especially at larger radii, limiting the range of interactions that can be explored and\nadditionally enforcing a cohort/study-specific approach.\nBelow, we show a grid plot of G-Cross curves between major cell types as target and reference types,\nshowing that this metric captures the vast majority of potential interactions at that choice of radius\ncalculation, and demonstrating that this metric containing characteristic shapes\nWe have clarified this choice in the manuscript:\nLines 310:\n“Since this feature captures nearest-neighbor distances, this radius was chosen to be above the scale of\nexpected local interactions, but less than the core radius to reduce biases by missing cells outside of the\ncore boundary”.\nIn addition, the distance unit reported in Line 310 is in pixels and this has been corrected to “150 µm,” to\nbe consistent with the correct value/unit that is stated in the methods section.\nAdditionally in methods:\nLine 715:\n“Internal testing indicated that 150um captured a majority of cell type proximity characteristics while\nrepresenting a biologically meaningful scale (Supplementary Fig 6b)”.\n12. What do the maximum FDR values (lines 368-369, 378, 379) signify in the context of Stabl-derived\nfeatures?\nIn the context of Stabl-derived features, maximum FDR refers to the upper bound on the false discovery\nproportions within that selected group of features, where the false discovery proportion is defined as the\nproportion of “uninformative” features that were selected. This provides a more rigorous method of\nfinding clinically relevant features as compared to univariate tests, as linear model fits to clinical\nendpoints can capture not just mean/median changes in single feature values or in clinical correlation\nbut can jointly capture multivariable correlations. By using this method, the number of features that\nhappen to align towards clinical relevance as a result of random statistical fluctuations in this high-\ndimensional feature space will be much less than that fraction of total selected features. The authors of\nthe Stabl software (Hedou et al., Discovery of sparse, reliable omic biomarkers with Stabl. Nature\nBiotechnology 2024) show that their method consistently estimates the false-discovery proportions\n(fraction of uninformative features) of selected features. We show in the paper that under traditional\nfalse-discovery correction methods that tend to be quite conservative, no features could be classified as\nclinically relevant when considering Mann Whitney Tests of the log-fold change of each of the features.\nHowever, our analysis shows that the estimated maximum false discoveries of selected features within\neach of the feature families are a small fraction of each family, demonstrating the power of this\nmodelling-based method for controlling the false-discovery rate of selected features in high dimensional\nspatial analyses.\nWe have clarified this in the paper by referencing “informative/uninformative features” in the “Clinical\nbenefit spatial feature selection” section of the paper, Lines 346-402, and in the Fig.6 figure legend line\n406.\n13. In Figure 6h, within the tumor region (top left pie chart), the label for the '1' component is missing.\nThis has been rectified in figure 6h, the label “Other” has been added to the top left pie chart of panel h to\nindicate the selected feature for the component in green with “1” within the tumor region.\n14. Regarding lines 406-432, were the features selected by Stabl in a supervised manner for survival\nanalysis? If so, this analysis might be biased.\nThe features selected by Stabl were derived from linear fits to the clinical benefit after 6 months (CB6)\nendpoint and independently validated on the survival information in k-fold fits using a train-test split. All\nfeatures were input into the model, and we have intentionally prioritized features that predict for CB6\nusing this approach and so some overlap with PFS is expected and desirable. The features are correlated\nwith the survival endpoint data (by implicit study design), however, the CoxPH and BoostedCox model fits\nto the set of normalized and selected features are performed to fit to the as-of-then unseen censored\nprogression-free-survival and overall-survival data only after the linear fits are performed for feature\nselection. The inevitable correlation between CB6 and the censored fit data may also manifest itself as\nthe feature selection and k-fold validation sets overlap in terms of the “seen data,” however, the applied\nlogic here is no different from many similar studies in which T-tests are done to find univariately-\nsignificant features, followed by Kaplan-Meier analysis of the features of interest. In contrast to those\nmethods, though, the combined modelling of the selected features allows for joint modelling of both\nunivariately-significant features and insignificant features, without constraints of conservative false-\ndiscovery correction methods. This allows for complementary looks at the selected features, both as a\ncohesive profile of the TME dynamics via an examination of all the selected features jointly, and\nindividually to search for univariately-significant features that may point to biological effects of defined or\nisolated mechanisms that may be individually prognosticative features. Additionally, the linear fits to\nCB6 are performed to random samplings of the two cohorts, meaning that it doesn’t merely select the\nmost univariately-significant features across the global dataset, but learns features that are more\ninformative for each sample, making it less biased towards picking random features that happen to\nfluctuate and align with both CB6 and PFS in the total dataset.\n15. To fulfill the objectives of TME characterization and biomarker identification for ICI treatment in lung\ncancer, more biological interpretation of the identified TME features is required.\nWe have identified a technical pipeline for spatially relevant biomarker identification within a NSCLC\ncohort as a discovery tool that might be applied to a variety of datasets. Our interpretation of features is\nbased on the available literature, where comprehensive mechanistic studies from identified features\nrequires experiments beyond the scope of this manuscript.\nWe have made several additions to extend the interpretation of our feature modelling:\nLines 521-558\n“… In addition to multivariate modelling of all samples to gain maximum statistical power, feature\nselection highlighted several properties that predicted outcome in both cohorts independently.\nTumor cells with high amino acid uptake proximal to macrophage cells within the tumor region\nwere predictive of poor outcomes for patients. Tumor cells in this state deplete the TME of\navailable amino acids, impairing immune cell functionality, and the ability for ICIs to reprime their\nactivation. Such a tumor cell metabolic state would promote M2-like macrophage\nimmunosuppression, through mechanisms such as kynurenine pathway activation from high\ntryptophan metabolism [32]. Metabolic stress also alters macrophage functionality impairing\nantigen presentation and cytokine signaling. T cells entering these tissue regions, despite having\nblocked PD-1/PD-L1 or CTLA-4, will succumb to dysfunction and create cold tumor regions\nunresponsive to immunotherapy alone. Further evidence of the impact of low nutrient availability\nis found in our data with concordant self-aggregation of IDO1+ macrophages in such low\nmetabolic tissue regions being predictive of poor outcomes. IDO1 expression by macrophages\nhas been reported to influence both M1 pro-inflammatory [33] and M2 anti-inflammatory [32, 34]\npolarized states. Moreover, IDO1 is responsible for tryptophan degradation and the subsequent\nproduction of its metabolite kynurenine, inhibiting T cell expansion while promoting Treg\ninduction, driving an overall immune suppressive TME [35] Together, these features may act in a\nfeedback system to minimise available resources for immune activation.\nAdditionally, granzyme B+ macrophages found in low metabolic activity neighborhoods predicted\npoor outcomes. The cytotoxic activity of granzyme B requires close cellular contact to enter\ntarget cells through perforin, primarily inducing caspase-mediated apoptosis. Univariate\nanalysis of broad tissue regions suggested their role in within tumor regions, while further\ndelineation of MBNs pointed to their Isolation within metabolically inactive neighborhoods,\nindicating a subset of tumor regions which might implicate areas of nutrient depletion or cellular\nquiescence where their cytotoxic effect is curbed.\nAmong positively associated features common to both cohorts, ICOS+ Tregs in proximity to\nfibroblasts in the tumor interface region associated with better ICI outcomes. ICOS expression\non Tregs promotes their immune suppressive activity [36], and their localization with fibroblasts\nin the tumor interface implicates the exclusion of these cells from the tumor gland in ICI benefit.\nThe proximity of Tregs with glycolysis+ tumor cells also associated with better ICI outcomes,\nperhaps implying a relationship between these cells that favors ICI sensitivity, such as feedback\nfrom Warburg effect [14] tumor cell metabolites and Treg directed immune suppression.\n”.\nReviewer #2\nIn this manuscript, Monkman et al. present a significant and original study characterizing the tumor-\nimmune microenvironment to predict immunotherapy response in NSCLC. The methodology is sound,\ncombining high-plex mIF, deep learning, and a spatial feature engineering pipeline with robust statistical\nselection. The development of a multivariate model predicting progression-free survival over 24 months\nwith high accuracy (AUC=0.8), based on spatial interaction and metabolic features, is a particularly\nnoteworthy achievement. The work is of high interest; the following points are offered to further\nstrengthen the manuscript for publication.\nMajor point\n1. This is a valuable study, and its impact could be maximized by making the acquired images and\nanalysis code publicly available. This practice greatly benefits the field by fostering transparency and\nreproducibility. The URL provided in the reporting summary was found to be inaccessible at the time of\nreview. Verification and updating of this link are requested during revision.\nWe thank the reviewer for their recognition of our work.\nWe have updated the link and the code is now available publicly at\nhttps://github.com/clinicalomx/metabolic-microenvironment-predictors-of-nsclc-immunotherapy-\nresponse.\nThe data is now openly available without login request available at https://doi.org/10.48610/73b218c.\nWe have included data and code availability statements in the main manuscript in:\nLines 774:\n“Data and code availability\nData used in this study is available at https://doi.org/10.48610/73b218c . The code used to undertake the\nanalysis is publicly available and has been deposited in the GitHub repository\nhttps://github.com/clinicalomx/metabolic-microenvironment-predictors-of-nsclc-immunotherapy-\nresponse”.\nMinor point\n2. The discussion on model robustness against tumor heterogeneity could be strengthened. A comment\non the model's performance consistency across different samples from the same patient (if such data\nexist) would be a valuable addition to substantiate the findings.\nWe appreciate the reviewer’s thoughts around the model robustness across multiple sampling within\none patient tissue block, however this sampling strategy was not captured in the study design by the\npathology team, especially for lung cancer where initial diagnostic biopsies are often quite small and\nlimited tissue is available. The study included single core biopsies per patient over two cohort collections\nfrom 2010-2019 from NSCLC presenting to Yale University School of Medicine.\nIntra-patient heterogeneity remains a limitation in these studies at this stage, however demonstrating\nconserved inter-patient reproducibility is more clinically relevant and thus the NSCLC TMAs were built to\ninclude one core per patient. In a clinical setting, single biopsy samples are routinely used in the first\ninstance to guide treatment. Our approach to biomarker discovery rests on feature selection models that\nprovide statistical power across diverse samples. To address this limitation, we have included a\nstatement in the discussion:\nLines 581-585:\n“Our study included one tissue sample per NSCLC patient; intra-patient heterogeneity poses a limitation\nin the study design. We instead sought to demonstrate inter-patient feature selection to reflect a\ndiscovery pipeline that may be more applicable in a clinical setting where single biopsy samples are used\nin the first instance to guide treatment”\nWe believe that multiple core sampling strategies may be useful in future studies and have started to\ndesign this prospectively. However, this comes at a technical cost across tissue real estate (imageable\narea on glass slides for spatial experiments), destruction of the whole block (by multiple coring), the\nnumber of patient samples one can get onto a glass slide using a multi-core strategy and the number of\nslides needed to capture the cohort (potentially this 2 TMA slide project could become a 5-10 slide TMA\nproject) which can increase technical variation, batch effects and compound the analysis (including\nexperimental costs).\n3. The Discussion would be enriched by a more detailed comparison with established predictive\nmethods (e.g., PD-L1 score, TMB, or other spatial analysis or omics-based approaches). This would allow\nreaders to better appreciate the advantages of the methodology presented in this work.\nPD-L1 IHC score, TMB and similar clinical measurements (e.g. MSI) are known biomarkers with limited\nutility clinically, often benefiting a very small subset of patients receiving immunotherapies. There\nremains no consensus as to the underlying mechanisms of resistance which remain to be clarified –\nhence the need for more predictive biomarkers for NSCLC and immunotherapy, likely informed by the\ncomplex tumor microenvironment.\nWhilst these are desirable biomarkers to have as part of a multi-omic clinical study, it was not possible to\nobtain sufficient tissue from TMA cores to isolate DNA to enable a TMB measurement (typically 5-10\nserial whole slide section are needed for mutational analysis, and our core sizes are a lot smaller than\nindividual whole cut sections). Rather, in our study, we sought to investigate spatially resolved proteins\nand cellular phenotypes from the TME. Spatial resolution of the TME, its interaction with innate and\nadaptive immunity, and its context specific state offer a more comprehensive characterization of the\ndynamic TME for cellular pattern recognition across cohorts, therapeutic groups and biomarker\nidentification. While other omics-based approaches hold value, we sought to identify the spatial context\nof cellular metabolic states through functional protein targets that might feasibly be translated and\nimplemented into clinical settings where this infrastructure is already in place.\nIn future studies and within translational endpoints for clinical trials, there is value in incorporating\nspatial TME assessments as part of the multiomics based readouts. However, providing the evidence\nthat spatial TME measurements are informative and a workflow to derive spatial TME biomarkers is\nimportant, as presented in this study. We envisage in the future that PDL1 score, TMB, interferon gene\nsignatures and spatial TME’s in addition to clinical data alone will be compared and contrasted for the\nability to predict responses.\nTo clarify this, we have amended the discussion:\nLines 488-500:\n“Current companion diagnostic assays for ICI vary in their accuracy and utility [24], where PD-L1 score\nalone is inconsistent with predictive response to ICI despite its current clinical use. In the KEYNOTE-001\ntrial, 10% of patients who didn’t express PD-L1 also responded. Tumor mutational burden (TMB) is also\nused clinically as a biomarker of response to immunotherapy, but like PD-L1, TMB based treatment\nstratification is not definitive. Non-synonymous mutations across different types of cancer and their lack\nof standardization make its sole use as a biomarker difficult to implement. Thus, current stratification\nassays suffer from relatively poor performance [25, 26]. Meta-analysis across ten tumor types indicated\nthat PD-L1 IHC (AUC 0.65) and TMB (AUC 0.69) are less accurate than mIF methods (AUC 0.79) in their\nability to predict ICI outcomes [27]. Hence, the development of a robust set of biomarkers for the likelihood\nof ICI response may shift the treatment paradigm towards a more personalized medicine approach for\nNSCLC patients. In our study, we sought to characterize the NSCLC TME for spatially informative\nbiomarkers through comprehensive cellular profiling.\n4. The readability of several figure panels requires improvement, as the text is currently difficult to\ndecipher. Specifically, the text within Figures 3e and 3f, Figure 5d, Figures 6b, 6d, and 6g, and Figure 7a is\nchallenging to read. Enhancing the clarity of these figures would greatly benefit the reader's\nunderstanding of the results.\nText in 3e is 5.4pt and cannot be enlarged further as the figure is full width.\nFigure 3f has been moved to supplementary figure 1 and figure indicators and main text adjusted\naccordingly.\nFigure 5d has been edited to increase the font size of the smaller text\nThe writing size has been increased as much as possible in Figures 6b, 6d, 6g\nFigure 7a has been amended to only show significant features and the features have been abbreviated\nand the font size increased. The full feature set has been moved to supplementary.\nReviewer #3 (Remarks to the Author)\nThis is an interesting work with specific target to reveal the spatial mechanisms of immunotherapy\nresponse in lung cancer based on MPX data. My detailed comments below:\nThe authors thank reviewer 3 for their valuable feedback.\n1) the overall writing is clear, but some parts are very wordy and convoluted. Please furthere proofread.\nWe have further proofread the manuscript and made extensive edits throughout the document to clarify\nthe language where appropriate. (in addition, in the legends of figures)\n2) The survival prediction approach first selects features based on binary labels, then filters them using\ncontinuous survival data via Cox regression. It would be interesting to compare this with a method where\nfeature selection is performed directly using time-to-event data, rather than relying on manual\ndichotomization.\nWe thank the reviewer for this very interesting suggestion. At the time of submission, the Stabl\nimplementation only supported binary comparison models. We have since developed an implementation\nof feature selection to model PFS directly instead of CB6 using CoxPH fits.\nThis tended to select more features than previously seen, so we combined several feature families to\nreduce the total number of selected features. Overall, out of the 88 features that we selected, 37 of them\nwere previously selected to model CB6, and we notice common themes in many of the remaining\nfeatures even if the selected features didn’t match exactly.\nThe performance boost of these selected features in modeling PFS was substantial, reaching\napproximately 90% across a 24 month period in k-fold fits, with additional gains in modeling OS. Using\nthese features to predict CB6 using logistic regression k-fold fits yielded a median ROC AUC of\napproximately 91%, demonstrating that fits to PFS likely yield higher predictive power as the\ndichotomization definition might induce noise in CB6, in addition to the reduced information available in\nthe binary response variable. We have included these results:\nLine 429-433:\nA supplementary model based directly on time to event (PFS) feature selection identified many similar\nfeatures (37/88), exhibiting common themes of metabolic pathways and interactions (Supplementary Fig\n3b-d) with similar predictive power (AUC 0.9). The binary model is shown here for consistency with\nprevious binary tests.\nSupplementary Figure 3b-c: Feature selection to model PFS. a) Hazard ratios of significant features\nthat were selected to model PFS using CoxPH fits during Stabl feature selection. Hazard ratios for OS are\nshown on the left and for PFS on the right. Features also selected to model CB6 are highlighted in bold.\nTime-dependent AUCs are shown for OS on the left and PFS on the right from k-fold fits. b) Prognosis\nscores formed by multiplication of feature hazard ratios with normalized feature values predict OS (left)\nand PFS (right). c) ROC curves using these features to model CB6 using logistic regression show excellent\nCB6 prediction power in k-fold fits.\n3) The exploration of spatial mechanisms in this work heavily depends on cell segmentation and\nclassification. However, technical details regarding the robustness and accuracy of these methods are\nlimited. Given the substantial variation and dynamic range typical of MPX images, more comprehensive\nevaluation of image processing quality is needed.\nWe appreciate the reviewer’s attention to the importance of a robust cell segmentation and classification\nmethod for these types of datasets. One major advantage of the commercial deep learning classifier\napplied here is that segmentation only provides a cell centroid, around which a tile is cropped for input\ninto the DL classifier, and so mean cell mask expression, which is often confounded by segmentation, is\nnot used in the classification of cells. We have improved conveyance of this in:\nLine 641-644\n“Cell segmentation provided cell centroids for downstream deep learning tile classification, however\nwas not used in cell classification. Subsequent spatial calculations were performed on cell centroids.”\nFurther notes on the classification pipeline:\nThe end-to-end pipeline has a cell-typing accuracy of over 91% and functionalization accuracy of over\n90%, with F1 scores on most cell-type classifications above 0.8, all tested on unseen data. It is a\nproprietary and state-of-the-art algorithm (https://doi.org/10.1101/2022.11.09.515776), and is trained on\nover 10,000 expert-annotated cells across a range of resolutions and bit depths. Issues associated with\nvariable dynamic ranges, which is typical of data such as these, are accounted for by a custom stain\nnormalization algorithm which both enhances the image signal-to-noise ratio and standardizes channel\nintensities across the cohort. The cell segmentation method is optimized for CODEX data, and like most\ncell segmentation algorithms, whole cell segmentation was based on the sum of all available membrane\nand cytoplasmic channels and DAPI, leveraging multiple markers to define cell boundaries using a\nnested convolutional neural network. Cell-typing was applied by supervised learning across thousands of\npathology-reviewed annotated patches, with accuracies far surpassing traditional expression-based\nclustering- approaches. Cells are labelled if their protein expression profiles matched those within the\ntraining data, reducing spurious classifications. This standardized algorithm for segmentation and cell\ntyping gives a robust depth of cell-typing and functional stratification uniquely suitable for our proposed\nmethod.\nWe show below the cell-typing accuracies stated for this pipeline, showing high accuracies across a\nrange of cell types that far surpass the cell-type accuracies found when employing traditional clustering-\nbased approaches.\nExamination of cell annotations on both our raw data and visualization of our heatmaps (now figure 1g)\nshowed that the algorithm worked very well for cell detection, cell typing, and functional classification.\nSince all downstream analysis components were dependent on the classifications and not the measured\nprotein expression within the segmentation masks, our analysis approach was relatively invariant to\nsegmentation inaccuracies, assuming that nuclei could be detected, which is not a substantial issue in\nCODEX data. Therefore, our study was designed to reduce any downstream propagation and\nreproducibility issues from the variability common in segmentation and cell-typing tasks, particularly for\nclustering-based cell-typing approaches, which is reflected in the translation of our spatial metrics\nacross two independent cohorts, with high prediction accuracies.\n4) The derivation of cell metabolism-related markers associated with survival and treatment response is\ncompelling. Nevertheless, these findings have not been validated on an independent cohort. Given the\nhigh-dimensional nature of MPX data, external validation is essential to confirm the reliability and\ngeneralizability of the results from the discovery cohort.\nWe thank the reviewer for their suggestion and agree that exhaustive cohort level validation is required to\nestablish clinically actionable biomarker targets. The work presented here aims to provide insights using\nstate of the art methods that may be expanded in future efforts. Despite the use of two cohorts on two\nTMAs generated independently of each other at separate timepoints over a nine year period from\ndifferent NSCLC patient cohorts, the available tissue limited the numbers required for such independent\ncohort level validation. We believe the target audience for this manuscript submission would understand\nthis limitation. We have made efforts in the manuscript to highlight this as well as the subsequent\nlimitations. Additionally, we have provided further exploratory data regarding consistency of cross cohort\nfeature selection, where we highlight features common to both independent cohorts as well as the\ncombined model.\nTo clarify this in the manuscript, we have added the following lines:\nLine 87:\n“Two TMAs were constructed from independent cohorts from NSCLC patient tumors presenting clinically\nbetween 2011 to 2017 (YTMA404) and 2017-2019 (YTMA471)”.\nResults lines\nLines 436:\n“Furthermore, we sought to define the cross-cohort predictive power by fitting each cohort to prediction\nmodels using both the CB6 and time to event PFS selection method (CB6 model; cohort 471 AUC 0.73,\ncohort 404 AUC 0.65, PFS model; cohort 471 AUC 0.87, cohort 404 AUC 0.93 )(Supplementary Table 4).\nWe additionally implemented a grouped bootstrap sampling method where only one cohort was utilized\nduring each subsampling fit performed during the feature selection. This method pulled more\nsignificant features by equally weighting the 404 cohort, which contained fewer patients but had lower average PFS\ntimes.\nNotably, seven of these features were common to all three feature selection models (combined, 404/471\nindependent) and were each significant by Kaplan-Meier tests, supporting generalizability of features\nacross two independent cohorts (Supplementary Table 4). Features that were predictive of benefit to ICI\nincluded the interaction between ICOS+ CD4 Tregs and fibroblasts in the tumor interface region,\nsuggesting Treg exclusion, as well as CD4 Tregs interacting with glycolysis+ tumor cells in stromal regions,\nimplicating a role for TME immune suppression of isolated but metabolically active tumor cells. Negatively\nassociated features implicated macrophages in several scenarios, where self aggregation of IDO1+\nmacrophages in low metabolic regions, and granzyme B+ macrophages in low metabolic regions were\nassociated with poorer ICI outcomes (Supplementary Table 4)”.\nDiscussion Lines 569:\n“Despite the use of two independently collected cohorts, a primary limitation of this study is the number\nof available tissues that passed quality control, allowing discovery and internal model cross validation\nonly. Thus, the work here provides a foundation for further rigorous biomarker discovery in spatial\nproteomics data that might evidence clinically actionable insights. While we were able to demonstrate the\nsuitability of this methodology across these two cohorts, feature selection applied using the cross-\nvalidation method can reduce sensitivity of selected features and increase the sensitivity of model\npredictions to cohort-specific effects, so we anticipate future studies utilizing larger cohorts across more\nthan two independently constructed TMAs to be able to apply feature cross-validation”.\nMethods Lines 603:\n“Two independently collected cohorts collected between 2011 and 2019”\nNotes to the reviewer:\nPublications which have used these study cohorts (DOI 10.1016/j.jtho.2022.04.009, DOI 10.1158/1078-\n0432.CCR-21-2649) acknowledge this as a study limitation for these types of exploratory work. While\nindependent cohort validation is useful to demonstrate signatures of response appropriate to the\nmaturity of the project, this is a discovery-stage cohort which provides groundwork to the scientific\ncommunity for rigorous analysis of datasets within spatial proteomics. Publication of novel\nmethodologies to analyze these datasets is essential to progress and optimize analysis pipelines.\nWhile clinical implementation is within the scope of our work, the requirement of independent validation\nfor discovery projects provides barriers to publication of early-stage discovery. There is no standardized\nworkflow for analysis of these types of data sets, and so we sought to dedicate resources in improving\nthis need. It presents an end-to-end pipeline for biomarker extraction and aims to provide a workflow for\nother labs to adopt to derive predictive spatially relevant biomarkers – a current unmet need in the field.\nThe study had appropriate statistical power to meet the discovery objectives, the study size is very\ncompetitive relative to similar impact studies, and the models were validated using k-fold fits to prove\ntheir translation to unseen data."}, {"round": 2, "reviewer_comments": {"reviewer_1": "1. Regarding the distance threshold selected for G-cross analysis, in the rebuttal letter, the authors indicated a threshold of 150um and this is contradicting to the statement \"Internal tests show that if the G-Cross function plateaus at a certain radius, for example 200μm, then the G-Cross AUC is typically the computed value up to 200um plus that additional “100 um”.\"\nIn the rebuttal letter, the authors indicated a threshold of 150um and this is contradicting to the statement \"Internal tests show that if the G-Cross function plateaus at a certain radius, for example 200μm, then the G-Cross AUC is typically the computed value up to 200um plus that additional “100 um”.\"\nIn addition, I don't see the edits in Line310: \"the distance unit reported in Line 310 is in pixels and this has been corrected to “150 μm,” to be consistent with the correct value/unit that is stated in the methods section.\"\nSimilarly, I cannot locate the new edits in Lines 346-402: \"We have clarified this in the paper by referencing “informative/uninformative features” in the “Clinical benefit spatial feature selection” section of the paper, Lines 346-402, and in the Fig.6 figure legend line 406.\"\n2. Regarding the predictive prognostic analysis, I cannot locate the following claim:\n\"Further cross validation predictive modelling of granzymeB+macrophages is provided later in the paper, where this feature within tumor regions is selected andpassed STABL feature selection criteria in both cohorts (Fig.3d-e)\"", "reviewer_2": "The authors have thoroughly addressed my concerns through detailed responses and constructive discussion. They have improved the quality of the figures and enhanced the reproducibility of their work by publicly sharing their analysis code and datasets, which will be valuable resources for the research community. In my opinion, the paper is now acceptable for publication.", "reviewer_3": "The author has addressed all my comments."}, "author_response": "REVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nWe thank the reviewer for their substantial effort to improve the presentation and depth of the\nmanuscript. We hope we are able to clarify any confusion in previous responses.\n1. Regarding the distance threshold selected for G-cross analysis, in the rebuttal letter, the authors\nindicated a threshold of 150um and this is contradicting to the statement \"Internal tests show that if\nthe G-Cross function plateaus at a certain radius, for example 200μm, then the G-Cross AUC is\ntypically the computed value up to 200um plus that additional “100 um”.\"\nWe apologise for confusion here. Our G-cross analysis summarises interactions on a scale of 150um,\nand we state this in both the manuscript, rebuttal and supporting supplementary figure 6b. Our further\nexample in the rebuttal to the reviewer is purely for descriptive purposes.\nWe have clarified this in methods to reference supp6b, lines 712:\nInternal testing indicated that 150um captured a majority of cell type proximity characteristics while\nrepresenting a biologically meaningful scale (Supplementary Fig. 6b)\nThe supporting figure provided in the rebuttal document is too small to be readable, and the\nsupposedly corresponding Supp Fig 6b is missing from the 'merged' file.\nThe plot for G-cross curves for all cell type pairs has been refined for visibility both here and in\nsupplementary 6b. High resolution vector graphics are provided for the journal.\nLines 848:\nSupp6 b) Globally observed G-Cross curves for cell type combinations sampled from 0 to 300\n(150um) pixels in increments of six pixels. Rows represent the query or \"from\" cell type, and columns\nrepresent the target or \"to\" cell type. Error bars indicate uncertainty bounds which represent ranges\nfrom all cores.\nIn addition, I don't see the edits in Line310: \"the distance unit reported in Line 310 is in pixels and this\nhas been corrected to “150 μm,” to be consistent with the correct value/unit that is stated in the\nmethods section.\"\nWe apologise for confusion here, as reference should have been to line 308 – original statement read:\nLine original text296-300:\nMore localized measures of the proximity of cells were captured by the “G-Cross” feature, which\nrepresents the area under the curve (AUC) of the nearest neighbour cumulative radial distribution of\ncross cell/cell plus functional/metabolic types, computed for radii up to 300 µm\nThis was pixel units, not um, and was corrected (as per reviewer’s point 11):\nLine new text 305-308:\nMore localized measures of the proximity of cells were captured by the “G-Cross” feature, which\nrepresents the area under the curve (AUC) of the nearest neighbor cumulative radial distribution of\ncell/cell or functional/metabolic types, computed up to a radius of 150 µm\nSimilarly, I cannot locate the new edits in Lines 346-402: \"We have clarified this in the paper by\nreferencing “informative/uninformative features” in the “Clinical benefit spatial feature selection”\nsection of the paper, Lines 346-402, and in the Fig.6 figure legend line 406.\"\nWe apologise for loose wording and confusion, this discussion with the reviewer regarding the meaning\nof FDR values references the paragraph headed : “Clinical benefit spatial feature selection”\nOriginal text lines 337:\nThe package allows for stable selection of relevant features while also estimating a bound on the false-\ndiscovery rate (FDR) by artificial feature injection (Supplementary Figure 1).\nLine new text 351:\nThe package allows for stable selection of relevant features while also estimating a bound on the false-\ndiscovery rate (FDR) of uninformative features by artificial feature injection (Supplementary Fig. 3).\nFig 6 legend, original text line 388:\nAn overview of the Stabl feature selection employed in this analysis, where feature families are tested\nfor stable selection across a range of pseudo-experiments that account for false discovery rates via\nnoise injection\nFig 6 legend,new text line 403:\nAn overview of the Stabl feature selection employed in this analysis, where feature families are tested\nfor stable feature selection across a range of pseudo-experiments that account for false discovery rates\nvia noise injection to derive informative/uninformative features\n2. Regarding the predictive prognostic analysis, I cannot locate the following claim:\n\"Further cross validation predictive modelling of granzymeB+macrophages is provided later in the\npaper, where this feature within tumor regions is selected and passed STABL feature selection criteria\nin both cohorts(Fig.3d-e)\"\nWe apologise for the mislabelling of where this statement occurs. In response to a separate reviewer,\nwe generated models for each ytma404/471 cohort independently, as well as for combined cohorts,\nand one such feature selected within both cohorts was the proportion of granzyme B+ macrophages in\nmetabolic low regions, which coincides with their significance in tumour regions in univariate analysis.\nWe acknowledge the reviewers concerns about independently prognostic features in the absence of\nlarger cohorts, and discuss this as a limitation.\nEdits to results regarding feature generalisability, lines 434-452:\nFurthermore, we sought to define the cross-cohort predictive power by fitting each cohort to prediction\nmodels using both the CB6 and time to event PFS selection method (CB6 model; cohort 471 AUC 0.73,\ncohort 404 AUC 0.65, PFS model; cohort 471 AUC 0.87, cohort 404 AUC 0.93 ) (Supplementary Fig\n3e). We additionally implemented a grouped bootstrap sampling method where only one cohort was\nutilized during each subsampling fit performed during the feature selection. This method pulled more\nsignificant features by equally weighting the YTMA 404 cohort, which contained fewer patients but had\nlower average PFS times.\nNotably, seven of these features were common to all three feature selection models (combined, 404/471\nindependent) and were each significant by Kaplan-Meier tests, supporting generalizability of features\nacross two independent cohorts (Supplementary Fig 3e). Features that were predictive of benefit to ICI\nincluded the interaction between ICOS+ CD4 Tregs and fibroblasts in the tumor interface region,\nsuggesting Treg exclusion, as well as CD4 Tregs interacting with glycolysis+ tumor cells in stromal\nregions, implicating a role for TME immune suppression of isolated tumor cells. Negatively associated\nfeatures implicated macrophages in several scenarios, where self-aggregation of IDO1+ macrophages\nin low metabolic regions, and granzyme B+ macrophages in low metabolic regions were associated with\npoorer ICI outcomes (Supplementary Fig 3e).\nDiscussion cohort limitations, lines 566-573:\nDespite the use of two independently collected cohorts, a limitation of this study was the available\ntissues that passed quality control, allowing discovery and internal model k-fold validation, with some\ncross cohort consistency in selected features. While we were able to demonstrate the suitability of this\nmethodology across these two cohorts, feature selection applied using the cross-validation method can\nreduce sensitivity of selected features and increase the sensitivity of model predictions to cohort-specific\neffects, so we anticipate future studies utilizing larger cohorts to be able to apply robust feature cross-\nvalidation."}, {"round": 3, "reviewer_comments": {"reviewer_1": "I appreciate the authors’ efforts to address the comments. I am satisfied with the revised version.", "reviewer_2": "The authors have thoroughly and satisfactorily addressed all the points I raised in my previous review. The manuscript is now suitable for publication."}, "author_response": ""}]

Decision letter