Introduction

In recent years, immunotherapy has emerged as the key promise to revolutionize cancer treatment1, delivering higher survival rates and lower toxicity than standard treatments. Although immunotherapy demonstrated significant potential in experimental models2, translation to clinical practice has shown a lesser impact. Traditional treatment approaches for head and neck cancer have centered on surgery, followed by risk-adapted radio/chemoradiotherapy3,4 or upfront chemoradiotherapy. The treatment landscape began to expand in 2016, with the FDA approval of anti-programmed cell death 1 (PD-1) immune checkpoint inhibitors, pembrolizumab and nivolumab, as second-line therapies. In 2019, pembrolizumab was approved as a first-line treatment for patients with metastatic or unresectable recurrent HNSCC. Despite these advances, clinical trials of anti-PD-1 monotherapy in the recurrent-metastatic setting have yielded only modest survival benefits, although they demonstrated a more favorable safety profile compared to standard chemotherapy5,6,7,8. These trials underscore the challenges of treating recurrent-metastatic HNSCC, where single-agent therapies often fail to address the diverse mechanisms of resistance4,9.

In this work, we investigate the relationship between the molecular composition of HNSCC tumors and the favorable response of HNSCC to anti-PD-1 therapy, either as monotherapy or in conjunction with glucocorticoid-induced tumor necrosis factor receptor-related protein (GITR) to enhance HNSCC diagnosis and treatment. GITR is a potential target in cancer immunotherapy, primarily due to its dual role in enhancing effector T-cell functions and attenuating the suppressive effects of regulatory T-cells (Tregs). GITR is expressed predominantly on Tregs and, to a lesser extent, on effector T-cells. Upon ligation, GITR lowers the activation threshold of CD8+ T-cells by enhancing CD28-mediated costimulation, which increases T-cell proliferation, survival, and cytokine production, notably of IL-2 and IFNg10. Additionally, agonistic anti-GITR antibodies have demonstrated the ability to provoke potent tumor-specific immune responses in preclinical models. These antibodies not only stimulate CD8+ T-cells but also reduce the suppressive activity of Tregs within the tumor microenvironment, facilitating enhanced tumor infiltration by effector T-cells and improving the efficacy of immunotherapy in advanced tumors11. Despite promising preclinical results, clinical trials involving patients with advanced solid tumors or advanced cancer have revealed that GITR monotherapy has limited antitumor efficacy. For example, early-phase trials with the GITR agonist TRX518 in patients with advanced cancer showed that, while Treg cells were reduced and effector T-cell to Treg ratios improved, these changes did not translate into substantial clinical responses. This has led to the exploration of combination therapies, particularly with PD-1/PD-L1 inhibitors. Combining GITR agonism with PD-1 blockade has shown synergistic effects in preclinical models, potentially overcoming resistance to monotherapies by reinvigorating exhausted T-cells and enhancing tumor cell killing12,13. However, even in combination therapies, the clinical activity of GITR agonists in patients with advanced solid tumors has so far been modest, with objective response rates similar to those observed with PD-1 monotherapy14. As such, GITR-targeted therapies, as a single agent or in combination, require further exploration to fully realize their potential in cancer immunotherapy15.

One problem in the clinical translation of immunotherapies is the wide variation in the immune environment of patient tumors. Biomarker expression levels, such as PD-1/PD-L1 and cytotoxic T-lymphocyte associated protein 4, have been used to predict immune response, however, with limited success in HNSCC patients16. This indicates that it is unlikely for any singular biomarker to comprehensively predict the efficacy of immunotherapy, given the intricacy of the immune response4,17. Therefore, there is a pressing need to develop a multi-index approach that integrates a broader understanding of the dynamic processes in HNSCC, rather than focusing on isolated biomarkers. We propose that this multi-index strategy will offer a more accurate prediction of HNSCC tumor response to immunotherapy17.

In this study, we employed explant models profiled using secretomes, HNSCC datasets available online, and computational information-theoretic analysis of HNSCC secretomes/transcriptomes to develop a predictive strategy for the immune response of HNSCC patients (Fig. 1A, B). We discovered new co-expression patterns of transcripts and cytokines associated with enhanced T-cell infiltration and immune-activation factors, indicating a more favorable response to immunotherapy. We verified our findings experimentally in patient-derived explants and two separate big datasets and demonstrated the potential of the co-expression patterns in predicting HNSCC outcomes following treatment with immunotherapy (Fig. 1C). This study emphasizes the importance of multiple-index characterization of HNSCC tissues in improving patient classification and better predicting immunotherapy efficacy.

