Introduction

Cancer is a complex disease driven by a sequence of genetic and microenvironmental changes. As a result, tumors comprise many different cell populations that are in constant evolution¹. This tumor heterogeneity has important biological and clinical implications as it differentially promotes key characteristics of individual cancer cells, including proliferation, invasion, and drug resistance ². Yet, we are far from fully understanding all the dynamic and bi-directional interactions in heterogeneous tumors, and how these interactions influence cancer cell behavior at the single-cell level. Emerging single-cell technologies and analysis tools provide a new opportunity to profile individual cells within tumors. For example, recent spatial transcriptomic and proteomic profiling provides evidence for the existence of multiple distinct microenvironmental niches within one tumor ³. However, these techniques can only provide static information and make use of thin tissue sections that are only one cell layer thick. To investigate the real-time behavior of individual cells intravital microscopy (IVM) has been developed to study tumor cell behavior evolution over time but also in their three-dimensional native tissue environment.

IVM is a technique used to capture both spatial and temporal information at a single-cell level in living animals. Over the past years, it has provided many new insights into the dynamics of tumor heterogeneity and the development of therapy resistance ^{4–14 6–10,12–15}. For example, intravital imaging of mouse breast cancer models revealed that only a small portion of tumor cells is motile ¹⁶, and that fast- and slow-migrating cells are present in distinct regions of the tumor ^7,17. Moreover, additional heterogeneity exists among migrating tumor cells that can, for instance, move randomly or directionally relative to other tumor cells and/or surrounding healthy tissue ⁷. However, traditional IVM analyses have considered cancer cells or the tumor microenvironment (TME) as a homogeneous population, thereby ignoring heterogeneity at the single-cell level and the influence of the local microenvironment, where interactions between neighboring cells may affect cellular behavior. To address these complexities, recent studies have used unbiased approaches using morpho-kinetic parameters derived from live in vitro¹⁸ or in vivo imaging to classify cells and explore the single-cell behavioral landscape as an additional omic layer ^19–21. In line with this concept of characterizing cellular dynamic properties for cell classification, we have previously developed an analytical platform termed BEHAV3D ^20,22 allowing to perform behavioral phenotyping of engineered T cells targeting cancer. This approach revealed substantial behavioral heterogeneity among the engineered T cells and identified, based on behavioral characteristics, a subpopulation of T cells with increased killing capacity.

However, the implementation of unbiased methods for cellular dynamic analysis, such as BEHAV3D, often demands computational expertise beyond the reach of many biomedical researchers conducting IVM experiments. To democratize the thorough understanding of heterogeneous cancer cell behavior and the microenvironmental components influencing it, there is a need for comprehensive, user-friendly, zero-code analysis tools. Such tools would facilitate unbiased analysis of diverse tumor dynamics, making advanced computational techniques accessible to a broader range of researchers.

Building upon our BEHAV3D pipeline for analyzing engineered T cell behavior²², we introduce BEHAV3D Tumor Profiler (BEHAV3D-TP) (Fig. 1), a new zero-code computational framework tailored to map the behavioral patterns of single tumor cells and their interplay with the TME. We used this novel tool to get more insights into the in vivo migration behavior of diffuse midline glioma (DMG), a highly aggressive pediatric brain tumor characterized by invasive growth and its capacity to spread to other regions of the brain and the spinal cord ²³. For the first time, by visualizing TME components either during or after IVM, we were able to show that DMG cell migration behavior correlates with specific microenvironmental factors, including tumor-associated macrophages and vasculature. In our study, we specifically focused on cancer cell migration as a specific feature of DMG, but our methodology can easily be extended to other cancer cell behaviors such as intravasation or metastatic colonization. Overall, BEHAV3D-TP provides an accessible computational approach for a better understanding of tumor-TME crosstalk and the functional implications of this communication in heterogeneous tumors.

Fig. 1:

Results

BEHAV3D Tumor Profiler: an accessible tool for mapping heterogeneity in tumor dynamics

