Laser interstitial thermal therapy and adjuvant pembrolizumab in recurrent high-grade astrocytoma: a Phase 1/randomized Phase 2b trial

Domain: chemistry
Published: 2026-02-26
Editorial outcome: minor_revision

Abstract
Immune checkpoint inhibitors (ICIs) show minimal efficacy in recurrent high-grade astrocytoma (rHGA). Laser interstitial thermal therapy (LITT), a minimally invasive cytoreductive approach, may prime rHGA for ICI response. A phase 1/randomized phase 2b trial (ClinicalTrials.gov: NCT02311582) was designed to test pembrolizumab in combination with LITT in patients with rHGA. Nine patients were enrolled in the phase I dose-escalation lead-in study. No dose-limiting toxicities were observed and 200 mg of pembrolizumab every three weeks was determined as the recommended phase 2 dose. The phase 2b study was initially designed to randomize (up to 45) patients 1:1 to either LITT followed by pembrolizumab (LITT+PEM) or non-LITT surgery followed by pembrolizumab (NLS+PEM). Phase 2’s primary endpoint was progression-free survival (PFS); secondary endpoints included overall survival (OS), safety, and immune signature. After 21 patients, based on an independent Data and Safety Monitoring Committee request of unscheduled interim review of accumulating efficacy data, randomization stopped as benefit from NLS+PEM appeared limited, and the subsequent 24 patients received LITT+PEM. The pre-specified study endpoints were achieved. Among 39 per-protocol patients, LITT+PEM (n=33) improved median OS (11.8 versus 5.2 months) and 18-month survival (42% versus 0%) compared to NLS+PEM (n=6) (hazard ratio [HR] 0.17; 95% confidence interval [CI], 0.06–0.49; P=0.0002). Median PFS was longer in LITT+PEM (4.5 versus 1.6 months; HR 0.21; 95% CI, 0.08–0.56; P=0.0006). In an intent-to-treat sensitivity analysis (n=21), OS (HR 0.29; 95% CI, 0.10–0.88) and PFS (HR 0.30; 95% CI, 0.10–0.87) again favored LITT+PEM (n=13). Treatment was well tolerated. LITT activated non-classical monocytes, and pembrolizumab unleashed CD8⁺ T cell proliferation, clonal expansion, and coordinated memory T-cell responses. Overall, LITT+PEM is safe and may overcome rHGA immunosuppression to generate anti-tumor immunity.

Introduction
Recurrent HGA, including glioblastoma (GBM) and IDH-mutant World Health Organization (WHO) grade 3 and 4 astrocytoma, remains one of the greatest challenges in neuro-oncology. After maximal resection, radiotherapy, and temozolomide chemotherapy, salvage options for rHGA offer only modest benefit, and prognosis is dismal1. These realities highlight the urgent need for new therapies that can meaningfully extend survival and enhance quality of life for these patients.

ICIs, such as the anti-PD-1 antibodies pembrolizumab and nivolumab, have transformed care in many solid tumors by unleashing cytotoxic T-cell responses2. Yet, in HGA, ICIs, whether added to frontline chemoradiation or used alone or combined with the salvage therapy bevacizumab in recurrence, have failed to extend survival beyond these established treatments3–7. Resistance is thought to stem from the blood–brain barrier (BBB), low tumor neoantigen burden, and a highly immunosuppressive tumor microenvironment (TME) that limits CD8⁺ T-cell infiltration and function8. More recently, neoadjuvant ICIs given before maximal resection were shown to transiently “heat up” the TME by boosting tumor-infiltrating lymphocytes, TCR diversity, and peripheral T-cell clonal expansion, but this has not consistently translated into a survival advantage over adjuvant ICI alone in rHGA9–12. These mixed results underscore the ongoing challenge of stimulating durable anti-tumor immunity in rHGA, where both intrinsic tumor features and limited neoantigen availability that is further depleted following surgical resection remain major obstacles. Despite these disappointing clinical results, PD-L1 is detectable in a substantial subset of high-grade astrocytoma and tends to increase with higher WHO grade and PTEN loss, a common alteration in these tumors that is linked to immune evasion13–17, providing a biological rationale for continued interest in evaluating PD-1/PD-L1 blockade in this setting.

LITT is a minimally invasive, MRI-guided thermal ablation technique for brain tumors18–22. It is used to achieve effective cytoreduction with lower surgical risk and has yielded survival outcomes comparable to open resection23. Moreover, LITT has been shown to induce transient BBB disruption, particularly in the peritumoral region, thereby enhancing CNS penetration of systemic agents that would otherwise be excluded21,24,25. Accordingly, the National Comprehensive Cancer Network guidelines have added LITT as a cytoreductive option for poor surgical candidates26. Beyond direct tumor debulking, the ablative injury appears to remodel the peri-tumoral TME by continuously liberating intracranial antigens, driving peripheral tumor-specific immune activation that peaks weeks to months later21,27,28. Leveraging these complementary mechanisms—antigen release, BBB disruption, and local immune priming—combining LITT with immune checkpoint blockade may overcome the profoundly immunosuppressive milieu and sparse neoantigen landscape of rHGA.

Here, we show in a Phase 1 lead-in and randomized Phase 2b study that LITT followed by pembrolizumab is safe and associated with improved PFS and OS compared with non-LITT surgery followed by pembrolizumab in rHGA. Longitudinal paired 5’ single cell RNA-seq (scRNA-seq) and TCR sequencing of peripheral blood mononuclear cells (PBMCs) reveals systemic immune remodeling after LITT and during PD-1 blockade, including early activation of non-classical monocytes and increased CD8+ T-cell proliferation with clonal expansion. Coordinated memory-to-effector CD8+ T cell responses characterize long-term survivors, whereas stagnant and exhausted phenotypes predominate in short-term survivors.

Results
Phase 1: Between May 2016 and May 2017, nine patients (four males) with bevacizumab- and immunotherapy-naïve recurrent WHO grade 3 (two patients) or grade 4 (seven patients) astrocytoma, unequivocal progression per the Response Assessment in Neuro-Oncology [RANO] criteria29, and Karnofsky performance status (KPS) ≥ 60% were enrolled in a standard 3×3 Phase 1 study with LITT followed within 1 week by pembrolizumab at three dose levels of 100, 150, and 200 mg intravenously every 3 weeks (Fig. 1a). Adverse events were recorded using the Common Terminology Criteria for Adverse Events (CTCAE v. 4.0). No dose limiting toxicities were observed at any of the three dose levels, and thus 200 mg every 3 weeks was the recommended Phase 2 dose (RP2D).

Phase 2: From January 2017 to December 2022, 45 patients with bevacizumab- and immunotherapy-naïve recurrent WHO grade 4 astrocytoma, unequivocal progression per RANO criteria29, and KPS ≥ 60% were enrolled. The study originally specified 1:1 randomization between LITT+PEM and NLS+PEM, which proceeded for the first 21 patients (intent-to-treat [ITT] population: 8→NLS+PEM; 13→LITT+PEM). In mid-2019, emerging data demonstrating limited benefit from NLS+PEM11 prompted an independent Data and Safety Monitoring Committee (DSMC)-requested, unscheduled interim review of accumulating efficacy data, which suggested a survival signal favoring LITT+PEM. Protocol Amendment 12 subsequently halted randomization and assigned all remaining enrollments (n=24) to the LITT+PEM arm. The trial then continued as a single-arm LITT+PEM extension designed to test whether PFS exceeded a conservative historical benchmark set at 2.9 months, compared with the 2.4-month mPFS reported for NLS+PEM in the neoadjuvant pembrolizumab trial in rHGA11. Although upfront surgical debulking before LITT was permitted, in practice, no patient in either the randomized or single-arm expansion cohorts underwent resection prior to laser ablation. Patients with suspected secondary WHO grade 4 astrocytoma (prior grade 2–3 disease treated with chemoradiation and new imaging suggesting transformation to grade 4) were provisionally enrolled pending histopathological confirmation of grade 4 on biopsy or resection; unconfirmed cases were replaced.

In the ITT cohort, two patients randomized to NLS+PEM withdrew consent after randomization to receive off-trial LITT+PEM, and three patients randomized to LITT+PEM without confirmed grade 4 histology were treated off-trial with chemotherapy. Among the 24 patients enrolled after Amendment 12, one was excluded due to insurance denial for LITT and three for unconfirmed grade 4 disease (two of whom were subsequently lost to follow-up). Ultimately, 30 Phase 2 patients received LITT+PEM (10 randomized, 20 post-amendment). For the per-protocol analysis, we added three Phase 1 patients treated at the RP2D (dose level 3) who met Phase 2 eligibility, yielding 33 evaluable LITT+PEM and six evaluable NLS+PEM patients (n=39 for per-protocol cohort; see Fig. 1a). The primary endpoint was PFS; secondary endpoints were OS, objective response rate, safety, and immunophenotyping, with MRI assessments per RANO at baseline and every 9 weeks (Fig. 1b).

In both the ITT and per-protocol cohorts, the LITT+PEM and NLS+PEM arms were well matched for key baseline characteristics, including age, sex, KPS, initial diagnosis, MGMT promoter methylation, IDH mutation status, prior recurrences, tumor location, and treatment center (Table 1). All patients had at least one prior recurrence, i.e., they were at second recurrence or beyond at trial entry, and several had reached a third (or greater) recurrence, underscoring the cohort’s advanced disease stage. After the early closure of the NLS+PEM arm, we established a contemporaneous historical control (CHC) of 13 recurrent, wild-type-IDH (wt-IDH) GBM patients treated at Washington University with LITT followed by non-ICI therapies (temozolomide, lomustine, carboplatin, or bevacizumab; see Supplementary Fig. S1a and “Methods” for cohort creation). This CHC was intended to bolster our comparison and isolate the effects of LITT without ICI. The CHC cohort (all wt-IDH; 61.5% female; 53.8% MGMT-methylated) closely matched our 29 wt-IDH GBM patients in the LITT+PEM arm, although CHC tumors were disproportionately frontal and included one cerebellar and one corpus-callosal lesion (Table 1 and Supplementary Table S1).

Safety
The majority of study treatment-attributed adverse events were mild (grades 1–2), with the most common being elevated alanine aminotransferase (six patients; 13%), diarrhea (six patients; 13%), muscle weakness (five patients; 11%), pruritus (four patients; 9%), nausea (four patients; 9%), and arthralgia (four patients; 9%) (Table 2). Grade 3–4 events requiring hospitalization included diarrhea (two patients; 4%) and muscle weakness (two patients; 4%) and resolved with treatment. No grade 3–4 event occurred in more than two patients. There was no grade 5 events.