Fig. 1: Overview of the study workflow.
figure 1

A Secretome data were collected from patient-derived fresh tumor explants, and RNAseq was obtained from two large HNSCC cohort studies. Multi-parameter immunofluorescence (mIF) of patient-derived explants was used to validate the results experimentally. B Surprisal analysis (SA) was carried out on the secretomes of explant supernatant and the RNA sequencing of the two large cohorts (Methods). In brief, proteins/genes whose expression levels deviate from the reference state in the same or opposing manners (upper right panel) are computed via the covariance matrix of the natural logarithms of the expression levels as dictated by the theory (middle panels and “Methods”). Transcripts/proteins that deviate from the steady state in a coordinated manner are grouped into unbalanced processes (lower panel). Several different protein–protein (mRNA-mRNA) unbalanced processes can be found in the entire population of cancer samples. C Patient-derived explants and available datasets of HNSCC patients were used to confirm the relationship between unbalanced processes associated with T-cell function and response to immunotherapies (IM) and patient survival.

Full size image

Results

Secretome analysis of HNSCC patient-derived tissues reveals predictive unbalanced processes strongly linked to enhanced T-cell infiltration and immunological response

Secretome profiles from 72 explant samples originating from six individual patients (18 untreated plus 54 treated, Supplementary Data 1, “Methods”) were investigated using surprisal analysis (SA). In brief, SA is a thermodynamic-like, information-theoretic approach that computes one or more co-expression protein–protein (transcript–transcript) altered patterns that reflect environmental and genetic constraints in a system under investigation. These constraints prevent a system (i.e., tissue) from achieving a steady state and are frequently associated with different tumor characteristics18,19,20. Since each tumor can be characterized by several characteristics (=phenotypes), such as migratory and highly proliferative phenotypes or immune-responsive and hypoxic features, SA calculates a set of tumor-specific constraints (=phenotypes) using measurements obtained in omics approaches. Proteins (or transcripts) that deviate from steady state in a coordinated manner are grouped into unbalanced processes that represent tumor-specific phenotypes at the molecular level18,19. SA assigns a set of unbalanced processes (=phenotypes) to each tumor based on the amplitudes (\({\lambda }_{\alpha }\left(k\right)\), importance) computed for each sample (k) (Methods).

SA examination of the secretome data revealed that the two most dominant unbalanced processes (responsible for the largest variability in the data, “Methods”) were associated with immune response. In process α = 1, immune-related proteins such as PD-L1, GZMH, and IFNg were induced and co-expressed in the patients labeled in blue (Fig. 2A, B). In contrast, proteins related to angiogenesis (TNFRSF12A, HMOX1, VEGFA) were decreased (Supplementary Data 1 includes all computed parameters as well as a detailed look at the biological meaning of each process based on gene ontology (GO) analysis of the co-expressed proteins in each process).

Fig. 2: Bulk SA analysis of secreted proteins in HNSCC tissues using thermodynamic methodology revealed unbalanced processes in HNSCC that are closely linked to T-cell infiltration and improved immunotherapy response.
figure 2

A, C Amplitudes \(({\lambda }_{\alpha }\left(k\right))\) of process 1 (Act+) and 2 (Inf+) were computed (Methods) for each patient (n = 3 samples per patient) and presented for each tissue. Patients 4698 and 4950 harbor both processes, while 4905 only Act+. B, E Cytokines with significant weights (Gi, “Methods”) were co-expressed and grouped into unbalanced processes using Gi values and STRING protein–protein connections. The circle radius corresponds to the Gi value of each protein i. The sign of \({\lambda }_{\alpha }\left(k\right)\) (A, C) determines correlation/anti-correlation, or separation, between the samples (Methods): Blue-labeled proteins are upregulated in patients with blue-labeled amplitudes \(\left({\lambda }_{\alpha }\left(k\right)\right)\), and red-labeled proteins are upregulated in patients with red-labeled amplitudes. D Inf phenotype differentiates between HPV+ and HPV− samples (T-test **p value < 0.01).

Full size image