To comprehensively map and analyze the dynamic heterogeneity of tumor cells within their TME we developed BEHAV3D-TP [https://github.com/imAIgene-Dream3D/BEHAV3D_Tumor_Profiler]. Aiming at enhancing accessibility and widespread adoption across scientists studying single tumor dynamics in live imaging, and inspired by previous work on democratizing advanced image analysis^24–26, we designed BEHAV3D-TP as a modular Jupyter Notebook featuring an intuitive graphical user interface (GUI) (Figure 1 and Supplementary Video 1). Operating directly within a web browser, the platform only requires a Google account for access. BEHAV3D-TP implementation within the Google Colaboratory (Colab) framework not only ensures efficient cloud-based computation with free and scalable resources but also enhances tool accessibility by eliminating the need for local software installation, a common barrier for inexperienced users. Regardless of programming proficiency, users can effortlessly navigate the computational pipeline to extract comprehensive insights into tumor dynamics (Figure 1). After the time-lapse image file processing and preparation (Figure 1a), BEHAV3D-TP includes a central module called the heterogeneity module (Figure 1b), enabling researchers to uncover intricate patterns of tumor behavior. Additionally, we provide two optional modules—the large-scale phenotyping module (Figure 1c) and the small-scale phenotyping module (Figure 1d)—which, combined with TME visualization during or after IVM, allow for correlating different tumor behavioral patterns with the composition of the TME. To demonstrate the application of BEHAV3D-TP, we utilized single-cell dynamic measurements obtained from intravital imaging of DMG-bearing mice.

Identification of highly heterogeneous in vivo cell behavior displayed by DMGs

To get insights into the behavioral heterogeneity among individual DMG cells, we performed IVM on mice implanted with a cranial imaging window to visualize cancer cells expressing H2BmNeonGreen marker for up to 5.4 hours (Figure supplement 1a). We tracked intravitally imaged DMG cells from different mice and in each imaged position we determined how each DMG cell’s position changes relative to the tumor edge (Figure 2a). Next, to untangle the complexity of tumor cell dynamics in an unbiased manner, we used the BEHAV3D-TP heterogeneity module, that implements multiparametric single-cell time-series classification allowing to identify distinct single-cell behavioral patterns (Figure 2 a, b). For this, we used kinetic parameters, including displacement, speed, persistence, and movement direction relative to the tumor edge, to measure similitudes between cellular behavior using the dynamic time warping algorithm²², followed by dimensionality reduction to classify cells. The BEHAV3D-TP heterogeneity module identified 7 DMG behavioral clusters: retreating, slow retreating, erratic, static, slow, slow invading, and invading (Figure 2b, c, d, Supplementary video 2). To accurately describe each of the clusters, we analyzed their most predominant features (Figure 2d) and back-projected the cluster information into the original time-lapse to visually inspect cell behavior in relation to the tumor (Figure 2a, c, Figure supplement 1b, Supplementary video 2). Of particular interest were cells from the invading and retreating clusters, characterized by fast movement away from and towards the tumor edge, respectively (Figure 2c-f, and Figure supplement 1c). In fact, in this particular DMG model, approximately 10% of the cancer cells were displaying invasive behavior (Figure supplement 1d, e), in line with the intrinsic infiltrative nature of DMG²⁷ . Interestingly, about another 10% of DMG cells were moving back to the tumor edge (retreating, slow retreating clusters) (Figure supplement 1d, e), behavior that has been previously observed for glioblastoma and suggested to be regulated by differential chemotactic signaling by the tumor edge and the surrounding brain parenchyma ⁷. We also observed erratic cells that display fast non-directed migration behavior, and static cells that are non-migratory (Figure 2c-f). Crucially, the distinct behavioral patterns were detected in all intravitally imaged mice (Figure supplement 1d, e), ruling out any sample bias that might lead to clusters consisting solely of cells from a single tumor or mouse. This underscores the effectiveness of the BEHAV3D-TP in revealing representative heterogeneous tumor behavioral patterns.

Fig. 2:

Mapping in vivo tumor cell migratory properties to distinct TME regions identified through ex vivo large-scale 3D imaging

To investigate whether the TME influences behavioral heterogeneity among cancer cells, we developed two complementary approaches to correlate the identified behavioral clusters (Figure 2) to the composition of the TME (Figures 1c, 3a and 4a).

Fig. 3:

Fig. 4:

The first approach consists of performing large-scale TME phenotyping and identifying regions with a specific cellular composition and architecture within the TME of intravitally imaged tumors that have been intravitally imaged and subsequently fixed. This is accomplished by using multispectral large-scale single-cell resolution 3D (mLSR-3D) imaging data ^28,29 and a spatial analysis framework called CytoMAP ³⁰ (Figure 3a). We selected microenvironmental components of interest, such as perivascular niches and glial cells, as they have been previously reported to play an important role in the invasive behavior of glioblastoma cells ^7,31–35. Specifically, intravitally imaged brains were collected immediately after the IVM imaging sessions and were fixed and stained for blood vessels (CD31), oligodendrocytes (Olig2), and tumor-associated microglia/macrophages (TAMM) (Iba1) (Figure 3b). Subsequent CytoMAP³⁰ spatial phenotyping analysis identified three different environmental niches: Void regions (less vascularized and with low glial infiltration), TAMM/oligo regions (enriched in oligodendrocytes, TAMM, and tumor cells), and TAMM/vascularized regions (highly vascularized regions enriched in TAMM) (Figure 3b, and Figure supplement 2a-c). To investigate the heterogeneity of tumor cell behavior in relation to the identified microenvironmental niches, BEHAV3D-TP large-scale phenotyping module mapped the different niches (Void, TAMM/Oligo, and TAMM/vascularized) identified by CytoMAP onto 3D intravital imaging data (Figures 1c, 3c). In the TME-defined IVM-imaged positions we compared the frequencies of the different behavioral clusters that we identified by with the BEHAV3D-TP heterogeneity module (Figure 3d). Interestingly, invading cells were more abundant in TAMM/vascularized regions compared to Void regions that contain a higher proportion of static cells, and TAMM/Oligo microenvironments contain more slow retreating cells compared to the other regions (Figure 3d). Finally, we compared the results obtained with BEHAV3D-TP to a more classical IVM analysis approach that relies on assessing single dynamic parameters. Interestingly, single parameters showed a restricted ability to identify differences in cellular behavior among various environments (Figure supplement 2e-g), underscoring the importance of multiparametric behavioral mapping in revealing more nuanced cellular behavior that cannot be captured solely by individual dynamic parameters.

Linkage of in vivo tumor cell behavior to the composition of the microenvironment at a single-cell level

With the goal of further refining TME phenotyping to better understand tumor cell behavior determinants, we implemented an alternative approach using the BEHAV3D-TP small-scale phenotyping module (Figure 1d, 4a). This module utilizes the detection of TME components during IVM to offer insights into their abundance and spatial relationship with DMG cells at a single-cell level. Together with the heterogeneity module that profiles the single-cell behavior of tumor cells (Figure 1), this additional module allows linking different DMG behavioral profiles with DMG cell small-scale niches (Figures 1d, 4a). For in vivo TME labeling, we injected mice with an anti-CD31-AF647 antibody to label blood vessels, and with either SR101, a fluorescent marker known to label oligodendrocytes and astrocytes ^36–39, or CD20r, a novel molecule described to label microglial cells ⁴⁰ (Figure supplement 3a, b). This allowed us to measure spatial TME characteristics of each single DMG cell, such as minimal distance to a certain TME component, the number of these cell types in a certain radius, or the cell density for this particular cell (Figure 1d). Interestingly, the BEHAV3D-TP small-scale phenotyping module revealed substantial TME heterogeneity across the different behavioral DMG clusters (Figure 4b). Firstly, compared to slow and slow-retreating cells, invading cells were usually found to be less densely positioned relative to other DMG cells (Figure 4b, and Figure supplement 3c), correlating with their invasive potential into the brain parenchyma. Second, and in line with the data from the large-scale TME phenotyping (Figure 3), our data suggests that invading DMG cells reside in a local environment that is enriched in TAMMs (Figure 4b). Indeed, invading cells are closer to TAMMs compared to erratic and slow migrating cells (Figure 4c), and move directionally towards TAMMs (Figure 4d, Supplementary video 3), suggesting that migration of this behavioral cluster is at least partially mediated by chemotaxis. Third, all fast-migrating behavioral clusters (invading, erratic, and retreating) were found to be generally closer to blood vessels (Figure 4b and e), suggesting that perivascular invasion is a more efficient movement route for these cells (Figure 4f), a mode of migration that has been previously observed for adult glioblastoma cells as well ^7,31,33–35. Lastly, particularly slow-retreating DMG cells were found in regions enriched for SR101⁺ cells (Figure 4b and Figure supplement 3d), perhaps reflecting oligodendrocyte-like (OC-like) tumor DMG cells that have been found near oligodendrocytes ⁴¹. Altogether, these data demonstrate the power of our BEHAV3D-TP pipeline in correlating single-cell tumor cell behavior to its local environment, with an important opportunity to identify novel environmental regulators of invasion.

Discussion

Here, we provide a new user-friendly platform tailored for IVM single tumor cell dynamics analysis, termed BEHAV3D-TP, aiming to better understand the mechanisms underlying tumor heterogeneity. In contrast to traditional pipelines that rely on the analysis of single parameters such as cell speed or displacement^42–44 or that use an artificial threshold to assign a specific behavior to cells¹⁰, BEHAV3D-TP is an unbiased approach based on the integration of multiple parameters measured by IVM. Using the heterogeneity module of our BEHAV3D-TP pipeline, we uncovered behavioral phenotypes that have not been described before and that could not have been identified using single parameter analysis. For example, our analysis revealed invading and retreating clusters, composed of cells with a similar cell speed but with an opposite migration direction, an important biological and clinical difference. Likely, bidirectional migration is mediated by different chemotactic cues and by the composition of the local TME⁴⁵. It is of vital importance for DMG and other cancer patients to better understand the biology underlying the different behavioral phenotypes, especially concerning invading cells that are difficult to therapeutically target ⁴⁶.

Previous studies suggest a major role for the microenvironment in promoting tumor cell spread into the healthy brain ^47–49. Indeed, using the complementary large- and small-scale TME phenotyping modules of BEHAV3D-TP, we identified substantial differences in the TME composition of the different behavioral DMG clusters. Our results also show that the small-scale TME analysis has a higher predictive value compared to the large-scale TME analysis, suggesting that a cell’s direct microenvironment is a more important regulator of its behavior compared to its macroenvironment. Future experiments directed towards identifying more TME components, for example using multispectral LSR-3D imaging²⁸, could therefore further refine the large-scale phenotyping module by detecting subregion heterogeneity.

IVM imaging allows for simultaneous labeling of various TME components and can be used to study live their association with tumor cell behavior. Indeed, we found that TAMMs and blood vessels were associated with specific behavioral phenotypes, showing the power of BEHAV3D-TP small-scale phenotyping module in identifying novel regulators of cellular behavior. This is exemplified by our finding that slow-retreating DMG cells were found in environments enriched for SR101⁺ cells, an observation that has not been made before. To take further advantage of BEHAV3D-TP, future work should be directed toward identifying more TME components, including cells of the adaptive immune system. For example, it would be interesting to use our pipeline on recently developed immunocompetent mouse models of DMG ^50,51, or in mice with a humanized immune system⁵².

To enhance the accessibility of BEHAV3D-TP for researchers keen on studying tumor cell behavioral profiles through live imaging, we have incorporated our software in a Google Colab notebook, allowing researchers with minimal programming skills to employ the tool online without the requirement to install specific software environments nor the need to directly interact with the code. The user can easily install the necessary dependencies, set up the working environment, process datasets and adjust parameters without requiring deep coding skills. Google Colab provides the necessary computational resources to be able to run the pipeline without demanding a powerful local machine. However, there are RAM usage and runtime limitations that could affect the analysis, given a big enough dataset in the free version of Google Colab. Given this and data privacy concerns, the workflow can also be run locally via a Jupyter Notebook, that can be downloaded from the Google Colab notebook. Aligned with the goals of BEHAV3D-TP, other methods utilizing the Google Colab framework are emerging, such as CellTracksColab²⁴, which is a general tool for tracking data analysis. While CellTracksColab offers valuable insights into tracking data, BEHAV3D-TP is specifically designed for IVM data and includes additional features. These features include directionality metrics crucial for invasion studies, detailed analysis of identified behavioral populations, data back-projection and correlation of these behaviors with both large-scale and small-scale TME features. This comprehensive approach provides a deeper understanding of the TME drivers behind intratumoral heterogeneity, making BEHAV3D-TP a robust tool for researchers investigating the complex behaviors of tumor cells in their natural environments.

In summary, BEHAV3D-TP allows to get insights into the heterogeneity of tumor cell migration behavior in an accessible way; however, it would be exciting to extend the application to other behaviors, such as therapy response to study TME-mediated drug resistance. Such an endeavor may benefit from including additional parameters beyond kinetic features. For example, others have additionally used morphological parameters ^19,21, and we envision that a multiparametric analysis based on a combination of morpho-kinetic features may further refine behavioral phenotyping.

Methods

Animal material

NOD.Cg-PrkdcscidIl2rgtm1Wjl/SzJ (NSG) mice were purchased from Charles River Laboratories (France). Experiments were conducted per Institutional Guidelines under acquired permission from the local Ethical Committee and as per current Dutch laws on Animal Experimentation. Mice were housed under sterile conditions using an individually ventilated cage (IVC) system and were fed with sterile food and water.

Human DMG sphere-forming culture

DMG007 (HSJD-DIPG-007) cells established from primary DMG material were grown in a base medium (TSM) consisting of 48% DMEM/F12 and 48% Neurobasal medium (Thermo Fisher) supplemented with GlutaMAX, 100mM Sodium Pyruvate, MEM non-essential amino acids and 1M Hepes buffer, at 1% each. Primocin (InvivoGen) at 50mg/ml was additionally added to the TSM. Working medium was prepared by adding 20ng/ml of human EGF and human basic FGF, 10ng/ml of PDGF-AA and PDGF-BB (Peprotech) and 2ug/ml of Heparin (StemCell Technologies). DMG cells were cultured at 37 °C in a humidified atmosphere with 5% CO₂.

Fluorescent reporter cloning and lentiviral transduction

Fluorescent reporter construct was based on the combination of Histone2B-mNeonGreen, P2A and mScarlet I-CAAX gifted from Dr. Hugo Snippert (UMC Utrecht, The Netherlands), and inserted by In-Fusion HD Cloning Plus (Takarabio) into a lentiviral vector including cPPT/CTS, WPRE, an EF1a promoter and IRES-puromycin-resistance cassette¹¹. Subsequently the above-mentioned elements were amplified by PCR and the amplicon was assembled in the final lentiviral backbone using Gibson Assembly (NEB). Sequences were validated using Sanger sequencing. DMG007 cells were transduced using a standard lentiviral transduction protocol and selected using puromycin to achieve a stable cell line.

Cranial imaging window (CIW) surgery and tumor cell injection

CIW surgery and tumor cell injection were performed on the same day as previously described ^7,10. Briefly, mice were sedated with Hypnorm (Fluanison [neuroleptic] + Fentanyl [opioid]) (0.4 ml/kg) + Midazolam [benzodiazepine sedative] (2 mg/kg) at a dose of 1:1:2 in sterile water and mounted in a stereotactic frame. The head was secured using a nose clamp and two custom-made 3D printed ear bars, made out of Polylactic Acid (PLA, Geeetech, 700-001-043). The head was shaved, and the skin was cut circularly. The mouse was placed under a stereo-microscope to ensure precise surgical manipulation. The periosteum was scraped, and a circular groove of 5 mm diameter was drilled over the right parietal bone. The bone flap was lifted under a drop of cortex buffer (125 mM NaCl, 5 mM KCl, 10 mM glucose, 10 mM HEPES buffer, 2 mM MgSO_4, and 2 mM CaCl₂, pH 7.4), and the dura mater was removed. Gelfoam sponge (Pfizer) was used to stop potential bleedings. Next, 1 × 10⁵ DMG007-H2B-mNeonGreen cells resuspended in 3 μl phosphate-buffered saline (PBS) were injected stereotactically using a 5 μl Hamilton syringe with a 2 pt style in the middle of the craniotomy at a depth of 0.5 mm. The exposed brain was sealed with silicone oil and a 6 mm coverslip was glued on top. Dental acrylic cement (Vertex) was applied on the skull surface, covering the edge of the coverslip, and a custom-made 3D printed ring made of biocompatible Polylactide acid (PLA, Geeetech, 700-001-0433) was glued around the coverslip to provide fixation during intravital imaging and to serve as a reservoir for a water drop required for the objective of the microscope. A single dose of 100 μg/kg of buprenorphine (Temgesic, BD pharmaceutical limited) was administered for post-surgical pain relief. Mice were closely monitored twice per week to assess behavior, reactivity, and appearance.

Intravital imaging

Mice were intravitally imaged as previously described ^7,10. In short, were sedated with isoflurane and placed face-down in a stereotactic frame. The Polylactide ring was used to fix the mouse’s head to the frame with custom-made Polylactide bars. Time-lapse images of the entire tumor volume were acquired for a maximum of 5.4 hours. The minimal time interval between serial images was set to 20 minutes. For tile scans, images of the complete z-stack of the tumor were acquired, with a step size of 2 µm. In a group of 3 mice, oligodendrocytes and astrocytes were imaged by intravenous injection of 50 µL SR101 (Thermofisher) at a concentration of 5 mM dissolved in PBS. In a different group of 3 mice, TAMMs were imaged by intravenous injection of 50 μL 100 µM CDr20 ⁴⁰. Simultaneously with the injection of either SR101 or CDr20, all mice received an intravenous injection of 10 ul CD31-AF647 antibody (ThermoFisher, A14716) to visualize the blood vessels.

Intravital imaging was performed on an upright FVMPE-RS multiphoton microscope (Olympus) equipped with a MaiTai DeepSee and an InSight DeepSee laser (Spectra-Physics), a 25x/1.05 numerical aperture water objective, and two GaAsP and two Multialkali photomultiplier tubes (Olympus). Images were acquired with simultaneous imaging with the MaiTai (960 nm) and Insight (1100 nm) lasers. Green, red, and far-red fluorescence was separated by 570 nm and 650 nm dichroic mirrors, and 610/70 nm and 705/90 nm filters (Semrock) and collected by GaAsP or Multialkali photomultiplier tubes. Scanning was performed in a bidirectional mode with a resonant scanner at 400 Hz and 12 bit, with a zoom of 1x, 512 × 512 pixels. For optimal detection of blood vessels (CD31-AF647), a single tile scan was made at the last timepoint with Insight laser tuned at 1250 nm. At the postprocessing step, this image was co-registered to the time-lapse movie.

Co-registration of IVM imaged blood vessel imaging

CD31-AF647 tile scan featuring optimal visualization of blood vessels during live imaging was co-registered to the last timepoint of the time-lapse movie. Both images (last timepoint of time-lapse tile-scan and blood vessel tile-scan) were visually inspected to determine and annotate the internal landmarks using Imaris (Oxford Instruments), versions 9.5 to 9.6. These landmarks served as a reference for the co-registration software to overlap both images. To perform the co-registration, a modification of elastix⁵³ was used; a python-based open-source software based on Insight Segmentation and Registration Toolkit (ITK) plus SimpleElastix⁵⁴ : a collection of different algorithms used on medical image registration. Due to the extensive number of potential parameter combinations, specific parameters were manually chosen. This selection was conducted using the Slicer3D visualization software, version 5.2.2.

Immunohistochemistry

3D immunohistochemistry was performed as described before ²⁹. After the live imaging session, mouse brains were immersed in 5 ml 4% paraformaldehyde (PFA) pH 7.4 overnight on ice. After fixation, brains were blocked in 5-10 ml Wash Buffer 1 (2 ml Tween-20, 2 ml Triton X-100, 2 ml 10% SDS, and 2 g BSA in 1 l PBS) for 5 hours. For multiplex immunolabeling, a two- or three-round staining protocol was used as previously described ^28,29. Washing and incubation steps were performed in Wash Buffer 2 (1 ml Triton X-100, 2 ml 10% SDS, and 2 g BSA in 1 l PBS). A thick (1 mm) slice of the cortex part that was under the cranial imaging window was dissected with a scalpel before proceeding to immunolabeling. The following primary antibodies were used: rabbit anti-Olig2 1:100 (AB9610, Merk), goat anti-Iba1 1:200 (NB100-1028SS, NovusBio) and CD31-AF647 1:250 (A14716, ThermoFisher). As secondary antibodies, donkey anti-rabbit AF405 1:250 (ab175651, Abcam) and donkey anti-goat AF633 1:250 (A21082, ThermoFisher) were used. After the last washing step, tissues were optically cleared by three stepwise incubations of an increasing concentration of FUnGI ⁵⁵ clearing agent (25%/50%/75%, diluted in PBS) for 1h at room temperature ²⁸. The final incubation step with 100% FUnGI was performed overnight at 4°C.

mLSR3D imaging

mLSR3D imaging ^28,29 of thick slices of fixed mouse cortex was performed using a Zeiss LSM880 system equipped with a 32-channel spectral detector. The signal from all fluorophores present in the samples was collected in a single acquisition and linear unmixing was used to obtain separate channels. Unmixing was carried out by the Online Fingerprinting mode of the microscope. Pre-acquired references of each single fluorophore were used to enable multispectral on-the-fly unmixing. Imaging was performed using a 25x multi-immersion objective 0.8 NA (Zeiss 420852-9871) with a working distance of 570µm. Tile scans were acquired using a voxel size of 0.332x0.332x1.200 µm, a pixel dwell of 2µs, and a 10% tile overlap.

Image processing

For 3D visualization, shift correction, object rendering, and cell tracking, time-lapse movies were processed with Imaris (Oxford Instruments), versions 9.5 to 10.1.

Cell tracking

The Channel Arithmetics Xtension was used for channel unmixing of overlapping signals. The Spots ImarisTrack module was used for object detection and semi-automated tracking of tumor cells (autoregressive motion). Tumor cell tracks were manually corrected to ensure accurate tracking. Around 50-100 cells per position were manually corrected and labeled as such for posterior selection of accurate tracks. To determine the direction of tumor cell movement relative to the tumor core, at the image edge closer to the tumor core a “Tumor edge” surface object was manually rendered using the contour tool. The Distance Transformation Xtension was used to measure the distance between tumor cells and the “Tumor edge”. For tracked tumor cells, time-lapse data containing the coordinates of each cell, the values of cell speed, mean square displacement, displacement delta length, displacement length, and distance to “Tumor edge” were exported. See https://qbi.uq.edu.au/files/40820/ImarisManual.pdf for statistical details.

Microenvironmental factors rendering

For rendering the microenvironmental factors detected during IVM or by post-IVM LSR-3D imaging, the Surface and Spots modules of Imaris were used. Surfaces were used for the SR101, CDr20, Olig2, DAPI, and Iba1 stainings, and Spots were used for CD31+ protruding or elongated structures, as previously described ³⁰. For rendered structures, positional data containing the coordinates of objects were exported.

Tumor large-scale spatial phenotyping with Cytomap

LSR-3D imaging data was processed for large-scale spatial phenotyping using the Cytomap Spatial Analysis Toolbox ³⁰. Briefly, neighborhoods with a 100 µm radius were generated. Neighborhoods were then clustered into 2 regions based on DMG cell positioning using the David Balwin algorithm. Using the gating tool, the tumor regions (containing DMG objects) were defined. Next, the neighborhoods were clustered again into 3 different regions based on the distribution of Iba1, Olig2, and CD31: Void, TAMM/Oligo, and TAMM/vascularized. By using the Manually Defined Regions option, only the tumor region was classified. Finally, for each region, fold-change values were calculated and exported using the Region Statistics tab. A 3D image of the cluster distribution was exported for each sample (Figure supplement 2a) and then used to visually assign phenotypical regions in the IVM movies (Figure 3c). For each mouse, one or two positions for each type of region were cropped and analyzed for tumor cell behavior as described above (Cell tracking).

BEHAV3D Tumor Profiler framework

BEHAV3D-TP was developed using the R Studio version 4.3.3. It features three modules: 1) Heterogeneity module with optional 2) Large-scale and/or 3) Small-scale phenotyping modules. All of the modules are included in the same script available in Github [https://github.com/imAIgene-Dream3D/BEHAV3D_Tumor_Profiler], but can be run independently from one another. To acquire this independency, the optional modules are directly combined to the necessary sections of the Heterogeneity module. Consequently, the user can use whichever module is best suited for their purposes. In order to further facilitate public access to the BEHAV3D-TP, a Google Colab notebook has also been created accessible through the same Github page [https://github.com/imAIgene-Dream3D/BEHAV3D_Tumor_Profiler]. This environment is friendly to anyone regardless of their coding expertise as it has a step-by-step follow-through execution through the pipeline. Nevertheless, a Wiki demo run and a video tutorial (Supplementary Video 1) are also provided through the GitHub page [https://github.com/imAIgene-Dream3D/BEHAV3D_Tumor_Profiler] for further clarification, if necessary. In this Colab notebook, the user can easily import data, modify analysis parameters, execute the pipeline and export the obtained results without any computational power demands or computational knowledge. Additionally, the original R version and R Markdown demo scripts are also provided.

Heterogeneity module

To account for missing time points and to create a time series with the same time interval for each tumor cell time series, linear interpolation was used to estimate the values of missing time points. To compare time series independently, tracks were cropped to the minimal common length among all tracks. In our case, to a 2.6 hour duration. For tracked tumor cells, time-lapse data containing the coordinates of each cell, the values of cell speed, mean square displacement, displacement delta length, displacement length, and relative distance to “Tumor edge” were used to calculate a cross-distance matrix between time series. First, we performed a principal component (PC) analysis and selected the number of PCs that explained at least 90% of data variance (3).

The previously calculated PCs were then used to compute the cross-distance matrix between each tumor cell multivariate time series using the dynamic time warping algorithm from the package “dtwclust”. To visualize heterogenous tumor cell behaviors in 2 dimensions, dimensionality reduction on the cross-distance matrix was performed by the Uniform Manifold Approximation and Projection method (“umap” package). Clustering was performed using the k-means clustering algorithm. For each cluster, the average values of distinct behavioral and environmental components were calculated and plotted as a heatmap (“pheatmap”).

Large-scale phenotyping module

The information obtained from mLSR3D and Cytomap about large-scale TME regions were combined with the DMG cells behavioral information to examine the large-scale features of tumor cells. For tumor cells, the mean speed, square displacement and raw movement were analyzed.

Small-scale phenotyping module

The positional coordinates of tracked tumor cells and detected microenvironmental factors were used to determine the small-scale landscape of tumor cells. For tumor cells, the average distance to the closest 10 neighbors was measured to determine tumor cell density. For each tumor cell and each environmental component, the number of microenvironmental objects in a 30-µm radius and the distance to the closest object were computed. For each behavioral cluster, the values of closest neighbors, minimal distance to each environmental factor, and number of environmental objects in the 30 µm radius were plotted and compared.

Statistical analyses

Statistical analyses were performed using R. Results are represented as mean ± s.e.m. unless indicated otherwise; n represents independent positions in 6 biological replicates. Two-tailed unpaired t-tests were performed between two groups unless indicated otherwise. To compare frequencies of different behavioral signatures between large scale regions identified, a Tukey Honestly Significant Difference (HSD) adjusted Estimated marginal means (emmeans) pairwise comparison test was applied. To compare the distribution of behavioral clusters between the environmental clusters, a mixed linear regression followed by a one-way ANOVA was performed. To compare the distribution of different behavioral and environmental features between the behavioral clusters, a one-way ANOVA followed by Tukey HSD correction was performed.

Data availability

Data will be made publicly available after publication.

Code availability

We provide the BEHAV3D-TP framework in a Google Colab user-friendly environment [https://colab.research.google.com/drive/1JI7ysqFf3tvdi6Df4YUsSZ8RbuXw8wba?usp=sharing], as well as in an R script format in Github [https://github.com/imAIgene-Dream3D/BEHAV3D_Tumor_Profiler]. Additionally, a Wiki [https://github.com/imAIgene-Dream3D/BEHAV3D_Tumor_Profiler/wiki] to follow through a demo dataset as well as a video tutorial (Supplementary Video 1) has been developed for further assistance.

Acknowledgements

We thank the Princess Máxima Center for Pediatric Oncology for technical support and the Hubrecht Institute and Zeiss for imaging support and collaboration. All imaging was performed at the Princess Máxima Imaging Center. We thank D. van Vuurden (Princess Máxima Center for Pediatric Oncology) and Dr. Ángel Montero Carcaboso (Sant Joan de Déu Barcelona Hospital) for providing primary DMG cultures; Dr Young-Tae Chang (Postech, Korea) for providing us the CD20r dye; Sandra Fernandez Archidona for her assistance with IVM data visualization; and members of the Dream3DLAB (Rios group) and imAIgene-lab (Alieva group) for offering critical feedback on the project. A.Z. was supported by a Veni fellowship from the Netherlands Organization for Scientific Research (09150161910076). R.C. and A.C.R were supported by an ERC-starting grant 2018 project (no. 804412). This work was produced with the support of a 2023 Leonardo Grant for Scientific Research and Cultural Creation, BBVA Foundation (LEO23-2-10305-BBM-BAS-144). M.A was supported by the Comunidad de Madrid (2022-T1/BMD-24021).