LITT followed by pembrolizumab is associated with extended survival in rHGA
Summaries of OS and PFS for all Phase 1 and Phase 2 enrollments, along with best response rates and representative brain images for the 39-patient per-protocol cohort, are shown in Fig. 2. In the randomized ITT population, LITT+PEM achieved a median OS of 21.3 months with 54% OS18, compared to 6.7 months with 12.5% OS18 in the NLS+PEM arm (HR 0.29; 95% CI, 0.1–0.88; log-rank P=0.021), and median PFS improved to 7.8 months versus 2.2 months (HR 0.30; 95% CI, 0.10–0.87; log-rank P=0.020) (Fig. 3a). Because randomization was halted early and the trial subsequently continued with single-arm enrollment, these between-arm comparisons are exploratory. In the per-protocol cohort, LITT+PEM extended median OS to 11.8 months with 42% OS18, versus 5.2 months with 0% OS18 (HR 0.17; 95% CI, 0.06–0.49; log-rank P=0.0002), while median PFS increased to 4.5 months compared with 1.6 months (HR 0.21; 95% CI, 0.08–0.56; log-rank P=0.0006) (Fig. 3b). Since the ITT population includes five randomized patients who did not complete protocol-defined treatment or lacked confirmed grade 4 histology and were excluded from the per-protocol analysis, differences in cohort composition and censoring likely contributed to the ITT median estimates, although no single extreme outlier drove the ITT survival curves. Finally, as a descriptive comparison to the conservative historical benchmark of 2.9 months specified for the amended single-arm extension, the per-protocol LITT+PEM cohort had an mPFS of 4.5 months (95% CI, 3.4–10.3), with the lower confidence bound exceeding 2.9 months.

Among wt-IDH GBM patients, LITT+PEM again outperformed NLS+PEM, with median OS of 11.1 versus 4.8 months (HR 0.16; 95% CI, 0.05–0.50; log-rank P=0.0003) and PFS of 3.9 versus 1.7 months (HR 0.29; 95% CI, 0.10–0.81; log-rank P=0.01) (Fig. 3c). Compared to the CHC treated with LITT and non-ICI agents, which was measured from the date of LITT rather than first pembrolizumab dosing as above, LITT+PEM maintained a significant survival benefit, with median OS of 11.3 versus 6.5 months (HR 0.27; 95% CI, 0.12–0.60; log-rank P=0.0006) and PFS of 4.1 versus 2.1 months (HR 0.30; 95% CI, 0.14–0.65; log-rank P=0.001) (Supplementary Fig. S1b); however, whether differences in tumor locations influenced these outcomes is unclear. To assess whether the survival advantage was confined to the early patients enrolled before the DSMC review, we compared Kaplan–Meier curves for the randomized and post-amendment expansion LITT+PEM cohorts; for both OS and PFS, these curves overlapped closely with one another and remained clearly separated from the NLS+PEM curve, indicating that the benefit was observed across early and later LITT+PEM patients (Fig. 4). Notably, MGMT promoter methylation did not affect survival in any cohort (Supplementary Fig. S2). In multivariable Cox hazard models adjusted for MGMT status, age, sex, IDH mutation, KPS, enrolling center, and number of prior recurrences, treatment arm was the only consistent independent predictor of both OS and PFS in the per-protocol population (Supplementary Table S2).

Best response analysis showed that in the LITT+PEM arm, nine patients (27%) had partial response (PR) and 14 (42%) had stable disease (SD) (Fig. 2b, c), while none achieved a complete response (CR), and ten patients (30%) had progressive disease (PD) including four with early progression (i.e., baseline MRI was the only obtained imaging). The resulting overall response rate (CR+PR) and disease control rate (CR+PR+SD) were 27% and 70%, respectively. The maximum duration of response was 32.8 months (from best response to progression) (Supplementary Table S3). In contrast, all six NLS+PEM patients experienced PD (Fig. 2).

In summary, patients treated with LITT plus pembrolizumab experienced significantly improved OS and PFS, as well as higher response and disease control rates, compared to those undergoing NLS plus pembrolizumab or those receiving LITT followed by a non-ICI therapy.

LITT is associated with early activation of non-classical monocytes
LITT-generated hyperthermia induces immunogenic cell death, endothelial injury, and BBB disruption, especially in the peritumoral region20,21,24. Because direct CNS sampling carries significant safety risks, we turned to peripheral immune profiling. Local reactivation of the GBM TME—such as via Tumor Treating Fields30,31, LITT21,27,28, or other localized immune-activating strategies32—has been shown to elicit robust responses detectable in PBMCs, which in turn serve as TME surrogates. Accordingly, serial PBMCs were collected from Phase 2 patients before NLS or LITT (Pre-Procedure), before cycle 1 (Post-Procedure), and before cycles 2, 4, 6, 9, 17, and 34 of pembrolizumab. For correlative analysis, we processed and analyzed only the University of Florida (UF) specimens to minimize batch effects and control costs and because five of the six evaluable NLS+PEM control patients were enrolled at UF. Seventeen UF patients (five NLS+PEM and 12 LITT+PEM) yielded 69 PBMC samples, generating 427,655 single cells by paired 5′ scRNA-seq and TCR α/β V(D)J sequencing. Single cells were clustered with Seurat33,34, visualized via UMAP35, resolving 32 immune cell types annotated as previously described30,31, and T cells interrogated for TCR diversity using Simpson’s index36 (Fig. 5a).

Because LITT-induced tissue injury primarily triggers innate immune activation, we first examined NK cells, dendritic cells, and classical and non-classical monocytes in PBMCs collected before and after LITT (before C1 pembrolizumab). Although overall subset frequencies were unchanged, non-classical monocytes declined slightly (Supplementary Table S4). To assess functional shifts of non-classical monocytes in response to LITT, we applied the optimal transport algorithm Earth Mover’s (Wasserstein) Distance Metric, a machine-learning approach adept at detecting subtle transitions in cell state37,38, and observed a marked post-LITT transition of non-classical monocytes into subcluster 10—a change that was absent after NLS (Fig. 5b). Non-classical monocytes, which generally acquire migratory, phagocytic, and antigen-presenting functions in response to tissue injury, including that caused by LITT, can infiltrate damage sites to promote repair and adaptive immune activation39,40. Among the 12 LITT+PEM patients showing a shift into subcluster 10, we computed log₂ fold-changes between pre- and post-LITT samples in non-classical monocytes, selected the 500 most differentially expressed genes ranked by false discovery rate (FDR), and grouped them into five functional pathways, i.e., adhesion/migration, macrophage differentiation, T cell activation, metabolic processes, and inflammatory response, to derive a composite non-classical monocyte activation score. This signature stratified patients into five with low activity and below-mOS [LITT+PEM (<mOS)] versus seven with high activity and above-mOS [LITT+PEM (>mOS)] (Fig. 5c). Expression of this program was concentrated in subcluster 10 (Fig. 5d) and enriched in longer survivors by gene set enrichment analysis (GSEA) (Fig. 5e), and in an age-adjusted Cox proportional hazard (CoxPH) model, higher composite activity remained independently correlated with improved survival (Fig. 5f). To assess which components of this program might underlie the survival association, we used the same five pathway grouping to derive pathway-specific activity scores and fit exploratory age-adjusted Cox models. Higher cell-adhesion, metabolic-process, and inflammatory-response activity were each and together significantly associated with improved survival, whereas T-cell activation showed a similar but non-significant trend, and macrophage differentiation was not correlated with outcome (Fig. 5f), suggesting that the survival signal captured by the 500-gene signature is largely accounted for by coordinated changes in adhesion, metabolic, and inflammatory programs within LITT-induced non-classical monocytes.

In short, local LITT specifically and acutely induces in the periphery a shift in non-classical monocytes to a new activated cell state, marked by high expression of pathways for migration, and metabolic and inflammatory processes—an effect not observed after NLS. This shift in cell state correlates with extended survival among LITT+PEM-treated patients and could serve as a candidate correlative biomarker of response.

LITT+PEM is associated with heightened CD8+ T cell activation and clonal expansion
Next, we investigated whether LITT-induced innate activation could influence the immune response to pembrolizumab therapy. To this end, we extracted single T cell transcriptomes and corresponding TCR α/β V(D)J sequences from PBMCs gated on CD3, then sub-clustered them using Seurat, UMAP, and Simpson’s diversity index33–36, identifying 77,947 CD3+ cells belonging to 11 immunologically recognized T cell subtypes (Fig. 6a). In analyses of pre- vs. post-pembrolizumab (C1 vs C2) samples, cycling effector CD8+ T cells (TYMS+GZMB+MKI67+) in the LITT+PEM group showed markedly increased expression of both the 556-gene T cell activation pathway (GO: 0042110) and the 761-gene adaptive immune response pathway (GO: 0002250) (Fig. 6b, c), as well as a trend toward greater relative frequency (Supplementary Table S5), compared to the NLS+PEM group. Although most other T cell subtypes (both CD4+ and CD8+) did not differ significantly between these two groups (Supplementary Table S5), CD8+ T cells overall exhibited enhanced TCR clonal expansion, a hallmark of antigen-specific T cell activation41, in LITT+PEM patients whose survival exceeded the median, distinguishing them from those with lower-than-median survival and from individuals in the NLS+PEM group (Fig. 6d). Notably, these findings align with the survival association previously observed for the LITT-activated non-classical monocytes, suggesting that early non-classical monocyte activation by LITT may foster an enhanced T-cell response upon subsequent pembrolizumab therapy.

LITT+PEM induces a coordinated transition toward CD8+ effectors in long-term survivors
To investigate how CD8+ T cells progress in response to pembrolizumab following LITT, we examined their distribution across CD8+ T-cell subtypes in pre- and post-pembrolizumab samples from the same three patient groups: NLS+PEM, LITT+PEM (<mOS), and LITT+PEM (>mOS). At baseline (Pre-Procedure and C1), CD8+ T cells displayed largely overlapping distributions in all three groups. However, upon pembrolizumab administration, while all groups generated cycling T cells, the LITT+PEM (>mOS) group quickly diverged by C2, producing effector memory T cells that persisted throughout treatment (Fig. 7a, red arrow). To follow the evolution of these CD8+ T-cell subtypes, we again applied the Wasserstein distance metric. In this analysis, arrows depict shifts between cell states, with arrow length and thickness representing the magnitude of movement and the number of cells in transit, respectively. In NLS+PEM patients, CD8+ T cells were relatively stagnant, showing minimal and uncoordinated evolution after NLS or pembrolizumab. In contrast, LITT+PEM patients exhibited more coordinated transitions into memory states, including central memory (CM), following LITT (Pre-LITT to C1), especially in those with survival exceeding the median. After initiating pembrolizumab (C1 to C2), both LITT+PEM subgroups began exiting memory states, but in longer-surviving patients, more cells progressed into the effector T cell cluster, whereas in shorter-surviving patients, cells shifted predominantly into exhausted states, including exhausted and stem-like exhausted clusters. With continued pembrolizumab exposure (C2 to C4), the shorter-surviving group showed further movement toward T cell exhaustion, whereas the longer-surviving group maintained a more robust effector phenotype throughout treatment (Fig. 7b).