A combination of proteins linked with this process is associated with immune system activation, thus, the process was termed the Act phenotype. Another process (process α = 2) involved co-expressed T-cell biomarkers such as CD8 and CD4, as well as decreased ligands involved in paracrine pro-survival cancer phenotype (EGF, PDGF). This combination suggests an increased infiltration of T-cells in the tumor independent of activation; thus, the process was termed the Inf phenotype (Fig. 2C, E). Not all examined tissues harbored Act+ or Inf+ processes. Patient tumors with IDs 4698 and 4950 included both Act+ and Inf+ processes (Fig. 2A, C), whereas tissue 4856 did not exhibit any of these phenotypes due to the opposite expression of biomarkers associated with Act or Inf phenotypes. Interestingly, process α = 2 could differentiate between HPV+ and HPV- patients (Fig. 2D), suggesting that Inf phenotype is associated with the distinct immunobiology of HPV+ tumors21.

Next, we investigated a potential relationship between the computationally discovered cytokine-related unbalanced processes, the experimental assessment of T-cell proximity to the tumor cells, and the ex vivo response of the HNSCC explants to immunotherapies. To assess infiltration, we employed quantitative multi-parameter immunofluorescence (mIF) labeling on tissue blocks from the same patients. The mIF panel incorporated CD3, CD8, PD-L1, and CK, as previously described22,23. This allowed us to identify infiltrating CD3 + CD8+ T-cells by the infiltration and immune-activation phenotypes. We performed a quantitative assessment of the number of CD8+ T-cell infiltrating cells per tumor, and identified a direct correlation between the computationally-derived Inf process and the number of infiltrating CD8+ T-cells. Interestingly, the number of infiltrating T lymphocytes was higher in HNSCC with both processes than in tumors with only the Act process (patient 4905) or without any process present. The number of infiltrating T lymphocytes was also highly correlated with the amplitude of the Inf phenotype (Fig. 3A, B).

Fig. 3: Multi-parameter immunofluorescence (mIF) of the patient-derived samples reveal that increased infiltration of T-cells and response to immunotherapy can be predicted by SA based Act+/Inf+ phenotypes.
figure 3

A CD8 and PD-L1 measurements using mIF in patients 4698 (Act+/Inf+) and 4321 (Act−/Inf−). B Average # of CD8+ T-cells per sample was plotted vs. the amplitude (\({\lambda }_{2}\left(k\right)\)) of Inf+ phenotype (process 2) of each sample k. (Pearson r correlation). C Schematic view of cell-proximity analysis performed using the R “phenoptr” package on the mIF images. CD8+ T-cells within the separation distance of 10, 20, 30, or 40 µm between the T-cells and PD-L1+ cancer cells were counted. D Quantification of CD8+ T-cells as shown in (C) (2-way ANOVA test ****p value < 0.0001). E Comparison of cytokine levels between individuals with activation-phenotype and patients without the phenotype (T-test P value *<0.05, **<0.01, *** <0.001, ****<0.0001). F “Before and after” treatment plots present the levels of cytokines associated with response to immunotherapy. For each patient, fold change of cytokine levels after immunotherapy (IM, anti-PD-1 gray, anti-GITR red, Combination orange) compared to the mean level of untreated explants (ctrl) was calculated. Patient with Act+ phenotype (4905) was compared to the Act+/Inf+ patients (4950 and 4698). (Wilcoxon signed-rank test and Mann–Whitney test were used, p value *<0.05, ***<0.001, ****<0.0001).

Full size image

Given that PD-L1 was a feature of the Act process, along with IFNg, we hypothesized that PD-L1 was upregulated in response to the proximity of tumor-specific CD8+ T-cells secreting IFNg. To test this, we used spatial analysis of the distance between PD-L1+ cancer cells (CK + PD-L1+) and CD8+ T-cells in the tumors. We found that in patients with both Act+ and Inf+ phenotypes, more CD8+ T-cells were clustered close to PD-L1+ cancer cells than in patients without both phenotypes (Fig. 3C, D). Thus, the combination of both Act+ and Inf+ phenotypes may indicate that not only are T-cells present in tumors, but they are also more likely to interact with the cancer cells. This is notable in the mIF specimen, where cancer cell PD-L1 expression is strongest at the stromal boundary, where these cells are more likely to interact with infiltrating CD8+ T-cells (Fig. 3A).

