Article https://doi.org/10.1038/s41467-024-47185-9 Microenvironmental reorganization in brain tumors following radiotherapy and recurrence revealed by hyperplexed immunofluorescence imaging Spencer S. Watson1,2,3,4 , Benoit Duc1,2,3,4, Ziqi Kang 5, Axel de Tonnac5, Nils Eling 6,7, Laure Font 1,8, Tristan Whitmarsh9, Matteo Massara1,2,3,4, iMAXT Consortium*, Bernd Bodenmiller 6,7, Jean Hausser 5 & Johanna A. Joyce 1,2,3,4,10 The tumor microenvironment plays a crucial role in determining response to treatment. This involves a series of interconnected changes in the cellular landscape, spatial organization, and extracellular matrix composition. How- ever, assessing these alterations simultaneously is challenging from a spatial perspective, due to the limitations of current high-dimensional imaging techniques and the extent of intratumoral heterogeneity over large lesion areas. In this study, we introduce a spatial proteomic workflow termed Hyperplexed Immunofluorescence Imaging (HIFI) that overcomes these lim- itations. HIFI allows for the simultaneous analysis of > 45 markers in fragile tissue sections at high magnification, using a cost-effective high-throughput workflow. We integrate HIFI with machine learning feature detection, graph- based network analysis, and cluster-based neighborhood analysis to analyze the microenvironment response to radiation therapy in a preclinical model of glioblastoma, and compare this response to a mouse model of breast-to-brain metastasis. Here we show that glioblastomas undergo extensive spatial reor- ganization of immune cell populations and structural architecture in response to treatment, while brainmetastases show no comparable reorganization. Our integrated spatial analyses reveal highly divergent responses to radiation therapy betweenbrain tumormodels, despite equivalent radiotherapy benefit. Glioblastoma is the most common and most aggressive primary brain tumor in adults, with an incidence rate of up to 5 per 100,000 people1. Standard-of-care treatment, consisting of surgical resection, ionizing radiation (IR), and temozolomide-based chemotherapy, leads to only a transient response. Consequently, median survival is just over 14 months, and the five-year survival is <5%2,3. This poor prognosis is due to tumor recurrence in nearly all cases, emphasizing the need for a deeper understanding of the mechanisms driving therapeutic resistance. Previous studies have investigated how glioblastoma responds and recurs following both IR therapy and macrophage-targeted immunotherapy4–10. One key mediator of resistance to treatment is tumor-cell heterogeneity and plasticity, resulting in IR-resistant subpopulations that can rapidly repopulate the lesion11–13. However, Received: 30 January 2023 Accepted: 22 March 2024 Check for updates A full list of affiliations appears at the end of the paper. *A list of authors and their affiliations appears at the end of the paper. e-mail: drspencerwatson@gmail.com; johanna.joyce@unil.ch Nature Communications | (2024) 15:3226 1 12 34 56 78 9 0 () :,; 12 34 56 78 9 0 () :,; http://orcid.org/0009-0003-2084-610X http://orcid.org/0009-0003-2084-610X http://orcid.org/0009-0003-2084-610X http://orcid.org/0009-0003-2084-610X http://orcid.org/0009-0003-2084-610X http://orcid.org/0000-0002-4711-1176 http://orcid.org/0000-0002-4711-1176 http://orcid.org/0000-0002-4711-1176 http://orcid.org/0000-0002-4711-1176 http://orcid.org/0000-0002-4711-1176 http://orcid.org/0009-0006-6648-7175 http://orcid.org/0009-0006-6648-7175 http://orcid.org/0009-0006-6648-7175 http://orcid.org/0009-0006-6648-7175 http://orcid.org/0009-0006-6648-7175 http://orcid.org/0000-0002-6325-7861 http://orcid.org/0000-0002-6325-7861 http://orcid.org/0000-0002-6325-7861 http://orcid.org/0000-0002-6325-7861 http://orcid.org/0000-0002-6325-7861 http://orcid.org/0000-0002-3321-7199 http://orcid.org/0000-0002-3321-7199 http://orcid.org/0000-0002-3321-7199 http://orcid.org/0000-0002-3321-7199 http://orcid.org/0000-0002-3321-7199 http://orcid.org/0000-0002-6332-2598 http://orcid.org/0000-0002-6332-2598 http://orcid.org/0000-0002-6332-2598 http://orcid.org/0000-0002-6332-2598 http://orcid.org/0000-0002-6332-2598 http://crossmark.crossref.org/dialog/?doi=10.1038/s41467-024-47185-9&domain=pdf http://crossmark.crossref.org/dialog/?doi=10.1038/s41467-024-47185-9&domain=pdf http://crossmark.crossref.org/dialog/?doi=10.1038/s41467-024-47185-9&domain=pdf http://crossmark.crossref.org/dialog/?doi=10.1038/s41467-024-47185-9&domain=pdf mailto:drspencerwatson@gmail.com mailto:johanna.joyce@unil.ch another critical aspect is how the tumor microenvironment (TME) responds to IR, as this can modulate the efficacy of IR through several mechanisms14–16. We previously revealed how different treatment modalities induce substantial changes in immune cell populations in murinemodels of glioblastoma4–7, which are related to eventual tumor recurrence. However, due to the substantial number of cellular and extracellular interactions in the TME, and consequently the highly focused nature of previous studies, many of the general mechanisms driving these pro-survival roles remain uncertain. It is thus increasingly evident that there is a need to investigate the TME as a whole, to interrogate the interplay of all themajor cell types with each other and with structural features in the tumor ecosystem. One challenge in studying the entire TME in all its complexity is the substantial amount of different cell types, and the high number of markers needed to delineate and identify each cell type in a single sample. The typical methods used to explore cellular heterogeneity in the TME include multiparameter flow cytometry and single-cell RNA- seq. However, both methods disassociate the tissue, thereby losing all spatial context of cellular interactions, in situ cellular morphology, and correlation with tissue architecture and ECM. For this reason, high-dimensional imaging techniques have become increasingly vital to the study of the TME. Imaging mass cytometry (IMC)17, tissue cyclic immunofluorescence (t-CycIF)18, Multiplexed Ion Beam Imaging by Time of Flight (MIBI-TOF)19, MALDI mass spectrometry imaging (MALDI-IMS)20, and oligo-conjugated antibody cyclic immuno- fluorescence (Ab-oligo cyCIF)21 are among the current approaches22,23 for highly multiplexed tissue imaging capable of assessing 20–100 markers on a single tissue section. While these techniques represent powerful and well-validated methods of examining the spatial topography of the TME, they are not universally suited to all applications. IMC has relatively small imaging areas of 1mm2 which have a spatial resolution of approximately 1 µm (equivalent to ~10Xmagnification).MIBI-TOF can achievemuch higher resolutions, up to 200nm,but this comeswith a significant trade-off in lengthy image acquisition times. MALDI-IMS can easily image whole- slide sections, but with a maximum resolution of only 10 µm24. This limits the ability to study tumors with structural features that exceed these dimensions, or with cellular interactions involving small den- dritic extensions, for example. In addition, commercial mass cyto- metry imaging and Ab-oligo cyCIF solutions have a substantial initial adoption cost, and require custom-conjugated antibodies, both of whichcanbeprohibitively expensive, or challenging to conjugate if the required markers are not commercially available. A significant addi- tional limitation of these approaches is the low throughput, with typical automated workflows allowing imaging of <10 slides per week. Immunofluorescence (IF)-based multiplexed approaches use off-the- shelf commercially available antibodies, but have many of their own challenges. Current highly multiplexed imaging approaches require repeated rounds of staining, imaging, and stringent stripping or quenching of markers which can damage tissue samples and the pro- tein epitopes targeted by antibodies. Moreover, while several of these antibody removal methods are well validated in formalin-fixed paraf- fin-embedded (FFPE) tissue, they can damage or destroy lightly fixed frozen tissue. In this study, we present a complete workflow for whole-slide highly multiplexed tissue imaging, which we term Hyperplexed Immunofluorescence Imaging (HIFI), to addressmultiple limitations of existing methodologies. By requiring no bespoke reagents or equip- ment, using commonly available reagents, and utilizing open-source analysis approaches we aim to address the issues of cost and accessi- bility of high-dimensional imaging for a broad scientific community. We apply this workflow to the study of TME response to focalized IR in a geneticmousemodel of glioblastoma. We employ a 45-marker panel to analyze a sample set of tissues collected prior to treatment, 7 days post-IR, or at the point of tumor recurrence to study the spatial microenvironment response to treatment. To investigate the extent of microenvironmental response depending on the tumor type, we also compare the glioblastomamodel to abreast-to-brainmetastasismodel at comparable volumes and timepoints. This comparison reveals substantial TME reorganization in glioblastomas in response to IR treatment, and the consistent generation of survival-promoting spatial niches at 7 days post-IR. Conversely, brain metastases do not show spatial reorganization or significant debulking in response to IR, despite both tumor types receiving equivalent survival benefit from the therapy. Results Overview of HIFI spatial analysis workflow We purposely designed HIFI to be straightforward to implement with open-source software and without the need for specialized lab equip- ment. For these reasons, all IF staining was performed manually using standard benchtop methodology and reagents, and all imaging was performed with conventional commercial slide-scanning microscopes. Figure 1 provides a graphical overview of the typical HIFI workflow, showing the process of cyclic immunofluorescence imaging, whole- slide image alignment and registration, machine-learning structural annotation, deep-learning cell segmentation, and clustering-based cell classification to generate highly annotated digital pathology images for spatial analysis. The workflow time course for a single experiment is dictated by the size of the antibody panel, how many markers can be multiplexed in each imaging round, and the speed of image acquisition. Multiple days are required to complete a full hyperplexed panel. However, as most modern slide-scanning fluorescence microscopes have capa- cities for 80–100 slides, the ability to scale up to high-throughput offsets the total time required for an experiment when compared to other highly multiplexed approaches. 40× and 60× magnifications across tissue sections are also feasible, however, this would come with a tradeoff of increased scanning time. A core component of the HIFI variant of cyclic-IF that we present here is the gentle antibody removal, coupled with the avoidance of crosslinkingbetweenfluorophores and tissue in the imaging stage. The most disruptive step in other cyclic-IF approaches is typically the elution of marker signal between imaging rounds. These methods use either stringent acidic buffers to strip antibodies25, or a combination of light and hydrogen peroxide to inactivate fluorophores18, both of which can potentially damage lightly-fixed and fragile tissue samples. HIFI utilizes a thiol-based elution buffer to reduce the disulfide bonds of bound antibodies, releasing the antibodies from the tissue without damaging tissue integrity. Thiol-based elution permits efficient yet gentle removal of standard off-the-shelf primary and secondary anti- bodies. Elution was further optimized for tissue sections by using an appropriate pH and a brief 3-min elution step. This method was pre- viously employed in cell culture using the iterative indirect immuno- fluorescence imaging (4i) approach26, and was adapted herein for use in tissue sections for the first time. The basis of 4i is the observation that high-energy light from fluorescencemicroscopy light sources can result in indirect oxidative crosslinking of fluorophores to the tissue, preventing efficient elution. Inhibition of this crosslinking through oxygen radical scavengers, or reduced light energy, can thereby improve antibody elution at less stringent conditions. To address this, we purposely optimized HIFI for low-power LED fluorescent light sources, to prevent crosslinkingwhilemaintaining robust fluorescence signals. These optimizations additionally allow for the use of standard glass coverslips and glycerol-based mounting media while imaging, rather than specialized mounting solutions. A further benefit of low- power imaging is reduced fluorescence spillover between channels when combined with appropriate bandpass filters. The workflow of immunostaining, imaging, coverslip removal, and elution are repeated until all rounds of markers are imaged. Article https://doi.org/10.1038/s41467-024-47185-9 Nature Communications | (2024) 15:3226 2 Following this non-destructive imaging, samples can be stained with hematoxylin and eosin (H&E) for future analysis and stored long-term, or used for additional experiments, such as proteomics. The post- processing, alignment, and analysis of HIFI-generated images are dis- cussed in the following use cases. Generation of IR treated brain tumor samples and antibody labeling panel Glioblastomas present with a high degree of heterogeneity, with spa- tial features that can cover substantial distances. Because of this, selecting individual regions of interest (ROIs) for high-dimensional image analysis can impart a substantial potential for selection bias. For this reason, to robustly investigate how the cellular and structural topography of glioblastoma responds to IR therapy, we endeavored to develop a high-dimensional spatial analysis pipeline for whole tissue sections that was not limited by sample size or preparation, and which could achieve subcellular resolution and semi-quantitative protein detection. We utilized the RCAS-hPDGF-B; Nestin-Tv-a; Ink4a/Arf KO geneti- cally engineeredmousemodel (GEMM) tomodel glioblastoma27,28. This model drives neoplastic transformation in nestin-expressing neural progenitor cells via overexpression of platelet-derived growth factor-B (Pdgfb), resulting in glioblastomas that mimic the human proneural phenotype, and with a fully intact immune system27,28. In addition, the transformation includes expression of green fluorescent protein (GFP) in Pdgfb-overexpressing cells, to fluorescently label and track the resulting tumor cells. These gliomas are termed PDGfp herein. For a comparative dataset we selected an orthotopic immune-competent model of breast-to-brain metastasis (BrM) utilizing a luminal HER2 + MMTV-PyMT-derived cell line5,29 injected at matched cranial coordi- nates as used for initiation of the PDGfp model. PDGfp or BrM tumors were initiated in mice at 5-6 weeks of age, and animalsmonitoredweekly bymagnetic resonance imaging (MRI) to screen for tumor formation and growth (Fig. 2a–c). PDGfp tumors were initiated in both male and female mice, while breast-BrM were injected only in femalemice. Once PDGfp tumors reached a volume of >20mm3, they were either harvested (n = 5), or treated with a single focalized 10Gy IR dose (n = 10). BrM were similarly monitored until tumors reached comparable volumes, and then either harvested (n = 3), or treatedwith a single focalized 15Gy IRdose (n =6) (Fig. 2d). Appropriate radiation doses for both tumor types were based on previous data that assessed efficacious IR doses in these models (ref. 5, and Wischnewski, [..], Joyce, manuscript in preparation). Critically, despite the lesser tumor volume reduction in BrM, both models ultimately showed equivalent survival benefit from these focalized IR doses. Comparisons between treated samples in this study and historical data for untreated subjects were completely in keeping with previous survival data for these models4,5,29 (Wischnewski, [..], Joyce, manuscript in preparation). At 7 days post-IR, MRI was performed on all mice, and tumor-bearing brains were harvested for n = 5 PDGfp and n = 3 BrM mice. MRI mon- itoring was continued bi-weekly for the remaining mice until tumor recurrenceswere detected. Brainswere subsequently harvested (PDGfp n = 5, BrM n = 3) and tissues embedded for cryo-sectioning. Volume measurements showed substantial reductions in tumor size in PDGfp tumors 7 days post-IR treatment, and a subsequent increase in volume at the point of tumor regrowth (Fig. 2b, c). BrM tumors showed more modest reductions in tumor volumes in response to IR, and a similar increase in volume at the point of regrowth (Fig. 2d). HIFI Cyclic- Immunofluorescence High-Dimension Whole-Slide IF Machine Learning Region Classification and Deep Learning Cell Segmentaiton Whole Tiled Image Alignment, Registration, and Autofluorescence Subtraction 1 mm Semi-Supervised Cell Classification Computer-Vision Digital Pathology Images Cell Community and Spatial Network Analysis Regional Cell Quantification Perivascular Niche Astrocyte_B Fibroblast TAM_A T_Cell Tumor_A 0% 20% 40% 60% Pe rc en tT ot al C el ls in R eg io n Astrocyte_B Neuron_A TAM_A Tumor_A Tumor_C 0% 20% 40% 60% Pe rc en tT ot al C el ls in R eg io n Untreated 7 Day Regrowth Tumor Leading Edge −10 −5 0 5 −10 −5 0 5 UMAP1 U M A P 2 Fig. 1 | Overview of Hyperplexed Immunofluorescence Imaging (HIFI) Work- flow. Experimental workflow of cyclic immunofluorescence staining, followed by image processing, alignment, and registration to create 45+ dimensional images across whole-slide sections. Specific domains within HIFI images of tumors were automatically annotated using trained machine-learning classifiers, and individual cells were segmentedwith deep-learning object detection. Single-cell objects were annotated as individual cell types with semi-supervised classification and mapped back onto images to create highly-annotated digital pathology images. Images were analyzed for regional cellular composition and spatial organizational analysis. The HIFI workflow is scalable to over 100 simultaneous sections for high- throughput spatial experiments. Article https://doi.org/10.1038/s41467-024-47185-9 Nature Communications | (2024) 15:3226 3 Following cryo-sectioning of brain tissue samples, at least three replicate sections spaced at regular intervals across the entire Z-depth of each sample were used to account for intratumoral heterogeneity. In total, 87 whole-brain tissue sections were prepared for the imaging pipeline. 45 markers were chosen to interrogate the PDGfp tumor samples at each treatment stage. The antibody panel was based on previous studies characterizing the TME of murine glioma models by flow cytometry and other methods4–6,30. Various markers were selected to label each of the predominant tumor, glial, vascular, and immune cell types (Table 1). In addition, multiple antibodies were selected to label ECM structures relating to the healthy brain, perivascular niche, and treatment-inducedfibrosis. An additional series ofmarkerswas used for the purpose of training machine learning models and for quality con- trol. A subset of 29 of these markers was selected to interrogate BrM samples, to which 4 new antibodies were added to specifically analyze breast cancer heterogeneity (Supplementary Table 1). A critical component in designing theHIFI panel was to determine which antibodies could be multiplexed together, as well as the ability to reliably label their target epitope following multiple rounds of staining, imaging and elution. Antibody validation was performed using existing methodology (Supplementary Note 1)18,26, which entails comparisons of marker specificity in previously cyclic-stained tissue versus unstained sections. To verify that markers are removed between imaging rounds, control experiments were performed to repeatedly label and image multiple tissues using a set of primary and secondary antibodies. Markers were eluted following imaging, and tissues were re-stained with only secondary antibodies and reimaged. This process was repeated multiple times to determine single-cell mean fluorescence intensity (MFI) for markers at each stage (Supple- mentary Fig. 1). These data also show which marker intensities are maintained over multiple elution rounds, which decrease, and which even increase over rounds, thereby informing the sequential order of markers in the panel. Markers were thus arranged into 12 individual multiplexed panels based on their labeling efficiency and species cross-reactivity (Table 1, Supplementary Fig. 1a). Whole-slide image alignment and registration Raw tiled image data from all rounds were stitched together for each image with affine transformation (Zeiss Zen software package) to create seamless whole-section images. Background fluorescent signal was removed using the rolling-ball algorithm. In cyclic fluorescence slide scanning, each round of imaging may not necessarily be in the same exact position or orientation due to stage drift and slide placement, so multiplexed images cannot be simply merged directly. Whole-section images on this scale could not be aligned with publicly or commercially available tools due to a b c Untreated 7 Days Post-IR Tumor Regrowth e N = 5 N = 5 N = 5 10/15 Gy Focalized Radiotherapy Initiate tumors (5-6 weeks of age) MRI - 5 w - MRI 7 Days Post-IRUntreated Untreated 7 Days Post-IR Regrowth 0 20 40 60 80 100 120 Le si on Vo lu m e (m m 3 ) Untreated 7 Days Post-IR Regrowth 0 20 40 60 80 100 120 Le si on Vo lu m e (m m 3 ) Untreated 7 Days Post-IR Regrowth PDGfp BrM d Tumor Regrowth Untreated 7 Days Post-IR Tumor Regrowth N = 3 N = 3 N = 3 Fig. 2 | Irradiation Treatment Sample Collection and Antibody Panel. a PDGfp and BrM tumorswere initiated inNTVA-Ink4a/Arf−/− and C57Bl6mice respectively, andmonitored byMRI. Experimental endpoints (red arrows) are indicated for each group. Biweekly MRI monitoring tracked the process of regression and recurrence for each mouse following IR for b, c PDGfp (n = 5 mice for each treatment group) and d, e BrM tumors (n = 3 mice for each treatment group. Upon tumor detection, mice were grouped into cohorts and collected as either (i) ‘untreated’ samples, (ii) treated with 10Gy (PDGfp) or 15 Gy (BrM) focalized irradiation (IR) therapy and harvested 7 days post-IR, or (iii) treated with 10 or 15 Gy IR and harvested upon tumor regrowth. White arrows indicate the tumor in both tumor types, red arrows indicate post-IR lesion in PDGfp tumors. Source data are provided as a Source data file. Article https://doi.org/10.1038/s41467-024-47185-9 Nature Communications | (2024) 15:3226 4 memory limitations and array-size limits. Therefore, we had to develop an automated software tool in Python to align and register all rounds of multiplexed imaging into a single hyperplexed image. The “HIFI Alignment” tool utilizes existing algorithms for pyramidal sub-pixel image alignment based on pixel intensity31, and implements them in Pythonusing an image handling strategy tobypassmemory limitations due to image array-size, bit-depth, or the number of channels (Fig. 3a). This strategy enableswhole-slide rigid-body or affine alignment of 45+ channel images, which is only limited by the physical RAM avail- able on a workstation (Fig. 3b, c). Whole mouse brain sections were aligned using consumer-gradeworkstations with 256GBof RAM,while image subsets of whole PDGfp tumors could be aligned on standard computers with 64 GB of RAM. HIFI Alignment includes a CZI reader and writer for Axio Scan.Z1 images, but can also accommodate images in multiple TIFF formats. An important feature of the HIFI Alignment tool is the ability to perform pixel-by-pixel removal of tissue autofluorescence from all channels to prevent any contamination of marker signals in sub- sequent rounds of imaging. Imaging of tissue samples in all channels prior to antibody labeling, as was performed with the PDGfp and BrM IR-treated samples, creates an image of overall tissue autofluorescence which is then subtracted from respective channels in each round of imaging. Autofluorescence subtractionwas then performed for all BrM samples, but not for PDGfp samples - so as to not remove the endo- genous GFP signal in the tumor cells. The HIFI methodology was specifically developed to work with challenging sample types that perform poorly with other highly mul- tiplexed imaging techniques. To validate that HIFI also works in more standard sample types, we performed additional HIFI experiments across multiple murine organ FFPE samples (Supplementary Fig. 2). Following imaging, sections were cropped to the area of interest. Whole-slide imaging facilitated selection of regions for analysis that encompassed the entire tumor area and large regions of the sur- rounding brain, thereby eliminating sample biaswithin the tumor. This workflow generated sixty PDGfp images in 45-dimensions, and twenty- seven BrM images in 33-dimensions (Fig. 3d, e). Untreated and 7-days post-IR PDGfp tumors were found to be consistent across replicates in terms of position and morphology, while regrowth gliomas showed extensive morphological heterogeneity (Fig. 3f, Supplementary Fig. 3a). Conversely, BrM tumorswere observed tobe highly consistent across all replicates and treatment conditions (Fig. 3g, Supplemen- tary Fig. 3b). Machine learning structural annotation and deep learning cell segmentation HIFI images were first reviewed for all channels to identify any imaging aberrations, such as dust contamination or tissue deformations from sectioning, for the purpose of excluding these areas fromdownstream analysis (Fig. 4a). During this critical quality control review,HIFI images that did not meet stringent criteria were excluded from the study, resulting in n = 81 images for subsequent in-depth analysis. Images were analyzed with the open-source digital pathology suite QuPath32. The ability to image an entire tumor area with subcellular resolution facilitated the correlation of individual cell typeswith larger features of tumor architecture. To analyze this topographical heterogeneity, we trained an AI pixel-classifier model with pathological annotations based on multiple cellular and ECM features, and then applied each model for each feature across all images (Fig. 4a). Tissue regionswithin each sample were automatically annotated to identify lesion bound- aries, tumor nests, vasculature, areas of fibrosis, and modified to delineate perivascular niches, the surrounding brain, and tumor-brain interface border. Table 1 | PDGfp HIFI marker panel Tumor GFP Pan tumor Olig2 OPC Sox9 Glial precursor Sox2 Glial precursor Immune CD3 Pan T cell CD8a Cytotoxic T cell Iba1 Macrophage CD68 Macrophage P2YR12 Microglia CD45 Pan leukocyte CD206 Meningeal Ly6b Neutrophil S100A8 Neutrophil activation MPO Neutrophil activation ECM Tenascin-C Fibrotic ECM Laminin Vascular ECM Periostin Non-structural ECM CSPG5 Parenchymal ECM Fibronectin Fibrotic ECM Collagen I Fibrotic ECM Collagen IV Fibrotic ECM Neurons FoxP1 Neural stem cell NeuN Pan neuronal NF-H Neurofilament Vascular ER-TR7 Reticular fibroblast PDGFRB Pan perivascular cells CD31 Endothelial αSMA Mural cells CD13 Pericyte VE-Cad Endothelial Glial Vimentin Astrocyte subset GFAP Astrocyte activation Podoplanin Reactive gliosis S100B Pan astrocyte Desmin Astrocyte subset NG2 OPC Cell State Ki67 Proliferation CC3 Apoptosis HIF1a Hypoxia Lamin AC Nuclear envelope Functional RedDot2 QC WGA Machine learning αTubulin Machine learning Phalloidin Machine learning List of multiplexed markers used specifically for PDGfp glioblastoma samples, indicating cate- gory of marker, marker name, and marker target. Article https://doi.org/10.1038/s41467-024-47185-9 Nature Communications | (2024) 15:3226 5 Tiled 5-channel 20X cyclic fluorescence images for multiple rounds over multiple slides Tile alignment, fusion, rolling- ball background subtraction a b Tunable pixel-by-pixel autofluorescence subtraction Channel selection and metadata propagation Pyramidal sub-pixel rigid body and affine transformation generated from reference channel for every image round Subsequent images loaded in memory, transformation applied then dumped from memory c Image channels merged into single HIFI image, exported into multiple formats 45 dimensional, whole-slide, 20X DAPI-Col IV-GFP-CD31-GFAP-NeuN d e DAPI-Col I-GFP-CD31-GFAP-P2YR12 DAPI-Olig2-CD31-NeuN-NF-H-GFAP Untreated Regrowth7 Days Post-IR DAPI-Col I-GFP-CD31-GFAP-NeuN DAPI-TNC-Ecad-CD31-GFAP-P2YR12 DAPI-Ecad-CD31-NeuN-GFAP-CD13-CD3 f g Untreated Regrowth7 Days Post-IR DAPI-TNC-EpCAM-CD31-GFAP-NeuN Fig. 3 | HIFI Alignment. a Image processing of raw image data consisted of tile stitching and flat-field correction, followed by rolling-ball background subtraction. Scale bar = 2mm.b Panel depicts the repeated steps of the “HIFI Alignment” tool to align cyclic IF images while avoiding memory limitations and removing endogen- ous tissue autofluorescence. c Aligned and registered images were then merged into a single 45+ dimension whole-slide image. Images for d PDGfp and e BrMwere cropped to encompass the entire tumor region and surrounding brain. Scale bars = 500 µm. Representative images of (f) PDGfp and (g) BrM samples from untreated, 7 days post-IR, and regrowth tumors. Scale bars = 400 µm. Article https://doi.org/10.1038/s41467-024-47185-9 Nature Communications | (2024) 15:3226 6 a Binary classifiers created for tumor nests, fibrosis, vessels, and adjacent parenchyma Pixel classifiers trained from 25% of image data Region annotations applied to all images, tissue deformations then removed from all ROIs Tumor - Brain - Fibrosis - Vessels - Perivascular - Border b Ground Truth Annotations StarDist Model Glioma-Trained Model Watershed Segmentation Cytoplasm Expansion DAPI Cell Classification 567 Cells 567 Cells 371Cells 379 Cells 559 Cells c d BrMPDGfp Tumor_A Tumor_B Tumor_C OPC-Like Neuron_A Neuron_B Astrocyte_A Astrocyte_B CD45+ Fibroblast Neutrophil T_Cell Vessel_A Vessel_B Vessel_C TAM_A TAM_B TAM_C TAM_D Tumor_A Tumor_B Tumor_C Tumor_D Neuron_A Astrocyte_A Astrocyte_B CD45+ Fibroblast Neutrophil T_Cell Vessel_A Vessel_B TAM_A TAM_B TAM_C TAM_D Fig. 4 | Region Annotation, Cell Segmentation and Classification. a Binary pixel classifiers were trained for multiple spatial regions (as indicated) in QuPath on a subset of image data. Tissue deformations and imaging errors were removed from regions of interest (ROIs) during quality control checks (white arrow indicates a representative example of a tissue tear). Scale bar = 500 µm. b Nuclei in murine glioblastoma samples were manually annotated to train StarDist nuclear detection models. Detection accuracy was compared to threshold-based watershed segmentation and the publicly available StarDist model for immunofluorescence (DSB_Heavy). Segmented nucleiwere expanded by 2.5 µmto capture cell cytoplasm and classified by semi-supervised cell classification. Scale bars = 50 µm. Region annotation, cell segmentation, andcell classificationwere applied independently to all (c) PDGfp and (d) BrMHIFI images to generate fully annotated digital pathology images. Scale bars = 400 µm. Article https://doi.org/10.1038/s41467-024-47185-9 Nature Communications | (2024) 15:3226 7 We then performed cell segmentation for single-cell measure- ments of MFI, cell size, and localization. However, glioma tissue contains densely packed irregular nuclei, andDAPI staining of nucleic acids can additionally create non-uniform labeling across nuclei (Fig. 4b). This combination of factors dramatically reduced the accuracy of intensity-based watershed algorithms for nuclei seg- mentation in these tissues. (Fig. 4b). For this reason, we instead employed the supervised deep-learning algorithm StarDist33 imple- mented in QuPath. StarDist utilizes U-Net34, a convolutional neural network (CNN) designed for biomedical image analysis. As such, StarDist requires annotated ‘ground truth’ images to train the CNN to learn to accurately predict object probabilities. Ground truth nuclear annotations were manually performed on 144 large images of mul- tiple murine tissues stained with DAPI. Our murine training dataset included multiple samples from PDGfp tissues, breast-BrM, lung- BrM, breast tumors, human breast tumor xenografts, and healthy tissue from brain, mammary glands, and lungs. Additional images were generated via image transformations using transposition, flip- ping of the axes, as well as random changes of signal intensity, to encourage the CNN to be robust to these transformations. Transfer learning from the published StarDist immunofluorescence model33 was also performed to further increasemodel accuracy. We used this approach to generate two robust deep-learning nuclei segmentation models; one specifically trained on only PDGfp tumors for enhanced accuracy, and another trained on the entire data set for broad applicability to any mouse tissue (Fig. 4b). The publicly released DSB_Heavy StarDist cell segmentation model was found to have an accuracy of 68.86% on our training data, while our glioma-specific model had an accuracy of 74.96% (Fig. 4b). Accuracy was further enhanced in the dataset by filtering all cells with a detection prob- ability score of <0.5. Segmentation was performed with StarDist using our murine- trained model to identify each nucleus as a vector object, and each object was given an additional maximum expansion of 2.5 μm to measure cytoplasmicmarker signals (Fig. 4b).Cell size, shape,MFI, and XY location were measured for every cell. Further spatial measure- ments were performed following single-cell MFI and morphometric measurements. Additionally, the distance of each cell to the nearest region annotationborderwas recorded tomeasure cell proximity to all structural features within the image. Clustering based semi-supervised cell classification A potential drawback of using cyclic immunofluorescence for high- dimensional imaging is the low dynamic range, or depth, of marker signals compared to mass spectrometry-based approaches. This can create challenges for standard clustering-based approaches of unsu- pervised cell classification, especially when some markers have low intensity relative to other markers. A further challenge to employing computationally intensive clustering algorithms is that applying them to very large single-cell HIFI datasets (such as the PDGfp dataset con- taining ~7 × 106 cells) is the time required for analysis, which can be >1 week for a single dataset. Therefore, to achieve distinct clusters for cell classification despite low signal dynamic range, we next optimized a computationally efficient FlowSOM approach35. MFIs for markers unique to specific cell types were measured for the whole-cell area or the nuclear area, depending on the cell biology of each marker (Supplementary 2). This enhanced cluster specificity and reduced signal spillover to neighboring cells. Each MFI used for classificationwas scaled from0 to 1, and clipped to the 99.7 percentile. Cells were filtered based on size to remove fragments, and each PDGfp and BrM dataset were clustered independently across all treatment types with FlowSOM to maintain uniformity. FlowSOM was set to produce 100 unique nodes for each dataset. Hierarchical clustering of marker MFIs was used to cluster nodes into cell type annotations (Supplementary Fig. 4a–c). Cellular annotation of unbiased clusters based on previous bio- logical knowledge of the TME resulted in semi-supervised cell classi- fication (Fig. 4c, d, Supplementary Fig. 5a, b). To maintain aspects of cellular heterogeneity we used generic nomenclature, such as Tumor_A or Tumor_B, to provide general classifications while still maintaining the unbiased heterogeneity that was revealed by cluster analysis. In particular, tumor cells, TAMs, astrocytes, and vasculature presented with consistent heterogeneity in both the PDGfp and BrM datasets. Four tumor cell clusters were identified in PDGfp samples, Tumor_A-C and the category OPC-Like which could not reliably be distinguished from normal oligodendrocyte precursor cells (OPCs) due to their known phenotypic similarity. These subtypes were stra- tified based on their expression of GFP, Olig2, Sox2, and Sox9 (Sup- plementary Fig. 5c). We further analyzed the proportion of Ki67 positivity, showing that tumor clusters A, B, and C had equivalent proportions of proliferating cells in untreated and regrowth tumors (Supplementary Fig. 5d). OPC-Like cells showed the least proliferation, but still sufficient levels to indicate that this population included neoplastic cells, despite their lack of the GFP tumor cell marker (Supplementary Fig. 6a). This is consistent with the reported ability of glioblastomas to recruit non-transformed cells into the TME36. Inter- estingly, the Tumor_A phenotype, which had the highest expression of GFP, was the most impacted by IR treatment in terms of proliferation, showing almost no proliferation at 7-days post treatment (Supple- mentary Fig. 5c, d). This phenotype is consistent with quiescent resi- dual glioma cells that are both radio- and chemo-resistant37. In BrM, we also identified four tumor cell clusters, termed Tumor_A-D, stratified based on expression of EpCAM, cytokeratin 8, cytokeratin 14, and E-cadherin (Supplementary Fig. 5e). Each BrM tumor cluster showed almost complete loss of proliferation following IR (Supplementary Fig. 5f). The percent of proliferating cells in regrowth BrM did not match that of untreated tumors, but this likely does not reflect long-term alterations in proliferative capacity caused by treatment. For animal welfare reasons, and tomaintain comparable tumor sizes, regrowth BrM were harvested when progression was measurable by MRI, not at the point that their growth curve matched that of untreated tumors. Tumor clusters in BrM showed consistent spatial localization across treatment types, with the type A phenotype comprising themajority of the tumor core, while type B was dominant in areas of dense luminal structure, and C and D was predominantly localized to the tumor border (Supplementary Fig. 6b). TAM populations were delineated in both the PDGfp and BrM datasets based on their expression of CD68, P2ry12, Iba1, CD206, and CD45 (Supplementary Fig. 5g). TAM_A cells showed higher overall expression of CD68 and CD45, resembling myeloid-derived macro- phages (MDMs). TAM_B showed higher expression of Iba1, which may represent both activated MDMs and resident microglia. TAM_C showed the highest expression of P2ry12, which along with their localization outside of the tumor mass, suggests they are pre- dominantly brain-resident microglia (Supplementary Fig. 6c). TAM_D cells showed the highest expression of CD206, localizing primarily to meningeal regions in similar distribution patterns to meningeal mac- rophages. The remaining immune cells that could not be reliably classified as known cell types were classified generally as CD45+, and represent a range of infiltrating myeloid cell types. Astrocytes also clustered into two distinct patterns, with Astro- cyte_A cells showing lower expression of GFAP and localizing to distant brain areas outside of the tumor, while Astrocyte_B cells showed higher expression of GFAP, and localized to the border regions of both PDGfp and BrM tumors (Supplementary Figs. 5h and 6d). Neurons stratified into two subtypes in PDGfp tumors based largely on their expression of FOXP1, which was specific to Neuron_B cells. FOXP1 was not included in theBrMantibodypanel, and soNeuron_B cellswere not identified in this dataset. Interestingly, blood vessels also separated intomultiple categories in both PDGfp and BrM tumors, despite a lack Article https://doi.org/10.1038/s41467-024-47185-9 Nature Communications | (2024) 15:3226 8 of diverse phenotypic markers for the endothelium in the antibody panelswe haddeveloped (Supplementary Fig. 5i). Vessel_A cellsmostly localized to brain regions outside of the tumor, while Vessel_B cells were associated with dysregulated vasculature within both PDGFp and BrM tumors. Vessel_C type cells were further stratified in PDGfp tumors based on high expression of αSMA (Supplementary Fig. 6c). Region and cell type quantification Quantification of PDGfp tumors across treatment types revealed pro- nounced changes in the cellular landscape of glioblastoma samples following IR therapy (Fig. 5a). Untreated gliomas were predominantly comprised of tumor cell types, neurons, and astrocytes, with a mixed population of TAMcell types, and very small populations of T cells and neutrophils. PDGfp tumors 7-days post-IR showed the expected depletion of tumor cell populations, and a notable increase in TAM_A, TAM_B, T cell, neutrophil, and fibroblast populations. TAM popula- tions more than tripled in percent-total, with marked increases in the TAM_A population in particular. Regrowth tumors largely recapitu- lated the cellular landscape of untreated tumors, with alterations in tumor type percentages, but with an increase in T cells compared to untreated tumors. Conversely, BrM tumors showed no widespread alterations in cellular landscape following treatment, or upon tumor regrowth, other than an increase in T cell populations (Fig. 5b). Intratumoral regions were also significantly altered by IR in PDGfp tumors. There was a substantial increase in the percent-total of ECM- rich fibrotic regions following treatment, which was subsequently reduced upon tumor regrowth (Fig. 5c). This correlates with a similar increase in the percentage of fibroblasts 7 days post-IR. Annotated structural regions in BrM tumors were much more consistent follow- ing IR treatment, with significant differences in only the ratio of tumor area to brain, consistent with the observations of reduced tumor volume following treatment, and subsequent increases at the point of regrowth (Fig. 5d). We further assessed the cellular composition of each of the annotated structural regions in both PDGfp and BrM tumors at each treatment point (Supplementary Fig. 7a–e). There were significant increases in TAMs, T cells, and fibroblasts in Border, ECM, and Peri- vascular regions of PDGfp tumors, which correlates with similar increases in areas of fibrosis (Fig. 5e). Conversely, BrM tumors only showed significant alterations in the distribution of tumor types and T cells in response to treatment (Supplementary Fig. 7d, e). Graph-based spatial network analysis To analyze spatial relationships of cell types in HIFI data we employed orthogonal graph network and clustering-based approaches. We used graph-based network analysis to directly assess and summarize con- sistent cellular distance relationships from combined spatial data for each treatment type (Fig. 5f, g). Proximity network graphs used the Fruchterman-Reingold algorithm to set edge lengths betweennodes as the mean distance between that cell type and the nearest neighbor of an alternate cell type. Edge weight was set as the inverse of standard deviation to reveal which interactions are more highly conserved between images. Edge weights that fell below a preset threshold were removed from the network. Node size for each cell type shows the relative percentage of those cells in the total population, binned into discrete node sizes. Network graphs were then clustered into sub- domains based on similar nodal connections. Network plots for PDGfp tumors revealed consistent cellular relationships within and between treatment types, along with pro- nounced changes between tumor and immune cell populations fol- lowing IR. Specifically, network plots revealed an increased association between surviving Tumor_A cells after treatment and increased populations of TAM_A, Fibroblasts, and Astrocyte_B cells. This corre- lates with the observed increase in fibrosis, and indicates a potential survival niche for radioresistant tumor cells (Fig. 5c, e). Additionally, theseplots show the global extent of organized patterning, with 7 days post-IR samples having far less architectural structure and compart- mentalization compared to untreated tumors. This was further reflected by cellular heterogeneity of each network; 2.35 mean Shan- non Diversity Index for 7 days post-IR, versus 2.01 for both untreated and regrowth samples (Fig. 5f). Consistent with the results discussed above, BrM network orga- nization did not change substantially following treatment, with similar cellular distance relationships, and clear cellular compartmentaliza- tion depicting the tumor, surrounding brain, and tumor-brain border (Fig. 5g). Cellular diversity was alsonot substantially different based on mean Shannon Diversity. Cell neighborhood analysis Cellular neighborhood analysis38 was performed with imcRtools39 by packaging HIFI data as SingleCellExperiment objects in R40. Treatment conditionswere independentlypooled for PDGfpandBrMdatasets, and 15 cellular neighborhoods (CNs) were identified in each (Fig. 6a, b). Metadata for each cell was annotatedwith theCN they belonged to, and the percent-total of all CNs was quantified for each image in each treatment condition (Fig. 6c, d). The total proportion for 13 of 15 CNs was found to change significantly following treatment in PDGfp tumors, while only 4 CNs were significantly different between untreated and 7 days post-IR in BrM tumors. This again demonstrates the absence of spatial reorganization of BrM lesions in response to IR. Of particular interest were the CN populations in PDGfp samples that were essentially unique to post-IR treated samples; CN2, CN5, and CN14. These 3 CNs clustered together in terms of cell type composi- tion, being predominantly comprised of T cells, Fibroblasts, Astro- cyte_B, Neutrophils, TAM_A, TAM_B, Vessel_B, Vessel_C, and small percentages of Tumor_A cells (Fig. 6a). Each of these CNs were mostly specific to regions of regressed lesions in samples 7 days post-IR treatment (Fig. 6e). Visual validation of regions enriched in these 3CNs (Fig. 6f–h) correlate with fibrotic regions identified by machine- learning annotation (Fig. 5c). Cell interaction analysis was performed to assess significant colocalization of cells within a 30 µmdiameter in untreated and 7 days post-IR PDGfp tumors39,41. The Tumor_A cell type was found to be significantly anti-correlated with Fibroblast, T cell, Neutrophil, and Astrocyte_B cell types prior to treatment. 7 days post-IR treatment, Tumor_A cells became significantly colocalized with each of these cell types, as well as TAM_A, and CD45+ cells. (Fig. 6i, j). These results corroborate the previous orthogonal spatial analyses: each of the cell types cluster together in proximity network analysis of 7 days post-IR samples (Fig. 5f), they comprise CNs 2, 5, and 14 (Fig. 6a), and each are observed to be increased in the fibrotic ECM niche of treated samples (Fig. 5e, Supplementary Fig. 7c). The presence of the non-proliferative Tumor_A phenotype in this fibrotic niche suggests these spatial superstructures represent a survival niche for dormant radioresistant tumor cells (Watson, Zomer, [..], and Joyce, manuscript in revision). The proximity network and cell neighborhood analyses of BrM tumors indicated that they do not respond to IR treatment in terms of spatial reorganization as PDGfp tumors do. Rather, the survival mechanism for IR-treated BrM samples appears to primarily rely on tumor cells entering into a quiescent state. Proliferating tumor cell populations are significantly reduced for all identified tumor types at 7 days post-IR (Supplementary Fig. 5f). Based onMRI volume data and image analysis, we did not observe significant reductions in tumor volume or cleaved-caspase 3+ apoptotic cells at the 7-day timepoint. These combined data indicate that BrM tumor resistance to IR therapy is driven predominantly by enrichment for lower proliferating cells, or by cell-state switching to quiescent states, rather than the formationof spatially protective niches. However, despite the difference in survival mechanisms, single dose focalized IR treatment was equally and transiently effective in BrM tumors as for PDGfp tumors. Article https://doi.org/10.1038/s41467-024-47185-9 Nature Communications | (2024) 15:3226 9 Discussion High-dimensional digital pathology is a powerful approach for deriving biologically meaningful data from complex tissues, such as the multicellular TME. Our goal herein was to create a workflow that was fully agnostic to tissue processing, and required no special reagents, conjugated antibodies, or expensive equipment, so as to lower the barrier of entry to highly-multiplexed image analysis. By optimizing the workflow for a challenging ‘worst-case scenario’, in this case, lightly-fixed cryo-embedded glioblastoma and brain metastasis tissues, we have developed a robust imaging and analysis pipeline that is non-destructive for the tissue, allows for whole-slide image sizes, uses standard off-the-shelf antibodies, and that is scal- able to high-throughput. This workflow includes the development of software tools that can handle large extracellular structural a c e Border FibrosisPerivascular f Untreated 7 Days Post-IR Regrowth 0% 25% 50% 75% 100% Untreated 7 Days Post-IR Regrowth Cell Types Tumor_A Tumor_B Tumor_C Neuron_A Astrocyte_A Astrocyte_B CD45 Fibroblast Neutrophil T_Cell Vessel_A Vessel_B TAM_A TAM_B TAM_C TAM_D Cell Types Tumor_A Tumor_B Tumor_C OPC_Like Neuron_A Neuron_B Astrocyte_A Astrocyte_B CD45 Fibroblast Neutrophil T_Cell Vessel_A Vessel_B Vessel_C TAM_A TAM_B TAM_C TAM_D Percent Total 0% 25% 50% 75% 100% Untreated 7 Days Post-IR Regrowth Tumor Brain Border ECM Vessel Perivascular 0% 20% 40% 60% 80% 100% Pe rc en tT ot al Ar ea Untreated 7 Days Post-IR Regrowth Tumor Brain Border ECM Vessel Perivascular 0% 20% 40% 60% 80% 100% Pe rc en tT ot al Ar ea Untreated 7 Days Post-IR Regrowth b d Untreated 7 Days Post-IR Regrowth Astr ocyt e_B Neuron_A TA M_A Tumor_A Tumor_C 0% 20% 40% 60% Pe rc en tT ot al C el ls in R eg io n Astr ocyt e_B Fibroblast TA M_A T_Cell Tumor_A 0% 20% 40% 60% Pe rc en tT ot al C el ls in R eg io n Fibroblast TA M_A T_Cell Vesse l_B 0% 10% 20% 30% 40% 50% Pe rc en tT ot al C el ls in R eg io n Tumor_A Tumor_B Tumor_C OPC-Like Neuron_A Neuron_B Astrocyte_A Astrocyte_B CD45+ Fibroblast Neutrophil T_Cell Vessel_A Vessel_B Vessel_C TAM_A TAM_B TAM_C TAM_D µH’ = 2.35 µH’ = 2.01µH’ = 2.01 g Tumor_A Tumor_B Tumor_C Tumor_D Neuron_A Astrocyte_A Astrocyte_B CD45+ Fibroblast Neutrophil T_Cell Vessel_A Vessel_B TAM_A TAM_B TAM_C TAM_D µH’ = 1.78 µH’ = 1.77µH’ = 1.78 Article https://doi.org/10.1038/s41467-024-47185-9 Nature Communications | (2024) 15:3226 10 features, endogenous tissue autofluorescence, and densely packed irregular nuclei. Some recent cyclic-IF approaches can now achieve between 60–100 markers42, which is especially useful for rare and precious patient samples where researchers seek to extract the maximum data possible. HIFI is envisioned as a platform to add to researchers’options when dealing with samples and research questions that are not well addressed by other current techniques, and with the important components of ease of use and accessibility to a broad research community. Combining large imaging areas, high magnification, and multiple markers of ECM components herein enabled a much broader envir- onmental context for exploring cellular interactions in the TME. In addition to directly interrogating changes in fibrosis, desmoplasia, and healthy brain ECM, thesemarkers facilitated improved classification of tumor domains by machine learning classifiers. Automated classifiers required minimal training due to the depth of our image data and provided valuable data points for spatially aware cellular analysis. These integrated metrics facilitated multiple orthogonal computa- tional analyses to assess both the spatial organization of glioblastoma and brain metastasis models, as well as how that microenvironment organization responds to radiotherapy.Weobserved significant spatial reorganization in our murine glioblastoma model in response to IR, both in terms of cell landscape, spatial relationships, and super- structure patterning. Of particular interest were spatial niches that correlatedwith treatment-inducedfibrosis specifically in glioblastoma, harboring quiescent tumor types. These niches indicated a potential survival mechanism whereby the local environment promotes tumor dormancy and survival, leading to subsequent glioblastoma recurrence. This is consistent with recent reports demonstrating ECM- mediated tumor cell dormancy43, and the identification of gene sig- natures related to fibrosis being predictive of more rapid tumor relapse and reduced overall survival in glioblastoma patients44. Conversely, the breast-to-brainmetastasis model showed no such organizational response to IR treatment. Instead, these tumors demonstrated reduced proliferation following treatment, with no significant debulking, or changes in cellular landscape other than increased T cells after treatment. Nonetheless, both glioblastoma and brain metastasis preclinical models received similar transient survival benefits from radiotherapy, indicating markedly distinct survival mechanisms between different tumor types in the same host tissue. In summary, our HIFI approach is non-destructive, allows for large-scale imaging, and is adaptable for high-throughput analysis. Critically, the workflow reported herein is low-cost and opensource, such that it can be easily adoptable and broadly accessible for the scientific community. While we have demonstrated the use of the HIFI strategy in cancer tissues, this pipeline can equally be used in the analysis of any tissue type. We have purposely designed the workflow to be independent of the sample preparation, making it versatile for various types of tissue samples and eliminating the need for specia- lized reagents or equipment. Additionally, we have incorporated software tools to address the challenges of analyzing the complexity of the TME, such as identifying large extracellular structures and irregular nuclei. Our approach thus offers researchers a new strategy for studying the TME during tumor evolution, and following therapeutic intervention, enabling a deep and comprehensive interrogation of multicellular regional interactions. Methods Cell lines Cell lines were cultured in DMEM + Glutamax (Gibco) containing 10% fetal bovine serum (FBS, Gibco) and 1% penicillin/ streptomycin (Gibco). DF1 chicken fibroblast and PyMT-BrM3 were grown under adherent conditions. The PyMT-BrM3 breast cell line was derived from the murine parental 99LN cell line, which was isolated from a meta- static lymph node lesion that arose in the MMTV-PyMT (murine mammary tumor virus; polyoma middle T antigen) breast cancer model (C57BL/6J background). This cell line was sequentially selected three times in vivo for brain-homing capacity, resulting in the PyMT- BrM3 variant used herein5. All cell lines were routinely authenticated for morphology and growth dynamics, and tested for mycoplasma contamination. Animal models, treatments, tissue processing All animal studies were approved by the Institutional Animal Care and Use Committees of the University of Lausanne and Canton Vaud, Switzerland (License numbers: VD3804 and VD3688). Mice were housed in theAgora In VivoCenter (AIVC) animal facility in individually ventilated cages, under a 12 h light/dark schedule at 22 °C and in the presence of 2–4 cage mates. Standard autoclaved lab diet and water were provided ad libitum. Murine genetically-engineered mouse models (GEMMs) of glio- blastoma were generated as previously reported27,45–47. Nestin-Tv- a;Ink4a/Arf−/− mice in the C57BL6 background were bred and main- tained at the Agora Cancer Research Center, University of Lausanne (UNIL), Switzerland. At 5–6 weeks of age, glioblastomas were induced in GEMMs by injection of DF-1 cells producing viral vectors containing both PDGF-B and GFP as described previously45, and monitored biweekly by MRI for tumor development. Female C57BL6 mice received intracranial injection of the PyMT-BrM3 breast-to-brain metastasis cell line29 at matched age and cranial coordinates for tumor initiation in the PDGfpmodel. Themaximum humane endpoint tumor size approved per out animal protocols was 1.5 cm3 for PDGfp and 0.2 cm3 for BrM, all mice were euthanized before reaching the maximum allowed tumor burden. Once tumors exceeded 20mm3, mice were randomly assigned without blinding to treatment groups (untreated, 7 days post-IR, or rebound). Mice were anesthetized by isofluorane and administered a single whole-brain focalized dose of ionizing radiation using the Precision X-Ray X-RAD SmART irradiator (10Gy for PDGfp mice, 15 Gy for BrM mice). MRI monitoring was continued following IR treatment, and mice were euthanized upon tumor recurrence as approved by the Institutional Animal Care and Use Committee. Animals were euthanized via pentobarbital injection, and per- fused intracardially with 10ml of PBS, followed by 10ml of PLP buffer. PLP buffer consisted of 1% paraformaldehyde (PFA), 0.2% NaIO4, 37.5% L-lysine and 37.5% P-buffer (containing 81% of Na2HPO4, 19% of NaH2PO4 diluted in water, pH = 7.4) (Sigma Aldrich). Brains were Fig. 5 | Spatial analysisof IR-treatedbrain tumors.Percentageof the total cellular composition quantified for each treatment group pooled across all images for (a) PDGfp and (b) BrM. Percentage of the area of tumor regions plotted for all groups for (c) PDGfp and (d) BrM. eCell composition of tumor border, perivascular niches, and fibrotic regions for selected cell types in each treatment group. Box-plots for c-e showpercent totals for each image (PDGfpUntreatedn = 19 images, 7 dayspost- IR n = 17 images, Regrowth n = 18 images, BrM n = 9 images for each treatment). p valueswere calculated using two-wayANOVA test. p values for c (from left to right): <0.0001, <0.0001, <0.0001, <0.0001, 0.0002, <0.0001, 0.0067, 0.0045. p values for d (from left to right): 0.0023, 0.0317, 0.0044. p values for e (from left to right): Border, 0.0002, 0.0256, 0.0002, 0.0411, <0.0001, <0.0001, <0.0001, <0.0001, <0.0001, <0.0001, Perivascular < 0.0001, 0.0003, <0.0001, <0.0001, <0.0001, <0.0001, 0.001, 0.0094, <0.0001, <0.0001, Fibrosis, 0.014, <0.0001, <0.0001, 0.0486, <0.0001, <0.0001. Cell adjacencygraphnetworkplots of pooled treatment types and associatedmean ShannonDiversity Index (μH’) for (f) PDGfp and (g) BrM microenvironments. Node sizes represent the binned range of percent-total cell populations, edge length is the mean distance between nearest neighbors, and edge width is the inverse standard deviation of mean distance. Gray outlines show clusters of nodeswith similar neighbor relationships. Source data are provided as a Source data file. Article https://doi.org/10.1038/s41467-024-47185-9 Nature Communications | (2024) 15:3226 11 excised and fixed in PLP overnight at 4 °C with gentle shaking. Tissue samples were washed in PBS, then transferred to 30% sucrose overnight at 4 °C with gentle shaking. Finally, mouse brains were flash- embedded in Optimal Cutting Temperature (OCT) compound (Tissue- Tek), then cryosectioned onto Fisherbrand Superfrost Plus slides at 10 µm thickness. All slides were stored at −80 °C until used for staining and imaging experiments. Antibody labeling and imaging Prior to study, all antibodies were validated for compatibility with the HIFI approach (Supplementary Note 1). Tissue section slides were thawed at room temperature (RT) and allowed to dry for 15min, then OCTwas removed in a PBS bath at RT for 5minwith gentle agitation. An optional step at this point is a post-fixation in 4% PFA on ice for 5min, followed 2 × 5-min PBS washes and quenching in 0.1M glycine for b c d a 12 6 7 13 10 11 4 8 1 15 3 9 5 2 14 Tum or_A Tum or_B Tum or_C O PC -Like N euron_A N euron_B Astrocyte_A Astrocyte_B C D 45+ Fibroblast N eutrophil T_C ell Vessel_A Vessel_B Vessel_C TAM _A TAM _B TAM _C TAM _D CN −3 −2 −1 0 1 2 3 3 1 11 12 9 15 6 8 10 14 5 7 13 2 4 Tum or_A Tum or_B Tum or_C Tum or_D N euron_A Astrocyte_A Astrocyte_B C D 45+ Fibroblast N eutrophil T_C ell Vessel_A Vessel_B TAM _A TAM _B TAM _C TAM _D CN 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 0% 20% 40% 60% Pe rc en tT ot al C el ls 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 0% 10% 20% 30% 40% Pe rc en tT ot al C el ls Untreated 7 Days Post-IR Regrowth e Untreated 7 Days Post-IRCN Cell Types 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 f g h CN2 CN5 CN14 DAPI-Col IV-GFP-CD31-GFAP-ER-TR7 Tu mor _A Astr oc yte _B TA M_A Fib ro bla st TA M_D CD45 + TA M_B All O the rs iHIFI Classified Cells T C ell s Astrocyte_A Astrocyte_B CD45+ Fibroblast Neuron_A Neuron_B Neutrophil OPC-Like T_Cell TAM_A TAM_B TAM_C TAM_D Tumor_A Tumor_B Tumor_C Vessel_A Vessel_B Vessel_C Astr oc yte _A Astr oc yte _B CD45 + Fibr ob las t Neu ron _A Neu ron _B Neu tro ph il OPC-Li ke T_C ell TAM_A TAM_B TAM_C TAM_D Tum or_ A Tum or_ B Tum or_ C Ves se l_A Ves se l_B Ves se l_C -10 0 10 sum_sigval Astrocyte_A Astrocyte_B CD45+ Fibroblast Neuron_A Neuron_B Neutrophil OPC-Like T_Cell TAM_A TAM_B TAM_C TAM_D Tumor_A Tumor_B Tumor_C Vessel_A Vessel_B Vessel_C Astr oc yte _A Astr oc yte _B CD45 + Fibr ob las t Neu ron _A Neu ron _B Neu tro ph il OPC-Li ke T_C ell TAM_A TAM_B TAM_C TAM_D Tum or_ A Tum or_ B Tum or_ C Ves se l_A Ves se l_B Ves se l_C -10 0 10 sum_sigval j Untreated 7 Days Post-IR Article https://doi.org/10.1038/s41467-024-47185-9 Nature Communications | (2024) 15:3226 12 20min at RT. Slides were then placed into humidified chambers, and tissue sections were outlined with a hydrophobic barrier using peroxidase-antiperoxidase (PAP) pens. Tissue sections were permeabi- lized with 0.2% Triton X-100 in PBS for 10min at RT, then washed twice in PBS at RT with gentle agitation. Slides were stained for nucleic acid with DAPI at 1:2000 dilution in HIFI Staining Buffer (HSB; 5% normal donkey serum (Merck) and 100mM NH4Cl (Sigma Aldrich) in PBS) for 10min at RT, then washed 3 times in PBS at RT with gentle agitation. SlowFadeDiamond (Invitrogen)mountingmediumand 22 ×22 cmglass coverslipswere used tomount slides for imaging. All slideswere imaged onaZeissAxio ScanZ1withColibri LED light source. TheLEDpowerwas set to the lowest possible settings based on the photonic energy of each wavelength. The 350 channel was set to 5% power, 488 set to 20%, 555 set to 50%, 647 set to 50%, and750 set to 100%. Exposure timeswere set to experimentally predetermined lengths to captureDAPI signal and endogenous tissue autofluorescence for all channels. Following ima- ging, slides were placed in a horizontal slide rack with the bottom removed, this rackwas then placed into a PBSbathwith gentle agitation until coverslips fell off without manipulation. Sections were blocked with HIFI Blocking Buffer (HBB; 10% normal donkey serum (Merck), 150mMMaleimide (SigmaAldrich), and 100mMNH4Cl (SigmaAldrich) in PBS) for 1 h at RT as previously described26. HBB was then removed and replaced with the primary antibody mix in HSB, and allowed to incubate for 1.5 h at RT on an orbital rocker. Slides were washed 3 times in PBS for 5min at RT with gentle agitation following primary antibody incubation, then secondary antibodymixwith conjugated fluorophores in HSB was added, and slides were incubated for 1 h at RT on an orbital rocker. All secondary antibodies were raised in donkey to optimize compatibility. Slides were washed 3 times in PBS for 5min at RT with gentle agitation following secondary antibody incubation, and directly conjugated antibodies in HSB were added, and slides were incubated for 1 h at RT on an orbital rocker. If no conjugated antibodies were included in that round of antibodies, slides proceeded to the next step. Slides received three final washes in PBS for 5min at RT with gentle agitation, andwere thenmounted for imaging. Following imagingof the first round of markers, coverslips were again removed, and antibodies were eluted by adding elution buffer (0.5M Glycine, 3M guanidine hydrochloride (Sigma Aldrich), 3M Urea (Sigma Aldrich), 40mM tris(2- carboxyethyl)phosphine (Sigma Aldrich), in deionized H2O) for 3min- utes with gentle agitation at RT. The above process was repeated for all rounds of antibody panels. HIFI staining was performed in two separate rounds of imaging,first for half of PDGfp samples (PDGfp1), then for the other half of PDGfp samples (PDGfp2) and all BrM samples. PDGfp1, PDGfp2, and BrM were each analyzed and clustered independently to avoid batch effects. PDGfp1 and PDGfp2 datasets were then merged at the point of cellular annotation. Image post-processing was performed for all 16-bit tile-scanned images using the Zeiss Zen software platform. Tile stitching and fusion was performed for all images. Background subtraction was performed with the rolling-ball subtraction method using a diameter of 75 µm. Image alignment, cell detection, and cell classification The alignment script was developed in Python to reassemble images from multiple imaging rounds into full-resolution OME TIFF images. For each HIFI image, the DAPI channel of the first imaging round was loaded to serve as the alignment reference. Subsequent DAPI channels were loaded in sequence. Rigid body (or affine) transformation aligned each DAPI channel to the reference by minimizing the squared dif- ference between the reference and the transformed channel using PyStackReg31. The transform was then applied to each image channel of that round, keeping only one channel in memory at a time. The transformed round was then saved to disk, and the next round was processed. Aligned channels from all rounds were concatenated into a single OME TIFF and compressed into a pyramidal image format. Metadata for pixel scale and proteinmarker identifiers of each channel were propagated to the final image. Cell segmentation was performed using the CNN-based StarDist algorithm33 implemented in Qupath. For improved accuracy, a model for deep learning segmentation was generated using manually seg- mented training data from the same image data set analyzed in this study. Images of 18 tissue samples were acquired which were derived from various healthy organs and tumors fromdifferentmousemodels. Image regions, with on average 305 cells per image, were manually selected to capture the wide range of variation in cellular appearance. The images were down-sampled from 2048 × 2048 to 512 × 512 for accommodating the receptive field of StarDist. Manual annotations were performedusingQuPath v0.3.x by eight volunteermicroscopists, who were assigned one image from each of the 18 tissue samples. The results in geojson annotation files were converted to label maps whereby overlaps were resolved by generating a distancemap for each cell and assigning pixels to a cell only when its value in the distance map was greater than for all other cells. This resulted in a manually annotated dataset of 144 images. This dataset was augmented by adding publicly available datasets including BBBC020, BBBC038v1, and BBBC039v1 from the Broad Bioimage Benchmark Collection48–50. For BBBC038v1 we used the fluorescence images of “stage1 train” only, from the unofficial fixes by Konstantin Lopuhin (https://github.com/ lopuhin/kaggle-dsbowl-2018-dataset-fixes). We also used the images from Coelho et al., which consists of hand-segmented nuclear images of 3T3 and U20S cells51. A probability score was generated for each nucleus predicted by StarDist to indicate nuclear segmentation confidence, and used for quality control checks. Each nucleus was expanded by 2.5 µm to approximate the surrounding cytoplasm. The expansion was con- strained by the size of the detected nuclei, so that a cell was not larger than 1.5 times the size of its nucleus. This produced fourmeasurement zones per object; nucleus only, cytoplasm only, whole cell, and cell membrane. The following measurements were taken for every single- cell object: X-Y object centroid coordinates of each nucleus, area, perimeter, and circularity of each nucleus and cell, distance of each object centroid to the nearest regional annotation border, and dis- tance of each object centroid to the nearest neighbor object of an alternate cell type. For subsequent analyses, either the nuclear MFI or the cell MFI were used, depending on the expression of the marker (Supplementary Table 2). Semi-supervised annotation was performed using FlowSOM, a tool that leverages self-organizing maps to interpret multiplexed flow cytometry data35. Self-organizing maps are artificial neural networks Fig. 6 | Differential Reorganization Response to IR Between Brain Tumor Models. Cellular composition column-scaled heatmaps of cell neighborhood (CN) analysis for (a) PDGfp and (b) BrM tumors. Percent total CNs across each treatment for (c) PDGfp and (d) BrM. Box-plots for c-d show percent totals for each image (PDGfp Untreated n = 19 images, 7 days post-IR n = 17 images, Regrowth n = 18 images, BrM n = 9 images for each treatment). p values were calculated using two- way ANOVA test. p values for c (from left to right): <0.0001, <0.0001, <0.0001, <0.0001, <0.0001, <0.0001, <0.0001, <0.0001, 0.0005, <0.0001, 0.0008, 0.0165, <0.0001, <0.0001, <0.0001, 0.0077, 0.0026, <0.0001, <0.0001, 0.008, <0.0001, <0.0001, 0.0116, 0.0003, <0.0001, <0.0001, <0.0001. p values for d (from left to right): 0.0013, 0.0014, 0.0003, 0.0008, 0.0122, 0.0471, 0.0006, 0.0079, 0.0003, 0.0027, 0.0004, 0.0052. e Representative positional plots for three untreated and three 7 days post-IR samples highlighting CNs 2, 5, and 14. f–h Representative images from7days post-IR samples showingHIFI (left) anddigital pathology (right) images of CN2, CN5, and CN14 respectively. Scale bars = 40 µm. Heatmaps of cel- lular colocalization calculated as the sum of two one-tailed permutation test p values (sum_sigval) for (i) untreated samples and (j) 7 days post-IR samples. Green boxes indicate significant colocalization and anticorrelation for Tumor_A cells, for example. Source data are provided as a Source data file. Article https://doi.org/10.1038/s41467-024-47185-9 Nature Communications | (2024) 15:3226 13 https://github.com/lopuhin/kaggle-dsbowl-2018-dataset-fixes https://github.com/lopuhin/kaggle-dsbowl-2018-dataset-fixes that perform unsupervised dimensionality reduction while maintain- ing the topological structure. We analyzed the two PDGfp batches and all BrM tumors separately, by building a self-organizing map with 10 × 10 grid dimensions, for a total of 100 nodes. Input data was the mean MFI intensity of marker proteins (nuclear or whole-cell accord- ing to thebiology) scaled to a0-1 range,with the highest values clipped to the 99.7th percentile. The assignment of nodes to clusters was gui- ded by: (i) the automatic optimizedmetaclustering by flowSOMon the minimal spanning tree, (ii) the heatmap of mean MFI intensities for marker proteins in each node, (iii) themapping of node annotations to QuPath objects and visual inspection of marker intensities. The annotation of clusters was defined across batches and mouse models based on the heatmap of the expression of marker proteins in clusters and visual inspection of cluster annotations inQuPath (Supplementary Fig. 4a–c). Machine learning object annotation Regional tissue annotation was performed in QuPath32 (version 0.4.4). Binary classifiers were generated to create ROIs for the entire lesion, vessels, areas of ECM, and tumor versus surrounding parenchyma. All classifiers used artificial neural network (ANN) machine learning, and were trained off of 25% of image data. Twoderivative annotations were created viamodifications of the above ROIs; the tumorborder, and the perivascular niche. The lesion classifier was trained using all available channels, holes and fragments less than 20,000 pixels were removed from ROIs, and the region was expanded 750 µm into the parenchyma to capture adjacent brain tissue. Manual removal of tissue deforma- tions resulting from tissue sectioning or misalignment was performed on all lesion ROIs. Vessel classifiers were trained using only CD31, VE- cadherin, and CD13 channels, holes and fragments less than 20 pixels were removed. Vessel ROIs were expanded by 15 µm to create peri- vascular annotations. The tumor versus brain classifier used for PDGfp samples was created using the DAPI, GFP, Sox9, Sox2, Olig2, NeuN, NF- H, CSPG5, WGA, and Lamin AC channels. The tumor versus brain classifier used for BrM samples was created using DAPI, E-cadherin, EpCAM, CK14, CK8, NeuN, and NF-H. Holes and fragments less than 4000 pixels were removed. The intersection between the tumor and brain was expanded by 70 µm to create the tumor border annotation. The ECM classifier was trained using the Col I, Col IV, CD13, WGA, CSPG5, Laminin, ERTR7, and Fibronectin channels. Holes less than 50 and fragments less than 300 pixels were removed from ROIs. Spatial analysis Cell neighborhood and interaction analysis were performed in R with RStudio using the imcRtools39 package version 1.0.2, and nearest neighbor distance network analysis used the igraph52 package. HIFI image measurements were repackaged into SingleCellExperiment objects for analysis. An expansion radius of 30 µm was used to gen- erate a spatial connectivity map of each image, connectivity maps for all images were pooled together and clustered to identify 15 distinct cell neighborhoods. Cell object identifiers were reassigned with their respective cluster number and quantified for each image and treat- ment type. Interaction analysis was performed using the ‘classical’ method, randomly reassigning each cell type with 1000 iterations to create a null distribution pattern for the statistical comparison of sig- nificant cellular interactions. Mean nearest neighbor distance measurements for each treat- ment typewere compiled into distancematrices ofmean and standard deviation, and distancenetworkswere generatedwith igraph using the Fruchterman–Reingold layout algorithm to set edge lengths as pro- portional to the mean distance between two cell types. Edge widths were modified based on the standard deviation of the mean, such that edge width was the inverse of variability (1/SD). This was then multi- plied by a constant esthetic value (100) for visualization of edge width heterogeneity. Edge widths below a threshold were removed from each network. The edge width threshold was determined as the med- ian standard deviation of a network multiplied by an experientially predetermined network heterogeneity cofactor constant (0.8) so that the network represents inherent heterogeneity while being compar- able to other networks in the experiment. Node size for each cell type was based on binned percent-total populations, so that cell type from <1% to >25% of the total population corresponded to node sizes between 3 and 70 respectively. Network plots were then clustered (gray outlines) using Louvain clustering with a resolution of 1.5 to determine cell type nodes with similar distance patterns. Statistics and reproducibility All relevant experiments and analyses had a minimum of 3 indepen- dent repetitions, and each dataset presented in the study had at least 3 technical replicates included. Statistical analyses and Shannon Diver- sity Index scores were performed in R with RStudio, or in GraphPad Prism (GraphPad9.0 software). Two-wayANOVA testswereused for all between-group comparisons, with the level of significance defined as P <0.05. Figure asterisks correlate to P-value thresholds: ns > 0.05, *≤0.05, **≤0.01, ***≤ 0.001, ****≤0.0001. Reporting summary Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article. Data availability The StarDist cell segmentation model and associated training data are available from the GitHub repository of the iMAXT consortium: https://github.com/TristanWhitmarsh/IMAXT_StarDist_Cellpose. All other data are available from the following Zenodo repository, DOI: 10.5281/zenodo.10778429. The file sizes of primary HIFI image data exceed hosting capacity, and so are available upon request. All remaining data can be found in the Article, Supplementary and Source data files. Source data are provided with this paper. Code availability The HIFI Alignment Tool is available from the GitHub repository of the Jean Hausser lab: https://github.com/jhausserlab/HiFiAlignmentTool. All other code is available from the following Zenodo repository, https://doi.org/10.5281/zenodo.10778429. References 1. Grech, N. et al. Rising Incidence of Glioblastoma Multiforme in a Well-Defined Population. Cureus 12, e8195 (2020). 2. Stupp, R. et al. Radiotherapy plus concomitant and adjuvant temozolomide for glioblastoma. N. Engl. J. Med. 352, 987–996 (2005). 3. Stupp, R. et al. Effects of radiotherapy with concomitant and adju- vant temozolomide versus radiotherapy alone on survival in glio- blastoma in a randomised phase III study: 5-year analysis of the EORTC-NCIC trial. Lancet Oncol. 10, 459–466 (2009). 4. Akkari, L. et al. Dynamic changes in glioma macrophage popula- tions after radiotherapy reveal CSF-1R inhibition as a strategy to overcome resistance. Sci. Transl. Med. 12 https://doi.org/10.1126/ scitranslmed.aaw7843 (2020). 5. Croci, D. et al. Multispectral fluorine-19 MRI enables longitudinal andnoninvasivemonitoringof tumor-associatedmacrophages.Sci. Transl. Med. 14, eabo2952 (2022). 6. Quail, D. F. et al. The tumor microenvironment underlies acquired resistance to CSF-1R inhibition in gliomas. Science 352, aad3018 (2016). 7. Zomer, A., Croci, D., Kowal, J., vanGurp, L. & Joyce, J. A. Multimodal imaging of the dynamic brain tumor microenvironment during glioblastoma progression and in response to treatment. iScience 25, 104570 (2022). Article https://doi.org/10.1038/s41467-024-47185-9 Nature Communications | (2024) 15:3226 14 https://github.com/TristanWhitmarsh/IMAXT_StarDist_Cellpose https://github.com/jhausserlab/HiFiAlignmentTool https://doi.org/10.5281/zenodo.10778429 https://doi.org/10.1126/scitranslmed.aaw7843 https://doi.org/10.1126/scitranslmed.aaw7843 8. Alexander, J. et al. Multimodal single-cell analysis reveals distinct radioresistant stem-like and progenitor cell populations in murine glioma. Glia 68, 2486–2502 (2020). 9. Berg, T. J. et al. The Irradiated Brain Microenvironment Supports Glioma Stemness and Survival via Astrocyte-Derived Transgluta- minase 2. Cancer Res. 81, 2101–2115 (2021). 10. Kievit, F. M. et al. Nanoparticle-mediated knockdown of DNA repair sensitizes cells to radiotherapy and extends survival in a genetic mouse model of glioblastoma. Nanomedicine 13, 2131–2139 (2017). 11. Ali, M. Y. et al. Radioresistance in Glioblastoma and the Develop- ment of Radiosensitizers. Cancers 12, 2511 (2020). 12. Gimple, R. C., Yang, K., Halbert,M. E., Agnihotri, S. & Rich, J. N. Brain cancer stem cells: resilience through adaptive plasticity and hier- archical heterogeneity. Nat. Rev. Cancer 22, 497–514 (2022). 13. Squatrito, M. & Holland, E. C. DNA damage response and growth factor signaling pathways in gliomagenesis and therapeutic resis- tance. Cancer Res. 71, 5945–5949 (2011). 14. Beach, C., MacLean, D., Majorova, D., Arnold, J. N. & Olcina, M. M. The effects of radiation therapy on the macrophage response in cancer. Front. Oncol. 12, 1020606 (2022). 15. McLaughlin, M. et al. Inflammatory microenvironment remodelling by tumour cells after radiotherapy. Nat. Rev. Cancer 20, 203–217 (2020). 16. Barker, H. E., Paget, J. T., Khan, A. A. & Harrington, K. J. The tumour microenvironment after radiotherapy: mechanisms of resistance and recurrence. Nat. Rev. Cancer 15, 409–425 (2015). 17. Giesen, C. et al. Highly multiplexed imaging of tumor tissues with subcellular resolution by mass cytometry. Nat. Methods 11, 417–422 (2014). 18. Lin, J. R. et al. Highly multiplexed immunofluorescence imaging of human tissues and tumors using t-CyCIF and conventional optical microscopes. Elife 7 https://doi.org/10.7554/eLife.31657 (2018). 19. Keren, L. et al. MIBI-TOF: A multiplexed imaging platform relates cellular phenotypesand tissue structure.Sci. Adv.5, eaax5851 (2019). 20. Norris, J. L.&Caprioli, R.M. Imagingmassspectrometry: anewtool for pathology in a molecular age. Proteom. Clin. Appl. 7, 733–738 (2013). 21. Kennedy-Darling, J. et al. Highly multiplexed tissue imaging using repeated oligonucleotide exchange reaction. Eur. J. Immunol. 51, 1262–1277 (2021). 22. Bai, Z., Su, G. & Fan, R. Single-cell Analysis Technologies for Immuno-oncology Research: from Mechanistic Delineation to Bio- marker Discovery.Genom. Proteom. Bioinforma. 19, 191–207 (2021). 23. Lewis, S. M. et al. Spatial omics andmultiplexed imaging to explore cancer biology. Nat. Methods 18, 997–1012 (2021). 24. Good, C. J. et al. High Spatial Resolution MALDI Imaging Mass Spec- trometry of Fresh-Frozen Bone. Anal. Chem. 94, 3165–3172 (2022). 25. Gendusa, R., Scalia, C. R., Buscone, S. & Cattoretti, G. Elution of High-affinity (>10-9 KD) Antibodies from Tissue Sections: Clues to theMolecularMechanismandUse in Sequential Immunostaining. J. Histochem. Cytochem. 62, 519–531 (2014). 26. Gut, G., Herrmann, M. D. & Pelkmans, L. Multiplexed protein maps link subcellular organization to cellular states. Science 361 https:// doi.org/10.1126/science.aar7042 (2018). 27. Hambardzumyan, D., Amankulor, N. M., Helmy, K. Y., Becher, O. J. & Holland, E. C. Modeling Adult Gliomas Using RCAS/t-va Technol- ogy. Transl. Oncol. 2, 89–95 (2009). 28. Ozawa, T. et al. Most human non-GCIMP glioblastoma subtypes evolve from a common proneural-like precursor glioma. Cancer Cell 26, 288–300 (2014). 29. Klemm, F. et al. CompensatoryCSF2-drivenmacrophage activation promotes adaptive resistance to CSF1R inhibition in breast-to-brain metastasis. Nat. Cancer 2, 1086–1101 (2021). 30. Bowman, R. L. et al.MacrophageOntogenyUnderlies Differences in Tumor-Specific Education in Brain Malignancies. Cell Rep. 17, 2445–2459 (2016). 31. Thevenaz, P., Ruttimann, U. E. & Unser, M. A pyramid approach to subpixel registration based on intensity. IEEE Trans. Image Process 7, 27–41 (1998). 32. Bankhead, P. et al. QuPath: Open source software for digital pathology image analysis. Sci. Rep. 7, 16878 (2017). 33. Schmidt, U., Weigert, M., Broaddus, C. & Myers, G. Cell Detection with Star-Convex Polygons; Lecture Notes in Computer Science (Including Subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics) (Springer International Publish- ing, 2018). 34. Ronneberger, O., Fischer, P. & Brox, T. U-Net: Convolutional Net- works for Biomedical Image Segmentation. Medical Image Com- puting and Computer-Assisted Intervention –MICCAI 2015. Lecture Notes in Computer Science, (Springer, Cham, 2015). 35. Van Gassen, S. et al. FlowSOM: Using self-organizing maps for visualization and interpretation of cytometry data. Cytom. A 87, 636–645 (2015). 36. Fomchenko, E. I. et al. Recruited cells can become transformed and overtake PDGF-induced murine gliomas in vivo during tumor pro- gression. PLoS One 6, e20605 (2011). 37. Bao, S. et al. Glioma stem cells promote radioresistance by pre- ferential activation of the DNA damage response. Nature 444, 756–760 (2006). 38. Schurch, C. M. et al. Coordinated Cellular Neighborhoods Orchestrate Antitumoral Immunity at the Colorectal Cancer Inva- sive Front. Cell 182, 1341–1359.e1319 (2020). 39. Windhager, J., Bodenmiller, B. & Eling, N. An end-to-end workflow for multiplexed image processing and analysis. bioRxiv https://doi. org/10.1101/2021.11.12.468357 (2021). 40. R Core Team. R: A language and environment for statistical com- puting. R Foundation for Statistical Computing. URL https://www.R- project.org/ (2021). 41. Schapiro, D. et al. histoCAT: analysis of cell phenotypes and inter- actions in multiplex image cytometry data. Nat. Methods 14, 873–876 (2017). 42. Kinkhabwala, A. et al. MACSima imaging cyclic staining (MICS) technology reveals combinatorial target pairs for CAR T cell treat- ment of solid tumors. Sci. Rep. 12, 1911 (2022). 43. Di Martino, J. S. et al. A tumor-derived type III collagen-rich ECM niche regulates tumor cell dormancy.Nat. Cancer3, 90–107 (2022). 44. Hoogstrate, Y. et al. Transcriptome analysis reveals tumor micro- environment changes in glioblastoma. Cancer Cell 41, 678–692 e677 (2023). 45. Pyonteck, S. M. et al. CSF-1R inhibition alters macrophage polariza- tion and blocks glioma progression. Nat. Med. 19, 1264–1272 (2013). 46. Dai, C. et al. PDGF autocrine stimulation dedifferentiates cultured astrocytes and inducesoligodendrogliomas andoligoastrocytomas from neural progenitors and astrocytes in vivo. Genes Dev. 15, 1913–1925 (2001). 47. Tchougounova, E. et al. Loss of Arf causes tumor progression of PDGFB-induced oligodendroglioma. Oncogene 26, 6289–6296 (2007). 48. Ljosa, V., Sokolnicki, K. L. & Carpenter, A. E. Annotated high- throughput microscopy image sets for validation. Nat. Methods 9, 637 (2012). 49. Caicedo, J. C. et al. Nucleus segmentation across imaging experiments: the 2018 Data Science Bowl. Nat. Methods 16, 1247–1253 (2019). 50. Caicedo, J. C. et al. Evaluation of Deep Learning Strategies for Nucleus Segmentation in Fluorescence Images. Cytom. A 95, 952–965 (2019). 51. Coelho, L. P., Shariff, A. & Murphy, R. F. Nuclear Segmentation in Microscope Cell Images: A Hand-Segmented Dataset and Com- parison of Algorithms. Proc. IEEE Int Symp. Biomed. Imaging 5193098, 518–521 (2009). Article https://doi.org/10.1038/s41467-024-47185-9 Nature Communications | (2024) 15:3226 15 https://doi.org/10.7554/eLife.31657 https://doi.org/10.1126/science.aar7042 https://doi.org/10.1126/science.aar7042 https://doi.org/10.1101/2021.11.12.468357 https://doi.org/10.1101/2021.11.12.468357 https://www.R-project.org/ https://www.R-project.org/ 52. Csardi, G. & Nepusz, T. The igraph software package for complex network research. InterJ. Complex Syst. 1695, 1–9 (2006). Acknowledgements Wewould like to thank Simona Avanthay for technical assistance during the initial development of experimental techniques, aswell as other past and present members of the Joyce lab for valuable scientific discus- sions. We would also like to acknowledge the important contribution of imaging experts who contributed to generating ground-truth micro- scopy training data: Iva Sutevski, Paola Guerrero Aruffo, Emilia Hurtado, Sabine Galland, and Michelle Ballabio. We also thank Philippe Thévenaz for his valuable advice. We are grateful to Véronique Noguet Brechbühl at the AGORA/FBM Plateforme technologique MPF for sectioning of pathological samples, and Sina Nassiri at the Swiss Institute for Bioin- formatics for valuable analytical advice. This study was greatly aided by scientific advice from members of the iMAXT Consortium, particularly Dario Bressan and JonasWindhager. Some figure graphics were created using BioRender.com. This researchwas funded by Cancer Research UK Cancer Grand Challenges – IMAXT Consortium (J.A.J., B.B.); the Swiss Cancer Research Foundation (KFS-5280-02-2021), Swiss National Sci- ence Foundation-Advanced Grant (TMAG-3_209224), Ludwig Institute for Cancer Research, University of Lausanne, and the Charlie Teo Foundation (to J.A.J.); The Brain Tumor Charity Future Leaders Fellow- ship, and The Brian Cross Memorial Trust (to S.S.W.); and the European Research Council (ERC) under the European Union’s Horizon 2020 Pro- gram under the ERC Grant Agreement no. 866074 (“Precision Motifs”) (to B.B.). JH acknowledges the support from the Swedish Cancer Fund (21 1731 Pj), the Swedish ResearchCouncil (2018-02530), SciLifeLab, and Karolinska Institutet. Author contributions Conception and study design, writing of the manuscript: S.S.W., J.A.J. Development of Methodology: S.S.W., B.D., A.dT., N.E., T.W., iMAXT Consortium, B.B., J.H. Acquisition and Analysis of Data: S.S.W., B.D., L.F., Z.K, M.M., J.H. Preparation of figures: S.S.W., B.D., Z.K. All others reviewedandapproved themanuscript. StudySupervision: S.S.W., J.A.J. Competing interests The authors declare no competing interests. Additional information Supplementary information The online version contains supplementary material available at https://doi.org/10.1038/s41467-024-47185-9. Correspondence and requests for materials should be addressed to Spencer S. Watson or Johanna A. Joyce. Peer review information Nature Communications thanks the anon- ymous reviewers for their contribution to the peer review of this work. A peer review file is available. Reprints and permissions information is available at http://www.nature.com/reprints Publisher’s note Springer Nature remains neutral with regard to jur- isdictional claims in published maps and institutional affiliations. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/ licenses/by/4.0/. © The Author(s) 2024 1Department of Oncology, University of Lausanne, Lausanne, Switzerland. 2Ludwig Institute for Cancer Research, University of Lausanne, Lausanne, Switzerland. 3Agora Cancer Research Center, Lausanne 1011, Switzerland. 4L. Lundin and Family Brain Tumor Research Center, Departments of Oncology and Clinical Neurosciences, Centre Hospitalier Universitaire Vaudois, Lausanne 1011, Switzerland. 5Department of Cellular and Molecular Biology, Karolinska Institutet and SciLifeLab, Stockholm, Sweden. 6Department of Quantitative Biomedicine, University of Zurich, Zurich, Switzerland. 7Institute for Molecular Health Sciences, ETH Zurich, Zurich, Switzerland. 8École Polytechnique Fédérale Lausanne, Lausanne, Switzerland. 9Machine Intelligence Laboratory, Department of Engineering, University of Cambridge, Cambridge, UK. 10Cancer Research UK, Cancer Grand Challenges iMAXT Consortium, University of Cambridge, Cambridge, UK. e-mail: drspencerwatson@gmail.com; johanna.joyce@unil.ch iMAXT Consortium Johanna A. Joyce10, Spencer S. Watson10, Tristan Whitmarsh10 & Bernd Bodenmiller10 Article https://doi.org/10.1038/s41467-024-47185-9 Nature Communications | (2024) 15:3226 16 https://biorender.com/ https://doi.org/10.1038/s41467-024-47185-9 http://www.nature.com/reprints http://creativecommons.org/licenses/by/4.0/ http://creativecommons.org/licenses/by/4.0/ mailto:drspencerwatson@gmail.com mailto:johanna.joyce@unil.ch Microenvironmental reorganization in brain tumors following radiotherapy and recurrence revealed by hyperplexed immunofluorescence imaging Results Overview of HIFI spatial analysis workflow Generation of IR treated brain tumor samples and antibody labeling�panel Whole-slide image alignment and registration Machine learning structural annotation and deep learning cell segmentation Clustering based semi-supervised cell classification Region and cell type quantification Graph-based spatial network analysis Cell neighborhood analysis Discussion Methods Cell�lines Animal models, treatments, tissue processing Antibody labeling and imaging Image alignment, cell detection, and cell classification Machine learning object annotation Spatial analysis Statistics and reproducibility Reporting summary Data availability Code availability References Acknowledgements Author contributions Competing interests Additional information