To confirm that CM CD8+ T cells play a pivotal role in regulating response to LITT+PEM and survival, we tracked single TCR in all seven CD8+ T cell subtypes that originated from CM clones and followed their progression. In LITT+PEM (>mOS) patients, CM-derived cells transitioned into many other CD8+ subtypes following pembrolizumab (C1 to C2), resulting in fewer exhausted cells and more effectors. By contrast, in both NLS+PEM and LITT+PEM (<mOS) patients, CM-derived clones failed to regenerate a full range of subtypes but instead predominantly exhausted and stem-like cells (Fig. 8a). To verify that CM CD8+ T cells themselves drive these effects, we followed the transitions of CM cells directly using the Wasserstein distance metric. In LITT+PEM (>mOS) patients, CM clones displayed sustained proliferative capacity, evidenced by a significant, coordinated shift that generated progeny across the entire CD8+ T cell spectrum. In contrast, CM clones in LITT+PEM (<mOS) patients primarily transitioned into the exhausted states, and transitions in the NLS+PEM group were uncoordinated and largely random (Fig. 8b).

The results so far indicate that the CD8+ T-cell profile in longer-surviving LITT+PEM patients began to diverge from those in shorter-surviving patients primarily after the initiation of pembrolizumab. To investigate whether major genomic alterations might explain these differing outcomes, we obtained whole-exome sequencing from available baseline tumor samples, including four NLS+PEM, three LITT+PEM (<mOS), and four LITT+PEM (>mOS), to assess tumor mutational burden (TMB), a metric known to correlate with ICI responses42. No significant differences emerged in functional TMB, stop-gain or stop-loss single-nucleotide polymorphisms, or insertions/deletions across the three groups (Supplementary Fig. S3a). Similarly, there were no notable differences in the frequency of mutations, deletions, or amplifications in known canonical glioma-related genes (e.g., CDKN2A/B, RB1, EGFR, PTEN, TP53)43 or in DNA excision repair and microsatellite stability genes (e.g., MLH1, MSH2, POLE)44 (Supplementary Fig. S3b). Planned analyses of intratumoral PD-L1 expression could not be performed in a robust manner because only five tumors (three NSL+PEM and two LITT+PEM) had sufficient residual tissue; in this small subset, RNA-seq–derived PD-L1 levels did not differ appreciably between arms, and the sample size was too limited for any formal association with outcome (Supplementary Fig. S3c).

Next, to determine whether differential immune checkpoint expression might be correlated to observed outcome differences, we analyzed mRNA levels of PD-1 and five additional checkpoint markers, including CTLA-4, TIM-3, LAG-3, TIGIT, and IDO-1, across all T cells at the first three timepoints (Pre-Procedure, Post-Procedure, or C1, and C2). All markers except IDO-1 were comparably expressed and distributed among the three patient groups. By contrast, IDO-1 initially showed elevated expression shortly after LITT within memory CD8⁺ T-cell populations (central memory, memory precursor effector, and effector memory) in the LITT+PEM (>mOS) cohort (Fig. 8c), potentially reflecting the higher activation and proliferative capacity of these cells in that group. IDO-1 encodes indoleamine 2,3-dioxygenase 1, an enzyme that regulates tryptophan metabolism; while transient increases in IDO-1 can be an early indicator of T-cell activation, chronic upregulation can foster immunosuppression by limiting effector T-cell function and promoting regulatory T-cell activity45,46. Notably, upon pembrolizumab treatment, IDO-1 expression was downregulated in most T cells, except for a small subset of cycling CD8⁺ T cells, in the LITT+PEM (>mOS) group, compared to its persistent expression across multiple T-cell clusters (particularly the exhausted subset) in the LITT+PEM (<mOS) and NLS+PEM groups (Fig. 8c, red arrow). These findings suggest that the improved responses in longer survivors could be associated with LITT+PEM-driven downregulation of alternate immune checkpoints such as IDO-1.

Taken together, LITT plus pembrolizumab was safe and well tolerated and demonstrated better survival outcomes than NLS plus pembrolizumab in rHGA. Immune profiling revealed early non-classical monocyte activation, enhanced T-cell responses, and reduced IDO-1 mRNA expression in those with prolonged survival. These findings suggest that localized thermal ablation may boost ICIs by creating a more favorable immune TME and mitigating key inhibitory pathways.

Methods
Study design and treatments
Approval. Human subject work was performed in accordance with an approved protocol by the respective internal review boards at the Washington University in St. Louis School of Medicine, the University of Florida College of Medicine, and the University of Southern California Keck School of Medicine, in accordance with the Declaration of Helsinki. A written informed consent was obtained from each human participant before study procedures and analysis were performed. All enrollments were completed at Washington University and the University of Florida. The correlative analysis was performed at the University of Southern California. The study’s ClinicalTrials.gov identifier number is NCT02311582.

Study design. The study began with a Phase 1 lead-in to assess pembrolizumab safety in rHGA, followed by a randomized Phase 2 at a ratio 1:1 comparing NLS+PEM to LITT+PEM in rHGA. Randomization was computer-generated using a blocked scheme generated in SAS (v9.3+) and conducted centrally via the electronic data capture system REDCap. After emerging data demonstrated the limited benefit of ICI monotherapy in recurrent GBM (ClinicalTrials.gov, NCT02337491)5, the independent Data and Safety Monitoring Committee (DSMC) requested an unscheduled interim review of accumulating efficacy data. On the basis of this review, the DSMC recommended closing the NLS+PEM arm and redirecting all subsequent enrollment to LITT+PEM (Protocol Amendment 12). Patients assigned to LITT+PEM could undergo optional surgical debulking prior to LITT; if resection was performed, the laser ablation then targeted only the residual tumor. To bolster our control framework and isolate LITT’s effects without ICI, we then introduced a contemporaneous historical control (CHC) cohort treated with LITT plus non-ICI therapies (LITT+Non-ICI; Table 1). The primary interventions were LITT or NLS followed by adjuvant pembrolizumab, initiated within 1-week post-procedure (Fig. 1b). Patients receiving at least one dose of pembrolizumab were evaluable for survival and response analyses.

Patient population. Between May 2016 and May 2017, nine patients with radiographically recurrent WHO grade 3 or 4 astrocytoma were enrolled in the Phase 1 dose-escalation (3×3) cohort, receiving pembrolizumab at 100 mg, 150 mg, or 200 mg IV every 3 weeks following LITT (three patients per dose level). In Phase 2, from January 2017 to December 2022, 45 patients with histologically confirmed recurrent WHO grade 4 astrocytoma enrolled: 21 in the randomized cohort and 24 assigned to LITT+PEM after Amendment 12 (Fig. 1a). Key inclusion criteria were unequivocal tumor progression by RANO, a minimum of 12 weeks from frontline chemoradiation, a candidate for LITT or surgical resection, and Karnofsky performance status (KPS) ≥ 60%. Patients with suspected secondary WHO grade 4 astrocytoma—those with prior grade 2 or 3 tumors treated with chemoradiation, now with radiographic changes suggestive of transformation to grade 4—were conditionally enrolled and received pembrolizumab only after histopathological confirmation of grade 4 post-LITT or post-NLS. Key exclusion criteria included dexamethasone > 4 mg/day at registration, prior immunotherapy (anti–PD-1/PD-L1 or anti-CD137), and bilateral multifocal or leptomeningeal disease. IDH mutation and multiple prior recurrences were permitted in both phases (Table 1).

Contemporaneous historical controls. We identified contemporaneous control patients treated at Washington University School of Medicine between 2016 and 2022 for histologically confirmed rHGA via retrospective chart review (see Supplementary Fig. S1a for the schema and criteria). Inclusion criteria mirrored those of the prospective study and required biopsy-confirmed recurrence. We excluded patients with newly diagnosed grade 4 astrocytoma, those enrolled in other LITT trials, and individuals receiving immunotherapy. Of 122 rHGA patients treated with LITT during this period, 13 met all eligibility criteria.

Laser interstitial thermal therapy. LITT is an FDA-cleared, minimally invasive technique for thermal ablation of brain lesions. In our study, we employed the Monteris NeuroBlate Laser Ablation System (Minnetonka, MN) under MRI guidance. Following a small scalp incision and cranial burr hole, an MRI-compatible laser probe is advanced into the tumor. The system then delivers thermal energy, raising the lesion core temperature to ≈70°C to induce coagulative necrosis, while peritumoral temperatures gradually climb to a peak of 40–45°C, thereby limiting collateral injury18–21.

Statistical analysis plan for clinical and survival comparisons. Phase I: A traditional 3+3 design was used to determine the recommended Phase 2 dose of pembrolizumab when combined with LITT for the treatment of rHGA.
Phase 2: The primary endpoint was progression-free survival (PFS); secondary endpoints were overall survival (OS), safety, objective response rate, and immune-response signatures. In both the randomized intent-to-treat and per-protocol cohorts, PFS and OS were measured from first pembrolizumab dose to radiographic progression (per RANO) or death. For the contemporaneous historical control (LITT+non-ICI), both endpoints were calculated from the date of LITT. Patients lost to follow-up or alive at data cutoff were censored at their last known alive date or at cutoff, whichever came first.
Initial Phase 2b calculations (one-sided α=0.15) assumed six-month PFS of 40% in the NLS+PEM arm—based on published rates (<10–30%) for common salvage regimens such as bevacizumab, temozolomide, carboplatin, carmustine, lomustine, etoposide, irinotecan50–54—and 65% in LITT+PEM, requiring 17 patients per arm for 80% power; we thus targeted 20 per arm. Randomization was not stratified by IDH mutation, MGMT promoter methylation, or other baseline factors because this planned small sample size did not permit reliable balance across multiple strata. After 21 randomized patients, emerging data showing limited ICI monotherapy benefit prompted the independent DSMC to request an unscheduled interim review of the accumulating trial data. This interim analysis was not prespecified in the original protocol and did not use a formal alpha-spending or stopping boundary. Descriptive comparisons at that time indicated longer PFS and OS in the LITT+PEM arm, and, together with the external evidence, led the DSMC to recommend closing the NLS+PEM arm and continuing the trial as a single-arm LITT+PEM extension in the spirit of minimizing harm, potential or perceived, to subjects and to complete the study in a timely manner. Under Amendment 12, the control arm closed, and the trial continued as a single LITT+PEM cohort. For this extension, we conservatively used a 2.9-month median PFS for the NLS+PEM arm (the upper 95% confidence bound rather than the 2.4-month median reported)11. Under that assumption, enrolling 27 additional patients (to achieve 20 evaluable GBM cases) provided 80.9% power (two-sided α=0.05, one-sample log-rank) to detect an increase in median PFS to 5.1 months versus the 2.9-month benchmark.
In the final analysis, we descriptively compared the PFS distribution of the per-protocol LITT+PEM group, estimated by Kaplan–Meier with its 95% confidence interval, to this 2.9-month benchmark. Two-sample log-rank tests comparing LITT+PEM with NLS+PEM in the randomized ITT and per-protocol populations, and with the contemporaneous LITT+non-ICI cohort, were performed as exploratory analyses to contextualize the magnitude of effect and do not preserve the nominal type I error. Kaplan–Meier curves for OS and PFS were generated with log-rank tests for these unadjusted comparisons. Univariate and multivariable Cox proportional-hazards models adjusted for age, sex, IDH mutation, MGMT methylation, enrolling center, baseline KPS, and number of prior recurrences provided hazard ratios (HR) with 95% confidence intervals and p-values. All tests were two-sided at α=0.05. Analyses were performed in SAS v9.4 (Cary, NC).