To assess the response of these tumors to immunotherapy, matched explants were treated with anti-GITR, anti-PD-1, or their combination. Tumor supernatants were harvested as in untreated tumors and quantified for secreted markers. Importantly, we noted that the response to immunotherapy was different between samples. Although HNSCC tissues with the Act+ phenotype had induced immune-response biomarkers and reduced signaling pathways, such as VEGFA (Fig. 3E), known as immune inhibitors24, only samples with both phenotypes (Act+/Inf+) responded significantly better to immunotherapy (IM: anti-PD-1, anti-GITR (activating T-cells), or anti-PD-1 + GITR combination) (Fig. 3F). This was demonstrated by a significant increase in biomarkers indicating a favorable response to immunotherapy, such as IFNg, GZMH, and FASL, and decreased levels of IL6, a negative immunomodulator (Fig. 3F).

These data imply that: (1) Immune response is not well predicted by one (or even a few) biomarkers; even a process that includes a collection of co-expressing biomarkers, such as process 1, when employed alone, has a lower predictive power; (2) HNSCC’s immune response and T-cell infiltration can be more precisely defined and predicted by complex co-expression patterns, such as process combinations (like the two processes found here); (3) Secretome data collected on explant platforms and evaluated by SA has the potential to be a successful technique for diagnostics and immune-response prediction.

Act and Inf phenotypes are found in a large independent cohort of HNSCC patients

To further evaluate the significance of the Act and Inf phenotypes discovered in the patient-derived explants, we analyzed transcriptomes from 518 primary HNSCC tumors retrieved from the TCGA database.

Using SA, we identified the Act phenotype, characterized by induced IFNG/GZMH and PD-L1 co-expression, in 74 patients and the Inf phenotype (characterized by CD8A/CD4 induced co-expression) in 116 patients (Fig. 4A, B, “Methods”). Additional co-expressed transcripts in Act and Inf phenotypes can be found in Supplementary Data 2. Both phenotypes were found in 21 tumors (Supplementary Data 2). Interestingly, individuals with both phenotypes showed higher levels of immune biomarkers (Fig. 4C), indicating an improved potential for response to immunotherapies (Fig. 3). Furthermore, decreased co-expressed transcripts found in Act and Inf phenotypes (Fig. 4D, E, Supplementary Data 2) were associated with immunoinhibitory signaling pathways, such as MAPK, NFKB, angiogenesis (Inf phenotype) or PI3K, angiogenesis, and MAPK (Act phenotype), similar to the results of the secretome analysis of the patient-derived explants (Fig. 2).

Fig. 4: Act and Inf phenotypes are found in 518 HNSCC patients from the TCGA dataset.
figure 4

A A two-dimensional graphic showing the Act and Inf phenotypic amplitudes (λ) for every patient under examination. Patients were classified as having those phenotypes if their amplitudes were greater than the error limits. The amplitudes’ signs indicate a spatial separation (Methods). Negative amplitudes are associated with induced Act and Inf phenotypes. B Selected co-expressed transcripts in the Act phenotype and Inf phenotype, respectively, are shown. A full list of additional transcripts can be found in Supplementary Data 2. C Expression levels of immune biomarkers IFNG, GZMH, and PD-L1 in patients classified as having one phenotype only (Act or Inf) or both phenotypes simultaneously (Act and Inf) (One-way ANOVA, p value *<0.05, ***<0.001, ****<0.0001). D, E Gene ontology enrichment for downregulated (red) transcripts in Act and Inf phenotypes, respectively. F Amplitudes (λ) of the immune-infiltration phenotype were determined and compared between p16+ (indicating HPV positive status) and p16− patients (total 102) (T-test ****p value < 0.0001). CD8A gene expression levels are shown in subgroups of individuals classified by HPV status (p16) (One-way ANOVA, p value *<0.05, **<0.01, ****<0.0001). Percentages represent the number of patients with or without the Inf phenotype.

Full size image

Likewise, the HPV+ HNSCC subgroup, as identified by the p16 biomarker25, had the Inf phenotype (Fig. 4F), as did the patient-derived explants with the Inf phenotype (Fig. 2). Our findings in HNSCC explants are validated by these data.

Higher survival rates for HNSCC patients who underwent immunotherapy correlate with the combined expression of the Act and Inf phenotypes

To further investigate the predictive potential of the Act and Inf phenotypes, we examined a dataset of 102 patient samples with advanced HNSCC (GSE159067, GEO database) treated with immunotherapy targeting PD-1/PD-L1 (as part of the clinical trial)26. The overall survival of patients treated with PD-1/PD-L1 inhibitors26 was related to the processes identified by SA analysis of RNAseq data. SA analysis was carried out on the samples that were collected prior to therapy (Supplementary Data 3).

The presence of Act and Inf phenotypes in HNSCC is also strongly supported by the SA results of this cohort (Fig. 5A, B). Similar to the other HNSCC datasets in this study, patients with both phenotypes have demonstrated a decrease in pro-survival cancer characteristics (Fig. 5C, Supplementary Data 4). Furthermore, clinical results clearly show that anti-PD-1/PD-L1 inhibitor treatment increases the survival rate of Act phenotypic patients. Individuals within this subgroup who had both the Act and Inf phenotypes demonstrated an even greater correlation with improved survival (Fig. 6).

Fig. 5: Pro-survival networks decreased in patients with Act and Inf phenotypes (GSE159067 dataset).
figure 5

A, B 2D plot depicting the amplitudes (λ) of the Act and Inf phenotypes in all examined patients and the selected co-expressed transcripts in the Act phenotype and Inf phenotype, respectively, are shown. The amplitudes’ signs indicate a spatial separation (Methods). Negative amplitudes are associated with induced Act and Inf phenotypes. The full list of additional transcripts can be found in Supplementary Data 3. C The transcripts identified in the Act and Inf phenotypes were used to perform the single-sample gene set enrichment analysis (ssGSEA, “Methods”). Among patients with both Inf+/Act+ or Inf−/Act− four biological categories (e.g., angiogenesis, EMT, hypoxia, and IL6/JAL/STAT3 signaling) enrichment scores are displayed (T-test, P value *<0.05, **<0.01, ****<0.0001).

Full size image
Fig. 6: Immunotherapy increases the survival rate of patients with Act and Inf phenotypes.
figure 6

Pearson-correlation plots between overall survival (months) post-treatment (with anti-PD-1/PD-L1 inhibitors) and the amplitudes of Act phenotype as computed for each patient (left panel). Increased Act phenotype (=more negative values) is associated with higher survival. The subpopulation of superior survivors—which is outlined in the left panel—was plotted further against the Inf phenotypic amplitudes (right panel).

Full size image

These findings validate the predictive power of the Act and Inf phenotypes and imply that the appearance of both should be a reliable indicator of how well immunotherapy will work for patients with HNSCC.

Discussion

Given the heterogeneous nature of HNSCC and the varied responses to immunotherapy, it remains crucial to develop strategies that can more accurately predict patient outcomes. A key objective is to characterize the immune microenvironment and identify resistance mechanisms at an individual patient level, thereby enabling the prediction of treatment efficacy before it is administered. Previously, we developed an explant system designed to predict treatment responses within a short timeframe27,28,29. This system demonstrated a correlation between the secretome profiles of explants and the immune environment in several HNSCC in vivo models, and in two human-derived samples27. Our current study expanded upon this methodology, focusing on combining anti-PD-1 and anti-GITR therapies. We aimed to correlate explant secretome patterns with mIF, TCGA, and GEO datasets from HNSCC patients.

Our study began by analyzing the secretomes of 6 patient-derived HNSCC explants (72 untreated/ treated samples), which were treated with anti-PD-1, anti-GITR, or a combination of these agents. Through SA, a thermodynamic-based information-theoretic approach, we identified two key unbalanced processes, the Act and Inf phenotypes, which were strongly linked to immune response and T-cell infiltration.

Although the discovery cohort is comprised of six patients, each explant profile contained 96 cytokines (>6900 measurements in total), providing a data-rich matrix for pattern extraction. Nevertheless, we consider these findings exploratory and avoid inferring mechanistic causality from this cohort.

These key unbalanced processes were identified despite the small number of patients in the study. The amplitudes we compute with the unbalanced processes provide an inherent measure of their importance, and this, combined with the multi-parameter analysis, can identify important linked processes with relatively small inputs. While we validate these linked processes on much larger patient populations, it is likely that the approach can be improved when performed on larger patient populations. While we show that we can successfully identify linked processes that incorporate known markers with predictive value in our small cohort, we believe that larger samples have the potential to discover additional novel associations and, importantly, identify a broad range of functional subpopulations to permit more nuanced profiling of patients than is currently available.

The Act phenotype was characterized by the co-expression of immune-activation markers, including IFNg, GZMH, and PD-L1, and the downregulation of angiogenesis-related proteins such as VEGFA. The Act phenotype suggests an environment primed for immune activation, essential for effective immunotherapy. The presence of this phenotype in patient tissues indicated a heightened potential for immune-mediated tumor regression.