Tumor tissue processing and analysis
All available tumor samples from the study participants that remained after clinical diagnostics were prioritized for whole-exome sequencing and, when sufficient material was left, bulk RNA-seq. Exploratory analyses of intratumoral PD-L1 expression by immunohistochemistry were prespecified but could not be completed in a sufficient number of cases because most biopsies were very small and were exhausted by diagnostic testing and other genomic assays. Only five tumors (three NLS+PEM and two LITT+PEM) had adequate material for bulk RNA-seq–based estimation of PD-L1 expression.

Single cell RNA-seq (scRNA-seq) analysis of PBMCs
Portions of the following scRNA-seq methods are adapted from our prior work30, published under a Creative Commons license, CC BY 4.0, and are included here in full for completeness:

Sample processing. Cryopreserved PBMCs obtained from the study participants were rinsed in PBS, and cell viability was assessed using Trypan Blue staining. Single-cell suspensions were then prepared and applied to the Chromium Single Cell Chip (10x Genomics) as per the instruction provided by the manufacturer. Subsequently, single-cell RNA-seq libraries were generated using the Chromium Next GEM Single Cell 5′ v2 (Dual Index). To ensure consistency, all patient samples and the corresponding libraries were processed simultaneously in a single batch. Sequencing of the single-cell libraries was performed on an Illumina NovaSeq system, utilizing an 8-base i7 sample index read, a 28-base read 1 for capturing cell barcodes and unique molecular identifiers (UMIs), and a 150-base read 2 for the mRNA insert.
Data processing. The main operations were performed using the Seurat R package (3.2.2)33,34, unless otherwise stated. When option parameters for the function deviated from the default values, we provided details of the changes accordingly. Most of the changes to the default options were made to accommodate and leverage the large size of the dataset. Cell Ranger Aggregation: The raw sequencing data were processed using cellranger mkfastq and cellranger multi commands of the Cell Ranger package as described in the TCR clonotyping section. Results from all libraries and batches were pooled together using the command cellranger aggr without normalization for dead cells, as it will be handled downstream. The filtered background feature barcode matrix obtained from this step was used as input for sequential analysis. Normalization of UMI: Using the global scaling normalization method, the feature expression for each cell was divided by the total expression, multiplied by the scale factor (10,000), and log transformed using the Seurat R function NormalizeData with method “LogNormalize”. Seurat aggregation and correction for batch effects: As the counts were from two different batches, to align cells and eliminate batch effects for dimension reduction and clustering, we adopted the multi-dataset integration strategy as previously described34. Briefly, “anchors cells” were identified between pairs of datasets and used to normalize multiple datasets from different batches. We chose a reference-based, reciprocal PCA variant of the method detailed in the Seurat R package33,34. First, the previously integrated dataset was split by batches using the Seurat function SplitObject. Next, for each split object, variable features were selected using the function FindVariableFeatures. Features for integration were selected using the function SelectIntegrationFeatures, and PCA performed for each split object on the selected features. The anchor cells were identified by using the function FindIntegrationAnchors with the reference chosen as the largest among two batches and the reduction option set to “rpca”. Finally, the whole datasets from the two batches were reintegrated using the function IntegrateData with the identified anchor cells.

Uniform manifold approximation and projection (UMAP) dimension reduction. The integrated multiple batch dataset was used as input for UMAP dimension reduction35. The feature expression was scaled using the Seurat function ScaleData, followed by a PCA run using the function RunPCA (Seurat) with the total number of principal components (PC) to compute and store option of 100. The UMAP coordinates for single cells were obtained using the RunUMAP function (Seurat) with the top 75 PCs as input features (dims=1:75) with min.dist=0.75 and the number of training epochs n.epochs=2000. Clustering of cells: We relied on a graph-based clustering approach implemented in the Seurat package, which embeds cells in a K-nearest neighbor graph with edges drawn between similar cells and partitions nodes in the network into communities. Briefly, a Shared Nearest Neighbor graph was constructed using the FindNeigbhors function with an option dimension of reduction input dims=1:75, error bound nn.eps=0.5. This function calculates the neighborhood overlap (Jaccard index) between every cell and its k.param nearest neighbors55. The graph was partitioned into clusters using the FindClusters function with different values for the resolution parameter. The differential expressed gene markers for each cluster were found using the FindAllMarkers function with the option of only returning positive markers and a minimal fraction of cells with the marker of 0.25. The default Wilcoxon Rank Sum test was used to calculate statistical differences in each cell cluster. Alternative 2D embeddings, including t-SNE, were evaluated using quantitative metrics of within-cluster compactness, between-cluster separation, k-nearest-neighbor purity, silhouette width, and kNN-graph modularity, and UMAP most faithfully preserved the predefined clusters and trajectories; UMAP embeddings were therefore used for visualization and as the feature space for downstream distance-based analyses.

Cell type annotation. To delineate specific cell types within the data, cell type labels were assigned manually to clusters emerging from UMAP analysis. This annotation process was guided by the expression profiles of a set of marker genes, which are characteristic of various cell types including T cells, B cells, natural killer (NK) cells, monocytes, dendritic cells (DCs), myeloid-derived suppressor cells (MDSCs), megakaryocytes, red blood cells (RBCs), CD34+ stem cells, granulocytes, lymphocytes, macrophages, basophils, eosinophils, and neutrophils. The marker genes utilized for this purpose encompassed a wide range of immune response and cell differentiation indicators such as CD3D, CD3E, ID3, IL7R, CCR7, ITGB1, CD95, TCF7, CD3D, CD3E, CD4, S100A4, CCR10, FOXP3, IL2RA, TNFRSF18, IKZF2, CTLA4, IL2, IL4, IL13, IL17A, CD3D, CD3E, CD8A, CD8B, CCL4, GZMA, GZMB, GZMH, GZMM, GZMK, IFNG, GNLY, TNF, PDCD1, LAG3, CD79A, CD79B, CD19, JCHAIN, GNLY, NKG7, CD16, NCAM1, KIR2DL4, SIGLEC7, CD14, LYZ, S100A8, S100A9, LGALS3, FCN1, FCGR3A, MS4A7, FCER1A, CST3, ITGAM, ITGAX, CLEC10A, CLEC9A, THBD, CD1C, LILRA4, CLEC4C, TLR7, TLR9, ITGAM, CD33, CD3D, CD3E, CD14, CD19, FUT4, CEACAM1, HLA-DRA, HLA-DRB1, HLA-DRB5, PPBP, PF4, ITGA2B, ITGB3, PEAR1, CD42D, CD59, HBG1, HBG2, HBB, CD34, CCR3, CD11b, CD13, CD18, CD229, CRACC, CD14, CD68, CD36, CD164, LAMP1, CD44, CD69, EMR1, MPO, CD62L, CD3D, CD3E, CD4, CD8A, CD8B, NKG7, GNLY, CD14, LYZ, FCER1A, CLEC10A, LILRA4, CLEC4C, CD79A, CD79B, HBB, PPBP, and PF4. T cells were further divided into clusters to annotate subpopulations: naive CD4, central memory CD4, central memory CD8, anergic CD4, activated CD4, Treg, exhausted CD4, stem-like CD8, NKT, exhausted CD8, effector CD8, naive CD8, cytotoxic CD4, and effector memory CD8 using the following marker genes: CD3D, CD4, CD8A, CTLA4, PDCD1, TIGIT, FOXP3, CCR7, GZMK, GZMB, GZMH, IL7R, CCL5, KLRB1, TRAV16, TRAV17, CX3CR1, CCL4, TRDC, CD69, FOS, BATF, IL2RB, TBX21, EOMES, PRDM1, and ICOS. Monocyte cells were further divided into classical and non-classical monocytes using the following marker genes: CD14, IRF8, CCR2, CX3CR1, and FCGR3A.

Differential expression analysis. To evaluate transcriptomic changes in non-classical monocytes between the Pre-LITT and C1 timepoints from the LITT+PEM cohort, a differential expression (DE) analysis was performed using an aggregated “pseudo-bulk” RNA-seq approach, followed by heatmap visualization of the top differentially expressed genes:
• Pseudo-bulk aggregation: the Seurat object was subset to retain non-classical monocyte cells at the Pre-LITT and C1 timepoints. For each patient and timepoint, the raw count data were summed across all cells in that group, yielding a pseudo-bulk expression profile (genes × samples). By using this aggregation strategy, cell-level biases are circumvented, and an approximation of bulk RNA-seq at the patient–timepoint level is achieved.
• Differential expression with edgeR:
a. A DGEList object was created from the aggregated counts.
b. Normalization factors were computed via calcNormFactors.
c. A design matrix was constructed (~Patient+TimePoint) to compare Pre-LITT vs. C1.
d. After dispersion was estimated, a negative binomial model was fit, and likelihood ratio tests (glmLRT) were performed.
e. The resulting DE table was filtered for significance, capturing log2 fold change (log2FC) and FDR values.
• Top-gene selection and visualization: The top 500 differentially expressed genes by FDR were selected for further inspection. To highlight variation among patients, the log2FC (Pre-LITT vs. C1) were transformed into z-scores, and these scaled values were visualized in a heatmap. The heatmap columns were ordered according to patient survival times, thus facilitating a direct comparison of expression patterns across different clinical outcomes.

Z-score calculation and GSEA ranking. In the non-classical monocytes of the Experimental group, gene expression differences between time point Pre-LITT and time point C1 were assessed using log2FC values. For each gene, the mean log2FC was calculated separately for the short-term group (xS)—LITT+PEM (<mOS)—and for the long-term group (xL)—LITT+PEM (>mOS). The sample standard deviations (SD) were estimated for each group, and the number of patients in each group was denoted as nS (short-term or <mOS) and nL (long-term or >mOS). A z-score was then computed to represent the difference between the two groups using the following equation:
z = (xL - xS) / sqrt(SD2L/nL + SD2S/nS), (1)
where xL - xS indicates the difference in mean log2FC [(>mOS) minus (<mOS)]L, and S the denominator represents the standard error. In this formulation, positive z-scores are indicative of higher expression in LITT+PEM (>mOS) survivors, whereas negative z-scores are indicative of higher expression in LITT+PEM (<mOS) survivors. In addition, genes with missing values or undefined z-scores were removed, and the remaining z-scores were compiled into a two-column rank file (gene name and z-score). The file was then sorted in descending order by the z-score and used as input for a GSEA.