The Inf phenotype was defined by the co-expression of T-cell markers CD8 and CD4 and the downregulation of ligands involved in cancer cell survival, such as EGF and PDGF. This phenotype was indicative of a higher capacity for T-cell infiltration into the tumor microenvironment, a critical factor for the success of immunotherapy. Notably, the Inf phenotype was more common in HPV-positive tumors, a finding consistent with clinical analyses showing superior objective response and overall survival to PD-1/PD-L1 blockade in HPV-positive HNSCC (HR 0.71, OR 1.79)30 and with immune-profiling data demonstrating greater B and T lymphocytes infiltration in the tumor microenvironment of HPV-driven oropharyngeal SCC31.

We emphasize, however, that this association does not demonstrate a direct causal link between HPV status and the Inf phenotype. Larger prospective cohorts will be required to test that relationship further.

Next, we examined whether there was an association among the computationally identified cytokine-related unbalanced processes, the experimental evaluation of T-cell proximity to tumor cells, and the ex vivo response of HNSCC explants to immunotherapies. The IHC data supported our computational predictions, showing that tissues with both the Act and Inf phenotypes exhibited higher levels of CD8+ T-cell infiltration compared to those with only the Act phenotype or neither phenotype. Moreover, spatial analysis revealed that CD8+ T-cells in Act+/Inf+ tissues were more closely clustered around PD-L1+ cancer cells, particularly at the stromal boundaries. This proximity suggests an enhanced potential for immune cell-cancer cell interactions, which is crucial for the effectiveness of immunotherapy.

To further substantiate the role of the Act and Inf phenotypes, we analyzed transcriptomes from 518 HNSCC tumors in the TCGA database. SA identified the same Act and Inf phenotypes in this independent cohort, with 74 patients exhibiting the Act phenotype, 116 the Inf phenotype, and 21 patients showing both. Patients with both phenotypes showed significantly higher levels of immune biomarkers, such as IFNg, GZMH, and PD-L1, consistent with an improved potential for response to immunotherapy. Additionally, GO analysis revealed that downregulated transcripts in both phenotypes were associated with immunoinhibitory signaling pathways, mirroring the results from the patient-derived explants.

Interestingly, the Inf phenotype was again closely associated with HPV-positive tumors, aligning with previous findings that HPV status may predict better responses to treatment.

The replication of Act and Inf across the TCGA database and, as described below, GSE159067 (combined n = 620) strengthens their general relevance. Yet, we interpret these cross-sectional data as associative rather than definitively causal.

The final validation of the Act and Inf phenotypes came from analyzing a cohort of 102 advanced HNSCC patients treated with PD-1/PD-L1 inhibitors. We found a strong correlation between the presence of these phenotypes and improved overall survival. Patients who exhibited both Act and Inf phenotypes before therapy had significantly better survival outcomes than those with either phenotype alone or neither. This finding suggests that the presence of both phenotypes serves as a predictor of favorable response to immunotherapy, highlighting the potential of these phenotypes as biomarkers for guiding treatment decisions.

Moreover, the combination of these phenotypes correlated with a reduction in pro-survival cancer characteristics, such as angiogenesis and immune inhibitory signaling, further supporting their role in enhancing the efficacy of immunotherapy.

The findings from this study underscore the complexity of the immune landscape in HNSCC and the limitations of relying on single biomarkers to predict treatment outcomes. The identification of the Act and Inf phenotypes provides a more nuanced understanding of the immune environment and its role in influencing the response to immunotherapy. By integrating these phenotypes into a multi-index predictive model, clinicians may be able to stratify patients better and tailor immunotherapy treatments to those most likely to benefit.

These data illustrate the efficacy of a novel, unbiased analytical approach for identifying multi-index predictive models with prospective utility. This work centered on the application of PD-1 and GITR-targeting immunotherapies, hence, we concentrated on associated markers pertaining to CD8+ T-cell immunobiology. This encompasses characteristics of the Inf phenotype that have established single marker associations with outcomes, such as CD8+ T-cell infiltration. The strength of our approach lies in the identification of additional biological categories within the processes, such as pro-survival and angiogenesis, which are less defined yet also classify patient malignancies. Furthermore, the found phenotypes may manifest individually or in combination, and this combination of processes enhances the predictive accuracy of responses to immunotherapy. The data indicate that although our methodology can detect recognized, verified biomarkers, there is a range of novel characteristics with unknown mechanisms and similar potentials that have not been identified as singular marker associations in previous investigations. We present an approach for multi-parameter profiling of novel gene/protein associations to classify responsive and unresponsive patient cancers. At the same time, we acknowledge that future research will be crucial to validate causal links, enhance the predictive model, and reveal further unbalanced processes.