GSEA for pathway identification and confirmation. Based on the z-scores computed as described above, a set of the top 500 differentially expressed genes by the lowest FDR was identified. This subset of genes, showing the most statistically significant differences between the Pre-LITT and C1 timepoints in non-classical monocytes in the LITT+PEM cohort, was treated as a custom gene set of interest. To further elucidate the biological pathways that may be driving survival differences, GSEA was performed using two inputs: (i) the ranked list of all genes (ordered by descending z-score), and (ii) the newly defined custom gene set of 500 genes. Following the standard procedure, GSEA was carried out by using the ranked list to calculate an enrichment score (ES), which reflects the degree to which the custom gene set is overrepresented at either the top (positive z-scores) or the bottom (negative z-scores) of the ranked list. Statistical significance was assessed by permuting the gene labels and recalculating the ES to generate a null distribution. This permitted estimation of the nominal P-value and the FDR q-value for enrichment.

Cox proportional hazard model for survival prediction using the identified GSEA gene set. The scaled and normalized log2FC of the top 500 genes by FDR was tested as a predictor of overall survival, beyond known clinical factors such as patient age, in a Cox proportional hazards model. First, the scaled and normalized log2FC of these 500 genes for each patient was summarized by taking the mean across the genes for each individual patient. Then, a binary covariate was created by categorizing patients into “high” or “low” expression groups based on the median of this summarized value for the cohort. This binary pathway-expression variable was combined with age (a continuous variable) as a covariate to determine whether patients with pathway-expression values above the median had a statistically different hazard ratio (i.e., risk of death or relapse) compared with those whose expression was below the median, after adjusting for age. The resulting hazard ratio and its confidence interval provided an estimate of how strongly pathway expression correlates with survival outcomes. The Schoenfeld residuals test56 was used to assess the proportional hazards assumption in the Cox regression model, and no covariates demonstrated significant violation, indicating that the assumption held for the model.

Pathway activity signal change comparison. To compare pathway activity changes between NLS+PEM and LITT+PEM arms across consecutive timepoints, T-cell subpopulations were first defined based on clustering annotations derived from the Seurat workflow. For each pathway of interest (specified in a GMT file), the relevant genes were extracted, and the average expression of those genes computed within each patient and time point. Specifically, the following steps were performed:
• Pathway Gene Extraction: A GMT file containing predefined pathways was parsed to identify the set of genes corresponding to each pathway.
• Subset and Normalization: For each T cell subpopulation, a subset of the Seurat object was created. Expression values were normalized using the relative counts (RC) method (normalization.method=“RC” with a scale factor of 10^6).
• Mean Pathway Expression: For each patient and timepoint, the mean expression of the pathway genes in the corresponding subpopulation was computed.
• Ratio-Based Change Computation: To quantify the relative fold change in pathway activity between two consecutive timepoints (A and B), the ratio of mean pathway expression at timepoint B to the mean pathway expression at timepoint A was used:
Signal Change = MeanExpr B / MeanExpr A (2)
• Statistical Assessment: For each comparison, the ratio values for the NLS+PEM and LITT+PEM arms were pooled, and a Wilcoxon rank-sum test was performed to evaluate differences between arms.
All comparisons were repeated across the T cell clusters of interest and for the two pathways, T Cell Activation (GO: 0042110) and Adaptive Immune Response (GO: 0002250), defined in the GMT file. Figures were generated to visualize both the per-patient longitudinal changes (line plots) and the distribution of ratio values by treatment arm (box plots).

TCR clonotyping
Portions of the following TCR clonotyping methods are adapted from our prior work30, published under a Creative Commons license, CC BY 4.0, and are included here in full for completeness:

Sample processing. Single-cell RNA-seq libraries were generated using the Chromium Next GEM Single Cell 5′ v2 (Dual Index) alongside the Human V(D)J Amplification Kit (10x Genomics), following the manufacturer’s protocols. To ensure consistency, all patient samples and the corresponding libraries were processed simultaneously in a single batch. Sequencing of the single-cell libraries was performed on an Illumina NovaSeq system, utilizing an 8-base i7 sample index read, a 28-base read 1 for capturing cell barcodes and unique molecular identifiers (UMIs), and a 150-base read 2 for the mRNA insert.
Data processing. The 5′ single cell TCR α/β V(D)J library data were first processed using the 10x Genomics Cell Ranger package (v. 7.0.0, with Java/9.0.1, bcl2fastq/2.20.0.422 dependencies). Command cellranger mkfastq was used to convert the raw sequencing data from the bcl to fastq format, and the cellranger multi command to align to the reference genomes GRCh38 (GENCODE v. 24) and single cell clonal identification. Clonal tracking plots were created using the Immunoarch R package v0.9.0 (https://cloud.r-project.org/web/packages/immunarch/index.html) with the function trackClonotypes.

TCR clonotyping, diversity, and evolution. A clonal tracking grid was established to map the presence and characteristics of TCRβ clones at each patient timepoint, focusing on aspects such as the number of cells, changes in cell type, and pathway activity.
• Clonal Diversity by Simpson Index: All unique TCRβ clones in the CD8 T cell subset were first identified for each patient. The Simpson diversity index, D, was then defined as:
D = 1 - sum(p_i^2), (3)
where p_i is the proportion of cells belonging to clone i, and n is the total number of unique clones within a patient’s T cell repertoire. Lower values of D indicate reduced diversity (i.e., fewer clones dominate the repertoire), reflecting a higher degree of clonal expansion. Conversely, higher values of D signify a more evenly distributed clonal landscape. Thus, this index measures the probability that two randomly selected individuals from the population will be of the same clone, with values ranging from zero (no diversity) to one (maximum diversity). This ratio provides a comparative measure of clonal diversity between two specified timepoints. A ratio greater than one indicates an increase in diversity, suggesting a diversification of the clonal population over time. Conversely, a ratio less than one reflects a decrease in diversity, pointing to a homogenization of the population. This analysis facilitates the understanding of clonal dynamics over the course of the study.
• Statistical Comparison Between Cohorts: After the Simpson index for each patient was computed, 2-sample t-tests were performed to compare the mean diversity between:
NLS+PEM vs. LITT+PEM (<mOS).
NLS+PEM vs. LITT+PEM (>mOS).
The corresponding P-values were calculated under the null hypothesis that the mean diversity of the two cohorts does not differ significantly. A low P-value (P < 0.05) indicates that one cohort exhibits a statistically different clonal diversity—and thus a different degree of clonal expansion—relative to the other.

Tracking T cells Clones Across Timepoints. An Optimal Transport (OT) algorithm37,38 was used to trace CD8 T-cell clones between consecutive timepoints across three cohorts [NLS+PEM, LITT+PEM (<mOS), and LITT+PEM (>mOS)]. The analysis quantified how individual T-cell subpopulations (e.g., effector CD8+, exhausted CD8+ T cells) transitioned at both single-cell and cluster levels from one timepoint to the next. The steps were as follows:
• Data Preparation and Subset Selection: For each cohort and T cell subpopulation, single-cell expression and UMAP coordinates were collected at each available timepoint (e.g., Pre-Procedure, C1, C2, C4, etc.). Clusters of particular interest (COI) were merged where necessary to reduce sparsity.
• Cost Matrix Construction for OT: For each pair of consecutive timepoints, a cost matrix was computed based on Euclidean distances in principal component (PC) space. Let x_i and y_j be the PC embeddings for cell i in the source timepoint and cell j in the target timepoint, respectively, with the cost defined as:
C_ij = ||x_i - y_j||^2 (4)
Uniform distributions a and b were imposed over the source and target cells, respectively, ensuring equal weights for all cells.
• Computing the OT Plan: The OT plan π* was solved via the Earth Mover’s Wasserstein Distance (ot.emd) function:
min_π sum_i,j π_ij C_ij subject to sum_j π_ij = a_i, sum_i π_ij = b_j, π_ij >= 0 (5)
The resulting π* yields a principled “cell-to-cell matching” under global cost minimization.
• UMAP-based Visualization: Using the final OT plan, each source cell was mapped to a corresponding target cell, and arrows were plotted in UMAP coordinates:
Single-Cell Arrows: Each source cell had an arrow drawn to its corresponding target cell. Down-sampling was carried out to ensure consistent cell counts for plotting.
Cluster-Level Arrows: By aggregating source and target coordinates within each cluster, representative arrows were drawn for each cluster of interest.
Aggregated Arrow: Summed cluster displacements produced a single arrow reflecting the global population shift between two timepoints.
• Stacked-Bar Distributions: For each COI, the fractions of its cells in the source timepoint that mapped to each possible target cluster were computed. These fractions were represented as stacked bar charts, illustrating cluster transition patterns over successive timepoints.
• Cohort Comparisons: This pipeline was repeated for all T cell subpopulations across the three cohorts [NLS+PEM, LITT+PEM (<mOS), and LITT+PEM (>mOS)], enabling direct comparisons of clonal dynamics. Variations in both the directional shifts (arrows) and target cluster memberships (stacked bars) offered insight into how clonal populations evolved under different experimental conditions.

Target distribution of T cell clones across timepoints. Following the determination of how each source T cell mapped to a target cluster, cohort-specific target distributions were compared to assess whether different cohorts yield distinct post-mapping cluster compositions. Specifically, for each cohort (e.g., NLS+PEM vs. LITT+PEM), the frequency with which source T cells, which were drawn from one timepoint and mapped to various clusters at the subsequent timepoint, was counted. These counts were then organized into a contingency table, where rows represented target clusters and columns represented cohorts.
Once the contingency table was constructed, a Chi-square test of independence was performed to determine whether the observed target cluster distributions differed significantly across cohorts. A low P-value (e.g., P < 0.05) indicated that the distribution of mapped target clusters was not independent of cohort membership, i.e., there was evidence to suggest that different cohorts lead to different cluster-level outcomes for T cells across timepoints. Conversely, a higher P-value implied that the observed distribution of target clusters could be attributed to chance alone, indicating no statistically significant difference between cohorts in their target distributions.

Immune checkpoint expression comparison
To visualize immune checkpoint expression at multiple timepoints across three cohorts [NLS+PEM, LITT+PEM (<mOS), and LITT+PEM (>mOS)], a set of genes of interest (e.g., PDCD1, CD274, IDO-1) was defined, and FeaturePlot- and density-based visualizations were then generated as follows:
• Data Preparation and Cohort Assignment: Patients were first classified into three cohorts [NLS+PEM, LITT+PEM (<mOS), and LITT+PEM (>mOS)] based on OS. A Seurat object containing T-cell data was filtered to include only those patients overlapping with the survival dataset. The metadata was updated with a SurvivalGroup label for each cell, enabling straightforward subgroup analyses.
• Timepoint Filtering and Down-sampling: Analyses were restricted to timepoints Pre-Procedure, C1, and C2, discarding cells from other timepoints. To avoid uneven sampling, each cohort was down-sampled at each timepoint to ensure approximately equal cell counts, thereby reducing biases introduced by disproportional representation.
• Split Panel Feature Plots: For each immune checkpoint gene, FeaturePlot was used to visualize expression on UMAP embeddings. The data were split by both cohort and timepoint (i.e., Group_TimePoint), resulting in subplots such as NLS+PEM_Pre-NLS, NLS+PEM_C1, and so forth. Each subplot included cells for the relevant subgroup and colored them according to the expression level of the target gene.
• 2D Density Plots for High Expressors: In parallel, 2D density plots were generated to highlight cells with high (e.g., top 10% or 50%) expression of each gene. The set of top-expressing cells were colored by cohort, allowing comparisons of the spatial distribution of strongly expressing cells across cohorts within the UMAP space. Specifically, geom_density_2d was applied in ggplot2, with expression values used as weights, thus reflecting how each subgroup’s high-expressing cells cluster or disperse.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This study provides the first prospective evidence that LITT combined with PD-1 inhibition can extend survival and elicit systemic immune activation in recurrent high-grade astrocytoma. Despite the small control cohort and partially amended design, the consistency across clinical and immunologic endpoints strongly supports a biologic synergy between local thermal ablation and immune checkpoint blockade.\nInterpretation of IDO-1 downregulation as mechanistic driver is intriguing but speculative.\n-Frame this as an association rather than a causal link\n-Discuss whether tryptophan–kynurenine metabolites were measured to support this.\nStrong correlative claims (e.g., monocyte activation predicting survival) need caution given small n.\nExplicitly state the DSMC’s criteria for halting randomization (safety vs efficacy signal).\nConsider adding sensitivity analyses (e.g., adjusted Cox models including center, KPS, recurrence number).\nClarify whether randomization was stratified by IDH or MGMT status.", "reviewer_2": "The authors of the study present a Phase 1/2 multiinstituitional study of LITT and immunotherapy for recurrent high grade gliomas. The study demonstrates that LITT and Immunotherapy improves overall survival compared to the surgical resection and pembro group. The study demonstrates and OS of 11 months in the protocol group with a OS of 21 months in the ITT group. The results provided are encouraging and signal some change in the plasma biomarkers in T cells in response to LITT+ pembro group. There are a few significant limitations that should be addressed listed below that would enhance the paper, below:\n1. The control surgery+ Pembro group outcomes are significantly worse than known contemporaneous control groups like Checkmate 143 (OS of 10 months). How is this explained? Would be good to compare their results to known trials in this setting to see if this is applicable.\n2. The OS of ITT group is very promising but the OS is nearly 10 months different from the per protocol group. An explanation of this difference would be helpful. If there is a big outlier in the ITT group, it would be preferable to just put the protocol OS in the abstract so it is clear what the benefit of pembro and LITT is, without confounders\n3. MRIs of the partial and stable response and treatment failures would be helpful. No imaging data or case data is presented in the main figures\n4. A discussion on the use and picking of the Immune checkpoint inhibitor would be helpful as well. Is there a use for newer agents or PD-1, CTLA4, LAG3 inhibitors?", "reviewer_3": "While this manuscript described exciting new efficacious treatment of LITT + Pem in recurrent high-grade astrocytoma. However, I do have some concerns on the validity of the conclusion. The primary concern is that the study design and the statistical analysis are not aligned.\n(1) Unplanned interim analysis: The original protocol did not include an interim analysis. Conducting an unplanned interim analysis can inflate the type I error rate because no pre-specified boundaries were defined. Without a pre-specified alpha-spending plan, it is unclear what p-value threshold was used to justify early stopping.\n(2) Potential instability of early findings: Early trends in treatment effect are often unstable. The apparent superiority observed at an unplanned interim analysis may not have persisted had the trial continued, particularly when a liberal or post hoc boundary was used for decision-making.\n(3) Mismatch between amended design and final analysis: After the interim analysis, the trial was amended to a single-arm design, referencing a historical median PFS of 2.9 months. The amended design and sample size calculation were based on a one-sample log-rank test. However, in the manuscript’s analysis, two-sample log-rank tests were applied to both the per-protocol and randomized populations. This represents a clear inconsistency between the design and the analysis. If the analysis method differs from the design basis, it becomes uncertain whether the nominal power and type I error rate are still properly controlled.\n(4) Was any analysis performed on intratumoral PD-L1 expression—either comparing expression levels between the two groups or evaluating the association between PD-L1 expression and PFS/OS? This analysis was described in the statistical analysis plan but was not reported or discussed in the manuscript.\n(5) When you consider using UMAP to perform feature extraction, have you considered using some other common methods such as t-SNE? In addition, there might be some strong or weak associations, and as you have used the “FindNeighbors” function, I wish to know if you have found some structural clustering from the original sequencing. For example, some genes are highly correlated so that they potentially “work together”. This is additional advice for further analysis. You may provide some insight if it takes too much time to re-run the analysis.\n(6) For the Cox proportional hazard model that predicts the GSEA gene set, is there a mediation pathway for these genes to work? These gene sets may affect a specific biomarker, and further these biomarkers affect the survival probability. If so, would you be able to perform a mediation analysis?\n(7) I am a bit confused with Figure 6 b, seems like only “pre-procedure to C1”, and “C2 to C4” is shifting to memory precursor, and “C1 to C2” have some shift to exhausted and effector memory, would you be able to comment on these three.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "We thank the reviewers for their careful evaluation and constructive comments, which have substantially improved our manuscript. In response, we have tempered causal interpretations (e.g., IDO-1 and monocyte activation), clarified key aspects of the trial design, interim analysis, and statistical methods (including new sensitivity and descriptive one-sample analyses), and explained the rationale for our choice of pembrolizumab and available biomarker data. We have also added representative MRI images, expanded and clarified our single-cell and computational analyses, and more fully acknowledged limitations and future directions.\nReviewers’ comments are shown in black, with our responses in blue. Please note that all line-number references pertain to the Track Changes version of the revised manuscript.\nReviewer #1 (Remarks to the Author):\n1. Interpretation of IDO-1 downregulation as mechanistic driver is intriguing but speculative.\n-Frame this as an association rather than a causal link.\n-Discuss whether tryptophan–kynurenine metabolites were measured to support this.\nOur response: Thank you for this helpful comment. In the revised manuscript, we now describe IDO-1 downregulation as an association rather than implying a causal role. We did not measure tryptophan–kynurenine metabolites, as this was an unanticipated finding and we lacked sufficient remaining T cells or PBMCs from many subjects for metabolite analyses. We have added to the Discussion that future studies should directly evaluate this metabolic axis in T cells following LITT+PEM treatment to test this potential mechanism (Lines 561–566).\n2. Strong correlative claims (e.g., monocyte activation predicting survival) need caution given small n.\nOur response: Thank you for this comment. We have revised the text to present monocyte activation as a correlative finding rather than a strong predictor of survival (Lines 366–369).\n3. Explicitly state the DSMC’s criteria for halting randomization (safety vs efficacy signal).\nConsider adding sensitivity analyses (e.g., adjusted Cox models including center, KPS, recurrence number).\nOur response: We have clarified in the main text and Methods that the DSMC’s decision to halt randomization was driven by an efficacy signal, not by safety concerns. Specifically, in the context of emerging external data suggesting limited benefit of adjuvant immune checkpoint inhibitors alone, the DSMC requested an interim efficacy analysis per the monitoring plan. This analysis demonstrated a statistically significant survival advantage for the LITT+PEM arm compared with the NLS+PEM arm, leading the DSMC to recommend stopping further randomization (Lines 172–179; 958–963; 1031–1040). We have also added adjusted Cox sensitivity analyses (adding center, KPS, and recurrence number), which yielded results consistent with the primary analysis (Updated Supplemental Table S2) (Lines 294–297).\n4. Clarify whether randomization was stratified by IDH or MGMT status.\nOur response: Randomization was not stratified by IDH or MGMT status. The initial Phase 2b sample size calculations (one-sided α = 0.15) assumed a 6-month PFS of 40% in the NLS+PEM arm, which was based on published rates (<10–30%) for common salvage regimens such as bevacizumab, temozolomide, carboplatin, carmustine, lomustine, etoposide, and irinotecan [1-5] and 65% in the LITT+PEM arm, requiring 17 patients per arm for 80% power. Given this relatively small, planned sample size, randomization was conducted without stratification factors, including IDH or MGMT status. We have added this clarification in the Methods sections (Lines 1031-1033).\nReviewer #2 (Remarks to the Author):\n1. The control surgery+ Pembro group outcomes are significantly worse than known contemporaneous control groups like Checkmate 143 (OS of 10 months). How is this explained? Would be good to compare their results to known trials in this setting to see if this is applicable.\nOur response: All patients enrolled in our Phase 2 trial were at least at their second recurrence, with 10–15% at their third recurrence, whereas CheckMate 143 included only patients at first recurrence [6]. Consequently, our cohort was more advanced and more heavily pretreated, which likely explains the shorter mOS observed in the NLS+PEM arm compared with the nivolumab arm in CheckMate 143. We consider the most appropriate comparison to be the post-surgical (adjuvant) PEM arm in the neoadjuvant pembrolizumab study [7, 8], which enrolled patients at first or second relapse. In that study, the adjuvant PEM arm had an mOS of 6.4 months, while our control arm showed an mOS of 5.3–6.2 months (PP vs ITT populations), indicating broadly comparable outcomes despite 15% of our patients being treated at a third, more advanced recurrence. We added this explanation to the Discussion section (Lines 614–634).\n2. The OS of ITT group is very promising but the OS is nearly 10 months different from the per protocol group. An explanation of this difference would be helpful. If there is a big outlier in the ITT group, it would be preferable to just put the protocol OS in the abstract so it is clear what the benefit of pembro and LITT is, without confounders.\nOur response: Thank you for this important observation. In the ITT population, 5 patients (2 in NLS+PEM and 3 in LITT+PEM) were included who were not evaluable in the per-protocol analysis because they did not complete protocol-defined treatment and/or have confirmed grade 4 astrocytoma. Given the small sample size, inclusion of these additional patients had a substantial impact and increased the mOS and mPFS estimates in both arms, which explains the ~10-month difference between the ITT and per-protocol results. There was no single extreme outlier driving the ITT findings. For transparency, we report both ITT and per-protocol outcomes and have clarified in the text (and abstract) that the per-protocol analysis represents the primary estimate of treatment effect, with the ITT analysis presented as a complementary sensitivity analysis. We also add this clarification to the Results section (Lines 54–59; 253–257).\n3. MRIs of the partial and stable response and treatment failures would be helpful. No imaging data or case data is presented in the main figures\nOur response: We have now added representative imaging data to Figure 2, panel c. Specifically, we include contrast-enhanced MRI examples illustrating (i) early progressive disease in a control-arm patient, (ii) prolonged stable disease, and (iii) a partial response as defined by RANO criteria.\n4. A discussion on the use and picking of the Immune checkpoint inhibitor would be helpful as well. Is there a use for newer agents or PD-1, CTLA4, LAG3 inhibitors?\nOur response: Thank you for this helpful suggestion. We have expanded the Introduction (Lines 118–122) and the Discussion (Lines 542–551) sections to clarify our choice of immune checkpoint inhibitor and to comment on newer agents. The study was designed in 2014, when CTLA-4 and PD-1/PD-L1 were the principal checkpoint pathways in clinical use or active clinical testing. We selected PD-1 blockade because several reports at that time showed that a substantial fraction of high-grade astrocytoma express PD-L1, that PD-L1 expression increases with WHO grade, and that PTEN loss (common in high-grade astrocytoma) further upregulates PD-L1 [9-13]. On this basis, and in consultation with Merck, pembrolizumab was chosen as the study drug; its every-3-week dosing schedule also reduced visit burden compared with every-2-week PD-1 regimens such as nivolumab and was considered favorable for patient adherence in this recurrent population. In contrast, ligands for CTLA-4 (CD80/CD86) are largely restricted to cells in the tumor microenvironment rather than tumor cells themselves, and subsequent data from CheckMate 143 demonstrated that adding ipilimumab to nivolumab did not improve efficacy but increased toxicity [14], so we did not pursue CTLA-4 combination therapy. At the time of trial development, newer checkpoints such as LAG-3 and TIM-3 were either in very early preclinical/clinical stages for brain tumors or not yet available. Looking ahead, if IDO-1 downregulation in T cells is validated as a key mechanism for durable responses to LITT+PEM, combining anti–PD-1 therapy with an IDO-1 inhibitor would be a rational strategy to explore in future trials (Lines 599–602).\nReviewer #3 (Remarks to the Author):\n1. Unplanned interim analysis: The original protocol did not include an interim analysis. Conducting an unplanned interim analysis can inflate the type I error rate because no pre-specified boundaries were defined. Without a pre-specified alpha-spending plan, it is unclear what p-value threshold was used to justify early stopping.\nOur response: We agree with this concern and now explicitly acknowledge that the interim analysis was not pre-specified and therefore does not control the overall type I error rate. The original protocol planned 1:1 randomization between MK-3475 alone (NLS+PEM) and MK-3475 after MLA (LITT+PEM). As external data emerged suggesting limited efficacy of post-surgical pembrolizumab alone (NLS+PEM), the DSMC requested an unscheduled interim review for ethical reasons. At that review, descriptive analyses showed higher PFS and OS in the LITT+PEM arm compared with NLS+PEM arm: in the LITT+PEM arm (2 Phase 1 and 10 Phase 2 patients), 6-, 9-, and 12-month PFS probabilities were 45.0%, 45.0%, and 33.8% with a median PFS of 5.49 months, and 6-, 9-, and 12-month OS probabilities were 75.0%, 62.5%, and 50.0% with median OS not estimable. In the NLS+PEM arm, the 3-month PFS probability was 44.4% with a median PFS of 2.73 months (95% CI: 1.41–3.88), and 3- and 6-month OS probabilities were 83.3% and 31.3% with a median OS of 4.77 months (95% CI: 1.74–8.78). Based on the magnitude and consistency of these differences, together with evolving external evidence, the DSMC recommended stopping randomization and converting the trial to a single-arm LITT+PEM study. No formal alpha-spending plan or pre-specified p-value boundary was used; the decision relied on DSMC judgment in the spirit of minimizing harm, potential or perceived, to subjects and to complete the study in a timely manner. We have clarified this in the Methods (Lines 958–963; 1031–1040; 1056–1065) and Discussion (Lines 605–612), noting that the unplanned interim analysis is a limitation and that these results should be interpreted as hypothesis-generating. Full details of the interim analysis and statistical output are provided in the provided amended protocol (pp. 58–59).\n2. Potential instability of early findings: Early trends in treatment effect are often unstable. The apparent superiority observed at an unplanned interim analysis may not have persisted had the trial continued, particularly when a liberal or post hoc boundary was used for decision-making.\nOur response: We share this concern and, as outlined in our response to Comment 1, we now explicitly acknowledge that the interim analysis was unplanned, without a pre-specified alpha-spending boundary, and therefore the apparent early treatment effect should be regarded as hypothesis-generating and potentially unstable over time. Full details of the interim analysis and statistical output are provided in the amended protocol (pp. 58–59). To partially address the stability of the survival signal, we divided the per-protocol Phase 2 population into three cohorts: randomized NLS+PEM, randomized LITT+PEM, and non-randomized LITT+PEM (treated as an expansion cohort of the randomized LITT+PEM arm), and present their Kaplan–Meier curves together in the new Figure 4. The curves for the randomized and non-randomized LITT+PEM cohorts largely overlap and remain visually separated from the NLS+PEM control curve, suggesting that the survival advantage is consistent across early and later LITT+PEM patients. A similar approach of comparing randomized and subsequent expansion cohorts using overlaid K–M curves has been employed in a neoadjuvant pembrolizumab study in GBM published in Nature Communications [8], which we now cite as precedent. Nevertheless, we emphasize in the Discussion that, given the unplanned interim analysis and modest sample size, these findings require confirmation in larger, prospectively randomized trials.\n3. Mismatch between amended design and final analysis: After the interim analysis, the trial was amended to a single-arm design, referencing a historical median PFS of 2.9 months. The amended design and sample size calculation were based on a one-sample log-rank test. However, in the manuscript’s analysis, two-sample log-rank tests were applied to both the per-protocol and randomized populations. This represents a clear inconsistency between the design and the analysis. If the analysis method differs from the design basis, it becomes uncertain whether the nominal power and type I error rate are still properly controlled.\nOur response: We agree with this point and have revised the manuscript to better explain how the amended design relates to the final analyses. After the DSMC-mandated conversion to a single-arm study, the amended protocol specified a one-sample log-rank comparison of PFS in the LITT+PEM cohort against a historical median PFS of 2.9 months and used this framework for the sample size calculation. In the revised Statistical Methods and Results, we now explicitly describe this historical 2.9-month PFS as the design benchmark and state that, in the final analysis, the PFS distribution of the per-protocol LITT+PEM cohort was compared descriptively to this benchmark using Kaplan–Meier estimates and their 95% confidence interval, without performing a formal one-sample log-rank hypothesis test (Lines 172–179; 1031–1045).\nThe two-sample log-rank tests comparing LITT+PEM with NLS+PEM are now explicitly described as exploratory, secondary analyses intended to contextualize the magnitude of effect within the randomized portion of the trial; we note that these post hoc comparisons were not part of the amended design and do not preserve the nominal type I error. As described in our response to Comment 2, we also include an overlaid Kaplan–Meier plot of three cohorts (randomized NLS+PEM, randomized LITT+PEM, and non-randomized LITT+PEM) in the new Figure 4 to illustrate the consistency of the survival signal across early and later LITT+PEM patients.\n4. Was any analysis performed on intratumoral PD-L1 expression—either comparing expression levels between the two groups or evaluating the association between PD-L1 expression and PFS/OS? This analysis was described in the statistical analysis plan but was not reported or discussed in the manuscript.\nOur response: Thank you for highlighting this, and we apologize for not describing our PD-L1 analyses in the initial submission. As specified in the statistical analysis plan, we intended to compare intratumoral PD-L1 expression between treatment arms and to explore its association with PFS/OS. In practice, this was not feasible because tumor tissue was very limited, particularly in the LITT+PEM arm, where patients underwent biopsy only. Many biopsies were small, with scant viable tumor and suboptimal quality, raising a high risk of under-sampling.\nGiven these constraints, we prioritized assays according to tissue requirements and sensitivity (DNA-based, then RNA-based, then protein-based). All usable samples after clinical diagnostic needs were processed for whole-exome sequencing; of these, only 5 tumors (3 NLS+PEM and 2 LITT+PEM) had sufficient material left for bulk RNA-seq, and none had adequate remaining tissue for PD-L1 IHC as originally planned. In these 5 RNA-seq samples, we did not observe statistically significant differences in PD-L1 expression between arms (new Supplemental Figure S3c), but the sample size is far too small for any definitive conclusion. We have now added a description of these limitations and limited results to the Results (lines 479–483), Discussion/Limitations (Lines 643–647), and Methods (Lines 1067–1073) sections, and clarify that the planned PD-L1–survival correlation analyses could not be performed. Of note, recent reports have suggested that PD-L1 expression in the TME has only a modest correlation with response to anti-PD-1 inhibitor [15, 16].\n5. When you consider using UMAP to perform feature extraction, have you considered using some other common methods such as t-SNE? In addition, there might be some strong or weak associations, and as you have used the “FindNeighbors” function, I wish to know if you have found some structural clustering from the original sequencing. For example, some genes are highly correlated so that they potentially “work together”. This is additional advice for further analysis. You may provide some insight if it takes too much time to re-run the analysis.\nOur response: Thank you for this thoughtful suggestion. We do routinely examine alternative non-linear dimensionality reduction methods, including t-SNE. For this study, we selected UMAP because the low-dimensional embedding is used not only for visualization but also as a feature space for downstream quantitative analyses (e.g., computing Earth Mover’s Distance to track subtype changes over time), where preservation of both local and broader-scale structure is important.\nTo substantiate this choice, we performed a systematic comparison of UMAP and t-SNE embeddings for all PBMCs, T cells, and non-classical monocytes using the same clustering and several quantitative metrics (within-cluster compactness, between/within distance ratio, silhouette score, kNN purity, and 2D kNN-graph modularity). Across all three datasets, UMAP produced markedly more compact clusters, equal or better relative separation between clusters, higher silhouette scores and kNN purity, and higher modularity, indicating better alignment between the 2D geometry and biologically defined clusters. These results are summarized in a new supplemental document for review only, entitled “Comparison of UMAP and t-SNE embeddings.” We also added a statement in the Method section the rationale of using UMAP (lines 1134–1138).\nRegarding the “FindNeighbors” function: In our current analysis, the FindNeighbors function was used to construct a cell-cell graph and identify clusters of transcriptionally similar cells, which we then used for cell type annotation as presented in the manuscript. This procedure does not directly assess gene-gene correlations or infer gene interaction networks; rather, it operates at the level of cells. We agree that investigating structural relationships among genes (e.g., co-expression or gene-gene interaction patterns) could provide additional biological insight and could a valuable direction for future extensions of this study.\n6. For the Cox proportional hazard model that predicts the GSEA gene set, is there a mediation pathway for these genes to work? These gene sets may affect a specific biomarker, and further these biomarkers affect the survival probability. If so, would you be able to perform a mediation analysis?\nOur response: In our Cox proportional hazards model, the predictor is a composite activity score derived from the top 500 genes associated with non-classical monocyte activation (Figure 5c–f), rather than a single GSEA pathway or individual biomarker. Conceptually, our biological model is that LITT induces an early shift in non-classical monocytes, and this activation program in turn relates to survival; in that sense, the monocyte program can be viewed as a putative mediator of the LITT–survival relationship. However, given the small cohort size and the strong collinearity between different summaries of this program, we did not feel that a formal causal mediation analysis with explicit estimation of direct and indirect effects of treatment on survival would be stable or interpretable.\nTo address the spirit of your question, i.e., through what biology the 500-gene signal may act, we performed an exploratory pathway-level decomposition of this signature. Starting from the same 500 genes shown in Figure 5c, we annotated them into five functional pathways (cell adhesion, T-cell differentiation, macrophage differentiation, metabolic processes, and inflammatory responses) based on known gene functions and derived a pathway activity score for each patient. In age-adjusted Cox models that included one pathway score at a time, three of the five pathways (cell adhesion, metabolic processes, and inflammatory responses) each reproduced the survival association observed with the full 500-gene score, whereas T-cell differentiation showed a similar but non-significant trend and macrophage differentiation was not associated with outcome (updated Figure 5f).\nWe then asked whether these pathway-level signals contributed information beyond the broader 500-gene program. To do this, we performed a formal mediation analysis using Cox-based causal mediation and Fisher-exact A–M tests. Binary pathway scores were strongly associated with the 500-gene exposure and behaved as prognostic markers for OS when considered alone with age. However, once the global 500-gene signature was included, the pathways no longer provided robust independent prognostic information, and the estimated proportion of the effect mediated by any single pathway was small and non-significant. In the special case of the cell-adhesion pathway, the mediator and exposure are perfectly correlated, so the pathway could potentially stand in for the full 500-gene signature, but we cannot statistically separate their effects in this small cohort. Together, these findings suggest that the pathways are biologically meaningful sub-modules of the 500-gene program rather than distinct mediating mechanisms. The detailed description of the formal mediation analysis is attached as a supplemental file for review only, entitled “Mediation Analysis.”\n7. I am a bit confused with Figure 6 b, seems like only “pre-procedure to C1”, and “C2 to C4” is shifting to memory precursor, and “C1 to C2” have some shift to exhausted and effector memory, would you be able to comment on these three.\nOur response: The pattern you note in Fig. 6b (now new Figure 7b) — with pre-procedure→C1 and C2→C4 transitions predominantly shifting toward memory-precursor states, whereas C1→C2 shows more movement into exhausted and effector-memory states — is consistent with a dynamic, wave-like response of CD8⁺ T cells to antigen release and PD-1 blockade. After LITT, tumor antigens are released into the circulation in intermittent “pulses” rather than as a single bolus, as we observed previously when measuring brain-resident antigens following LITT [17]. Each pulse of antigen exposure in theory can recruit and activate both pre-existing memory T cells and newly generated clones.\nWith each subsequent cycle of pembrolizumab, these antigen-experienced T cells expand further. Previously exhausted cells in particular can be reinvigorated and transiently acquire effector or effector-memory phenotypes, which explains the increased flux from central/stem-like memory into effector and effector-memory clusters in the C1→C2 interval. As antigen levels decline, a proportion of these cells then contract back toward memory and memory-precursor states, which is what we predominantly observe in the pre-procedure→C1 and C2→C4 transitions. Thus, the three intervals you highlight likely capture different phases of repeated activation–expansion–contraction cycles.\nUsing our Earth mover’s distance metric–based approach, we can quantify the net direction of these shifts among T-cell clusters. On average, in the LITT+PEM long-term survivors, CD8⁺ T cells cycled more between memory and effector states and transitioned less frequently into the exhausted state than in the shorter-surviving cohort across corresponding time points.\nSpecifically, during the C1→C2 transition, arrows from central memory, effector-memory precursors, and stem-like cells tended to terminate within effector/effector-memory clusters in long-term survivors, whereas in short-term survivors, arrows from central memory and stem-like clusters preferentially converged on the exhausted cluster. Similar repeated waves of T-cell expansion and return to memory, coupled to ongoing TME activation and anti-PD-1 therapy, have also been reported with other TME-activating modalities, such as Tumor Treating Fields plus pembrolizumab [18] and with immune checkpoint blockade in other tumor types [19, 20], although those prior studies did not directly analyze cell-state transitions as we do here.\nReviewer #4 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nOur response: Thank you for your time and helpful comments.\nReferences for the responses to reviewer’s comments:\n1. Franceschi, E., et al., Phase II trial of carboplatin and etoposide for patients with recurrent high-grade glioma. Br J Cancer, 2004. 91(6): p. 1038-44.\n2. Reardon, D.A., et al., Phase I trial of irinotecan plus temozolomide in adults with recurrent malignant glioma. Cancer, 2005. 104(7): p. 1478-1486.\n3. Brandes, A.A., et al., How effective is BCNU in recurrent glioblastoma in the modern era? A phase II trial. Neurology, 2004. 63(7): p. 1281-4.\n4. Kappelle, A.C., et al., PCV chemotherapy for recurrent glioblastoma multiforme. Neurology, 2001. 56(1): p. 118-20.\n5. Chua, S.L., et al., Phase 2 study of temozolomide and Caelyx in patients with recurrent glioblastoma multiforme. Neuro Oncol, 2004. 6(1): p. 38-43.\n6. Reardon, D.A., et al., Effect of Nivolumab vs Bevacizumab in Patients With Recurrent Glioblastoma. JAMA Oncology, 2020. 6(7): p. 1003.\n7. Cloughesy, T.F., et al., Neoadjuvant anti-PD-1 immunotherapy promotes a survival benefit with intratumoral and systemic immune responses in recurrent glioblastoma. Nature Medicine, 2019. 25(3): p. 477-486.\n8. McFaline-Figueroa, J.R., et al., Neoadjuvant anti-PD1 immunotherapy for surgically accessible recurrent glioblastoma: clinical and molecular outcomes of a stage 2 single-arm expansion cohort. Nature Communications, 2024. 15(1): p. 10757.\n9. Wintterle, S., et al., Expression of the B7-related molecule B7-H1 by glioma cells: a potential mechanism of immune paralysis. Cancer Res, 2003. 63(21): p. 7462-7.\n10. Wilmotte, R., et al., B7-homolog 1 expression by human glioma: a new mechanism of immune evasion. Neuroreport, 2005. 16(10): p. 1081-5.\n11. Parsa, A.T., et al., Loss of tumor suppressor PTEN function increases B7-H1 expression and immunoresistance in glioma. Nat Med, 2007. 13(1): p. 84-8.\n12. Yao, Y., et al., B7-H1 is correlated with malignancy-grade gliomas but is not expressed exclusively on tumor stem-like cells. Neuro-Oncology, 2009. 11(6): p. 757-766.\n13. Bloch, O., et al., Gliomas Promote Immunosuppression through Induction of B7-H1 Expression in Tumor-Associated Macrophages. Clinical Cancer Research, 2013. 19(12): p. 3165-3175.\n14. Omuro, A., et al., Nivolumab with or without ipilimumab in patients with recurrent glioblastoma: results from exploratory phase I cohorts of CheckMate 143. Neuro-Oncology, 2017. 20(5): p. 674-686.\n15. Davis, A.A. and V.G. Patel, The role of PD-L1 expression as a predictive biomarker: an analysis of all US Food and Drug Administration (FDA) approvals of immune checkpoint inhibitors. J Immunother Cancer, 2019. 7(1): p. 278.\n16. Lu, S., et al., Comparison of Biomarker Modalities for Predicting Response to PD-1/PD-L1 Checkpoint Blockade: A Systematic Review and Meta-analysis. JAMA Oncology, 2019. 5(8): p. 1195-1204.\n17. Leuthardt, E.C., et al., Hyperthermic Laser Ablation of Recurrent Glioblastoma Leads to Temporary Disruption of the Peritumoral Blood Brain Barrier. PLoS One, 2016. 11(2): p. e0148613.\n18. Chen, D., et al., Efficacy and safety of adjuvant TTFields plus pembrolizumab and temozolomide in newly diagnosed glioblastoma: A phase 2 study. Med, 2025.\n19. Wang, K., et al., Combination anti-PD-1 and anti-CTLA-4 therapy generates waves of clonal responses that include progenitor-exhausted CD8(+) T cells. Cancer Cell, 2024. 42(9): p. 1582-1597.e10.\n20. Huang, A.C., et al., T-cell invigoration to tumour burden ratio associated with anti-PD-1 response. Nature, 2017. 545(7652): p. 60-65.\nPoint-by-point response to the reviewers’ comments (NCOMMS-25-76760B)\nWe thank the reviewers for their careful evaluation and enthusiasm for the manuscript. Reviewers’ comments are shown in black, with our responses in blue.\nReviewer #1 (Remarks to the Author):\nNo comments provided\nOur response: Thank you for your time.\nReviewer #2 (Remarks to the Author):\nThe authors have successfully addressed all of my comments. I believe this is worthy of publication. Congratulations.\nOur response: Thank you for your time and helpful comments.\nReviewer #3 (Remarks to the Author):\nThe authors gave extensive responses and addressed our concerns.\nReviewer #3 (Remarks on code availability):\nReadME file is easy to understand and follow. I didn't install or run the code since it is too time consuming.\nOur response: Thank you for your time and helpful comments.\nReviewer #4 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nReviewer #4 (Remarks on code availability):\nI have reviewed the updated code without running everything due to the anticipated time and effort needed. I have reviewed the structure of the guidelines for the analysis, which is well-written and easy to follow. I have reviewed most of the R scripts, and there are no obvious mistakes. The only suggestion from me is that the author may consider providing a small diagram in the GitHub page that illustrates the overall working process (trial design, data processing, sequencing, then UMAP, survival analysis, etc.), like what they did in the manuscript. This is not required for publication, just for the ease of future readers who wish to use the method.\nOur response: Thank you for your time and helpful suggestion. We have added a diagram in the GitHub page as suggested."}, {"round": 2, "reviewer_comments": {"reviewer_2": "The authors have successfully addressed all of my comments. I believe this is worthy of publication. Congratulations.", "reviewer_3": "The authors gave extensive responses and addressed our concerns.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": ""}]

Decision letter