Future research should focus on further validating the Act and Inf phenotypes across larger and more diverse patient cohorts. Additionally, exploring the mechanisms underlying these and other phenotypes could provide insights into how they might be modulated to enhance immunotherapy efficacy. The development of clinical tools that incorporate these phenotypes could ultimately lead to more personalized and effective treatment strategies for HNSCC patients, potentially improving survival outcomes in HNSCC.

Methods

Secretome measurements

Secretome profiles from 72 explant samples originating from six individual patients (18 untreated plus 54 treated, Supplementary Data 1) grown as fresh tumor explants27 immediately post-resection from HNSCC patients were obtained using Olink Target 96 Immuno-Oncology Panels. Secreted proteins were quantified in explants left untreated, treated with anti-PD-1, anti-GITR, or the combination (anti-PD-1+anti-GITR). Basal secretome networks in each tissue were investigated using SA (see below).

Sample collection from this study was approved by the institutional review board at Providence Portland 440 Medical Center, Portland, OR, with the study ID number PHS: 06-108. All patients provided informed 441 consent for participation in this study. This research was performed in accordance with the Declaration of Helsinki.

Explant analysis

Fresh tumor samples were delivered on ice to the laboratory, and processing of the tissue into explants began immediately. The method is as previously described27,28. Briefly, the tumor was dissected into cuboidal 1.5–2.0 mm fragments, and each fragment was placed into 1 well of a 96-well plate with media alone or media containing an immunotherapeutic agent. Tumor fragments were left untreated, treated with 10 µg/ml anti-GITR (INCAGN Incyte), 25 µg/ml anti-PD-1 (MGA12 Incyte), or the combination of anti-GITR and anti-PD-1. Three replicate fragments per treatment condition were collected. After 24 h of incubation at 37 °C, the supernatant media were collected and stored at −80 °C for later bulk analysis.

Multi-parameter immunofluorescence (mIF) analysis of infiltrating immune cells

mIF analysis was performed as recently described32. Briefly, human FFPE tissue blocks were cut into 4 μm thin sections, deparaffinized, and rehydrated. Antigen retrieval was performed using 1X Tris-EDTA (10 mM Tris Base, 1 mM EDTA, pH 9.0). Stripping was performed using 0.1 M Sodium Citrate pH 6.0, in a microwave oven. Endogenous peroxidase was blocked by 3% H2O2 for 15 min at room temperature (RT). Tissue sections were blocked with blocking/antibody diluent (ARD1001EA, PerkinElmer/Akoya) for 10 min at RT. For Opal fluorescent staining, tissue sections were incubated with primary antibodies at RT in a shaking stain tray. Incubation times for anti-PD-L1, anti-CD8, anti-CK, and anti-CD3 were 30, 30, 20, and 15 min, respectively. After a brief wash, tissue sections were incubated with MACH2 Rb HRP-Polymer (RHRP520H, Biocare Medical) or MACH2 M HRP-Polymer (MHRP520H, Biocare Medical) for 10 min at RT. A quick wash in 1X TBST was followed by incubation with Opal690 (1:200), Opal650 (1:200), Opal570 (1:400), and Opal540 (1:200), respectively, for 10 min at RT. Counter stain was done with DAPI (1 drop of DAPI solution into 0.5 ml of TBST, PerkinElmer/Akoya) for 5 min at RT. After a quick wash, the slides were mounted with Prolong Diamond Antifade Mountant (p36970, Thermofisher). p16 staining occurred as part of standard diagnostic pathology.

Imaging was performed on a Vectra 3 automated quantitative pathology system as recently described32. A 10× objective lens was used for whole slide scans, while a 20× objective lens was used for the multispectral images. To obtain images for establishing the spectral library, snapshots were taken from each single biomarker-stained slide. The snapshots for DAPI-stained and unstained slides were taken through the DAPI filter cube. Ten to fifteen multispectral images/ROIs were taken from each patient slide. The spectral library was generated using the snapshots of the single biomarker-stained slides. The spectra were selected based on the fluorophores. Composite images were generated by extracting unmixed signals using the spectral library. Snapshots of unstained slides were used to remove the background caused by tissue fixation. Quantification analysis was performed by following the standard procedure, which includes Tissue Segmentation, Cell Segmentation, Phenotyping, Scoring, and Export.

Cell-proximity analysis

The analysis was performed using R v4.2.2. Cell segmentation files were combined into a single aggregated cell segmentation file. Cell phenotyping was done by selecting for certain positive marker combinations. CK+ or CK+/PD-L1+ cells were considered cancer cells. CD8+, CD8+/PD-L1+, CD3+/CD8+, or CD3+/CD8+/PD-L1+ cells were considered CD8+ T-cells. Aggregated cell segmentation files were read into R via the phenoptr package (v 0.3.2, https://akoyabio.github.io/phenoptr/). CK+/PD-L1+ cells within the tumor were designated PD-L1+ cancer cells, and CK+ cancer cells were designated PD-L1− cancer cells. CD8+ T-cells within 10, 20, 30, or 40 µm of PD-L1± cancer cells were counted using the “count_within” function in phenoptr. Dunn’s test from the rstatix package (v 0.7.2, https://rpkgs.datanovia.com/rstatix/.) was used as a non-parametric method to test for differences in means between using PD-L1− tumor cells as centroids and PD-L1+ tumor cells as centroids in the “count_within” function.

Surprisal analysis

A thermodynamic-like information-theoretic approach was applied in this study, as detailed in refs. 18,33,34. Briefly, the analysis is based on the premise that biological systems reach a balanced state when the system is free of constraints18,33. However, when under the influence of environmental and genomic constraints, the system is prevented from reaching the steady state. Each constraint can induce a change in a specific part of the protein/transcript network in the cells. The subnetwork that is altered due to the specific constraint is termed an unbalanced process. The system can be influenced by several constraints, thus leading to the emergence of several unbalanced processes (phenotypes) in each tumor.

Technically, the experimental levels of each protein/transcript are decomposed into steady-state values (values unaffected by constraints) and deviations due to constraints:

$$ln\,{X}_{i}\,(k)={\rm{l}}{\rm{n}}{X}_{io}(k)-\sum _{\alpha =1}{G}_{i\alpha }{\lambda }_{\alpha }(k)$$
(1)

where i is the protein/transcript of interest, \(ln{X}_{{io}}\left(k\right)\) is the protein/transcript level in a sample k when the sample is at the steady state free of constraints. \({\sum }_{\alpha =1}{G}_{i\alpha }{\lambda }_{\alpha }\left(k\right)\) reflects the total of the different constraints’ effects on the protein i’s expression level deviations18,33,35.

The term \({G}_{i\alpha }\) denotes the degree of participation of the protein/transcript i in the unbalanced process α. Proteins/transcripts with significant \({G}_{i\alpha }\) values (“G values”; Supplementary Data 1–3) are grouped into unbalanced processes (Figs. 1, 2B, E, 4B and 5A, B) which are active in the dataset. The term \({\lambda }_{\alpha }\) represents the importance of the unbalanced process α in the tumor k. The detailed description of how \({G}_{i\alpha }\) and \({\lambda }_{\alpha }\) values are calculated, including calculation of error limits, can be found in the refs. 18,33,34.

In all datasets, the two most dominant unbalanced processes, e.g., with the largest \({\lambda }_{\alpha }\) values, were selected (“Lambda values” Supplementary Data 1–3).

Generation of functional subnetworks: The functional subnetworks, representing unbalanced processes, presented in Figs. 2B, E, 4B and 5A, B were generated by combining the Gi values (represented by the size of circles) and the STRING database (used to add functional protein–protein connections). The representative subnetworks were designed using Cytoscape software.

The biological meaning of the unbalanced processes was deciphered as follows: transcripts with significant Gi values18 (“G values”; Supplementary Data 1–3) were grouped into biological categories according to GO using the David database. Representative categories are shown in Fig. 4D, E.

Single-sample gene set enrichment analysis (ssGSEA)

ssGSEA was performed on the GEO dataset using the publicly available server GenePattern to obtain MSigDB-based hallmark gene set scores. Enrichment scores were then used to compare samples with/without Act/Inf phenotypes as determined by SA as described above (Fig. 5C, Supplementary Data 4).

Statistical analysis

Significant differences between experimental conditions and experimental reproducibility were determined using the Student’s t-test (two tails, two samples equal variance) or ANOVA as indicated in figure legends; p values of ≤0.05 were considered statistically significant. All data were represented as the mean ± S.E. (standard error of the mean) if not indicated otherwise. All experiments were performed minimally in biological triplicate if not indicated otherwise.