A prognostic human brain network for diffuse midline glioma

Patient cohorts

This study was approved by the institutional review board of GOSH NHS Foundation Trust (24/HRA/4335) with informed consent waived for retrospective analyses. Primary data collection was also approved by the institutional review boards of each participating study site before study commencement in accordance with the Declaration of Helsinki as amended⁶¹. Three independent cohorts of children diagnosed with primary pontine DMG/DIPG or thalamic DMG were identified and analysed retrospectively: (1) a discovery cohort from GOSH, UK; (2) an independent, multicentre external validation cohort from the Children’s Hospital Colorado (CHCO), USA; University of São Paulo, Brazil; Institute of Neurosurgery Dr. Alfonso Asenjo, Chile; and the HERBY clinical trial (NCT01390948)²⁴; (3) a second, independent cohort of children with biopsied pontine DMG, H3K27-altered, enrolled on PNOC clinical trials. Participating centres are major international institutions with recognized subspecialty expertise in paediatric neuro-oncology. PNOC trial inclusion criteria have been previously reported. Inclusion criteria were identical across cohorts 1 and 2, with children reported between January 2000 and January 2024 inclusive:

1. (1)

   Consensus clinical–radiological diagnosis of DIPG, defined as the rapid (<3 month) onset of cranial nerve palsies, long-tract signs or cerebellar signs accompanied by MRI identification of an expansile, diffusely infiltrative mass arising from and involving ≥50% of the pons and which is T1-hypo- or iso-intense, T2-hyperintense and lacks or has minimal contrast enhancement^20,21,22.

   or

2. (2)

   Neuropathological (tissue biopsy) diagnosis of a primary pontine or thalamic DMG, H3K27-altered, consistent with World Health Organization 2021 criteria¹.

   and

3. (3)

   Child (aged under 18 years).

4. (4)

   Treatment-naive brain MRI to include a three-dimensional volumetric T1-weighted sequence with high spatial resolution and two-dimensional T2-weighted or T2-weighted fluid-attenuated inversion recovery (FLAIR) sequences, both with a slice thickness of ≤5 mm.

5. (5)

   Complete clinical data available for analysis, defined as treatment covariates and time to last follow-up or death, as appropriate.

Neuropathological diagnosis was required for the inclusion of long-term survivors of DMG (defined a priori as an overall survival ≥18 months from diagnosis), irrespective of tumour primary pontine or thalamic location²³. Individuals with brain MRI of insufficient quality for analysis (for example, partial brain coverage or severe artefact(s) due to metal or motion) were excluded. Individuals lost to follow-up within 18 months (before the threshold for long-term survivorship) were excluded due to indeterminate outcome. Inclusion was based exclusively on the availability of the above data with no exclusions based on race, ethnicity, sex or other characteristics.

Clinical definitions

Age at diagnosis was defined as age at first brain MRI. Extent of resection was categorically stratified as unoperated, biopsy only, subtotal or gross-total by a board-certified paediatric neuroradiologist on postoperative brain MRI performed within 48 h of surgical resection. Gross-total resection was defined as the removal of all contrast-enhancing tumour volume while subtotal resection was defined as a greater than 10% but less than 90% reduction in tumour volume^62,63. Further data were also collected for treatment modality, including but not limited to adjuvant radiotherapy and chemotherapy. Overall survival was defined as the time from diagnosis to death. Surviving patients were censored at the date of database closure (20 September 2024). Long-term survivors who did not participate in follow-up were censored at the date of last clinical follow-up. The time to last clinical follow-up was defined as the time from diagnosis to last clinical encounter.

Neuroimaging acquisition, segmentation and preprocessing

All individuals in the discovery cohort underwent brain MRI on the Siemens MAGNETOM Avanto (1.5 T) or Prisma (3 T) scanner. Given the retrospective multicentre nature of the study, there was significant heterogeneity in terms of scanner manufacturer and magnet field strength across the external validation cohorts. Patient brain MRI data were downloaded from institutional picture archiving and communications systems as digital communications in medicine (DICOM) files. Tumours were manually segmented in three planes (axial, coronal and sagittal), with reference to all available MRI sequences acquired at diagnosis, by researchers with expertise in neuroimaging (J.S., V.L., F.S., B.S.P. and R.S.O.) using ITK-SNAP (v.4.2.0)⁶⁴. This yielded a binary tumour mask which was reviewed and, if necessary, corrected in consensus review with three board-certified paediatric neuroradiologists (K.M., A.B. and S.S.) consistent with a standardized, pan-institutional protocol for neuro-oncology MRI segmentation. After segmentation, the DICOM files were converted to Neuroimaging Informatics Technology Initiative (NIfTI) format using the nifti2dicom Python library (v.1.2.6). Brains were skull-stripped and extracted using SynthStrip and then registered to the Montreal Neurological Institute (MNI) MINC1 60MB paediatric template in MNI space (1.0 × 1.0 × 1.0 mm) using ANTs symmetric image normalization (SyN) registration^65,66,67. Cost-function masking was applied to tumour masks to remove them from computations⁶⁸. Tumours were then transformed to MNI space using forward SyN transforms and nearest-neighbour interpolation⁶⁹. All registered tumours were manually compared to original, patient-specific neuroimaging and reviewed for neuroanatomic accuracy in MNI space by a board-certified paediatric neuroradiologist (K.M.). All segmentations were performed blinded to patient identity and clinical outcome.

Human connectome data acquisition and pre-processing

Integrative analyses of fMRI and dMRI data were performed to clinically translate the previously reported neural integration of DMG through tumour-specific whole-brain connectivity. High-resolution normative structural (three-dimensional T1-weighted) and resting-state fMRI were obtained for 1,000 healthy children aged 9.0 ± 0.2 years from the Adolescent Brain Cognitive Development Study (ABCD1000)³⁵. High-resolution normative structural (three-dimensional T1-weighted) and advanced dMRI (multiple b values and directions) were obtained for 497 healthy children aged 13.0 ± 2.9 years from the Human Connectome Project (HCP) Lifespan Development dataset^36,37. Cohort demographics and imaging parameters for both normative connectomes have been previously reported^35,36,37. fMRI data preprocessing was performed using the Computational Brain Imaging Group preprocessing pipeline (https://github.com/bchcohenlab/BIDS_to_CBIG_fMRI_Preproc2016)^70,71. dMRI data were preprocessed in line with the HCP minimal processing pipelines (https://github.com/Washington-University/HCPpipelines.git)^72,73 and the resulting distortion-corrected b = 0 image was registered to MINC1 60MB through ANTs SyN registration^74,75.

Tumour location mapping

Identifying tumour voxels associated with patient overall survival

Tumour location frequency maps were plotted in MNI space as the sum of all binary tumour masks using ANTs iMath functions. To establish whether tumour anatomic location alone was associated with patient overall survival, we first implemented a univariate VLSM approach. A voxelwise general linear model with permutation testing was applied to the data to compare voxels between short- and long-term survival groups in NiiStat (https://github.com/neurolabusc/NiiStat)^29,76. Analyses were limited to voxels involved in ≥5% tumours, controlling for tumour volume and correcting for multiple comparisons^76,77,78,79. The significance at each voxel was assessed using 2,000 Freedman–Lane (shuffled) permutations. Given the conservative nature of voxel-based permutation testing⁸⁰, we replicated univariate VLSM by further implementing threshold-free cluster enhancement in FMRIB Software Library (FSL) permutation analysis of linear models (PALM)^77,81. Given the significantly shorter overall survival in children with pontine versus thalamic DMG, we also performed univariate VLSM across pontine DMG/DIPG and thalamic DMG separately to control for the prognostic effect of primary tumour location. In each of these analyses, we limited VLSM to voxels involved in at least seven patient tumours to avoid single case-driven effects, consistent with best-practice recommendations^76,78,79. We replicated analyses implementing multivariate support vector regression (SVR)-VLSM, which can identify complex dependencies between neighbouring voxels that traditional univariate VLSM approaches cannot^82,83. SVR-VLSM was implemented in SVR-LSM⁸⁴, a MATLAB toolbox for multivariate and cluster-based lesion symptom mapping, and limited to voxels involved in ≥5% tumours, controlling for tumour volume using direct total lesion volume control and correcting for multiple comparisons. Significance was assessed using 10,000 permutations. The results of univariate VLSM and multivariate SVR-VLSM were considered to be significant at a two-tailed FWE-corrected P < 0.05 and P < 0.005, respectively, as in previous work²⁹.

Interrogating the impact of tumour VLSM map scores on patient overall survival

The resulting univariate and multivariate VLSM maps represented unified, group-level maps of the prognostic strength of each voxel across overlapping tumours. We examined whether VLSM maps could determine patient overall survival based on tumour location alone, irrespective of the presence or absence of voxels significantly associated with overall survival in the discovery cohort, using data from an independent, multicentre external validation cohort, with identical inclusion criteria and identical methodology for tumour segmentation and brain MRI pre-processing. Each binary tumour segmentation in the external validation cohort was weighted by overlapping voxels in each VLSM map. A mean ‘tumour VLSM map score’ was then computed as the average of these weighted values for univariate and multivariate analyses, respectively²⁵. Patients were stratified into two risk groups using a median cut-off for each VLSM map: high risk (defined as median or above-median tumour VLSM map score) and low risk (defined as below-median tumour VLSM map score). Kaplan–Meier models were then used to model the overall survival for each VLSM map score-defined risk group, with the difference between groups evaluated using a two-tailed log-rank test. We also examined whether the tumour VLSM map score could act as an independent predictor of patient overall survival through inclusion in both univariable and multivariable Cox proportional hazards models incorporating clinical covariates for patient demography (age at diagnosis and sex) and clinical management (extent of resection, and completion of adjuvant radiotherapy and chemotherapy). In all instances, the proportional hazards assumption was verified before analysis and the results were expressed as HRs with 95% CIs and corresponding P values.

Tumour network mapping

Modelling tumour networks from tumour locations using normative fMRI

To define patient-specific tumour network maps, every tumour location was used as a ‘seed’ from which to compute resting-state functional connectivity with every other voxel in the brain across all individuals in the ABCD1000 normative fMRI dataset²⁶. The correlated time series (R-map) from seeded fMRI data per patient was then Fisher z-transformed and compared using a two-tailed, one-sample t test, thereby unifying each R-map into patient-specific voxelwise T-maps: a statistical representation of each tumour’s resting-state functional connectivity with all other brain voxels^26,85,86.

Identifying functional connections associated with overall survival

Resting-state functional connections associated with patient overall survival were identified by performing whole-brain general linear model permutation testing, correcting for multiple comparisons, and implemented in FSL PALM, adapting methods from lesion network mapping^77,81. First, a voxelwise one-sample t test was performed across short-term survivor patient T-maps (as previously defined, n = 106), correcting for tumour volume as a covariate with 2,000 Freedman–Lane sign-flip permutations⁷⁷. These analyses generated a unified, group-level statistical spatial map of significant voxels with greater and more consistent positive or negative resting-state functional connectivity across short-term survivor tumours—the DMG network. Next, to identify functional connections significantly different between short-term and long-term survivors, an independent voxelwise Aspin–Welch test was performed across all patient T-maps (n = 125), controlling for tumour volume as a covariate²⁹. Hereafter, we refer to this methodological framework as tumour network mapping. Regions of peak prognostic significance in the DMG network were identified using an unsupervised, cluster identification algorithm via FSL Cluster⁸⁷ using a minimum voxel count threshold ≥20⁷⁵. DMG network cluster peak significance values, centre of gravity coordinates and neuroanatomical location, hereafter termed network nodes or peaks, were reported⁸⁸.

Examination of DMG network robustness

To investigate the robustness of the DMG network to statistical choices, we compared voxelwise as opposed to threshold-free cluster enhancement and FWE as opposed to FDR approaches. We also examined the robustness of the DMG network to methodological choices and clinical covariates using the following experiments conducted in PALM: covarying patient overall survival with patient age at diagnosis and sex; restricting to pontine DMG/DIPG only; restricting to biopsied pontine DMG only; restricting to thalamic DMG only; and increasing the threshold for long-term survival from ≥18 to ≥24 months, consistent with both clinically reported thresholds^23,38. While all tumour segmentations normalized to MNI space were manually checked for neuroanatomic accuracy by a board-certified paediatric neuroradiologist (K.M.) against patient-specific imaging before seeding in ABCD1000, to further rule out the possibility that tumour-induced neuroanatomical distortions might influence tumour-specific normative functional connections, we also conducted an additional experiment covarying patient overall survival with a measure of mean peritumoral anatomic distortion. While also potentially reflective of patient head size, age and sex, the mean voxelwise SyN deformation fields required to normalize each patient brain to MNI space, excluding the tumour cost-function mask, were used as a measure of mean peritumoral anatomic distortion. Furthermore, given the brainstem-predominant location of DMG, we sought to exclude the possibility of the DMG network being derived from or confounded by fMRI signal originating within surrounding cerebrospinal fluid (CSF) (for example, secondary to heartbeat or CSF pulsations)⁸⁹. We verified this by implementing a voxelwise partial correlation analysis when computing whole-brain-to-tumour connectivity for each patient, controlling for temporal signal derived from a fourth ventricular CSF mask⁹⁰. Tumour network mapping was then repeated in PALM across these patient tumour functional connectivity maps with CSF signal partialled out⁵⁴. We also tested the convergence of our results across three normative connectomes with different preprocessing methods: ABCD1000³⁵, Yeo1000³⁹ and GSP1000⁴⁰. DMG network signal origin was also examined using increasingly strict, stepped statistical thresholds to visualize and eliminate noise. Finally, we also replicated analyses in PALM using a continuous, rather than binary, contrast for patient overall survival.

The correspondence in brain-wide spatial topography between the DMG network and each control network was assessed using both label permutation^32,91 and spin permutation^91,92 tests. First, the whole-brain spatial Spearman’s correlation, or ‘spatial similarity’, between each control network and the DMG network was compared against a null distribution of 10,000 permuted network results generated by randomly permuting short- and long-term survival labels 10,000 times across patients in the discovery cohort, yielding an empiric P value reflecting the significance of each observed spatial similarity. Next, spatial autocorrelation-preserving permutation testing, or ‘spin testing’, was achieved by parcellating the DMG network into paediatric MINCI 60 MB SynthSeg masks⁶⁵ and randomly rotating the spherical projections of the parcellated network 10,000 times. Each control network was parcellated into the same brain regions and its spatial similarity to each randomly rotated DMG network computed. The P value for each spin test was computed as for each label permutation test. Both spin and label permutation tests were performed and reported given the widespread use of spin testing despite reports that spin permutation tests probably inflate the statistical significance of spatial similarity between network maps, in contrast to label permutation-based approaches. Finally, the convergence of brain regions identified across all control analyses was also compared through the visualization of control network overlap maps plotted as overlapping binarized voxels surviving statistical thresholds.

Examination of network convergence through tumour network mapping on a normative structural connectome

Fibre orientations in preprocessed diffusion data were estimated using GPU-accelerated FSL Bayesian estimation of diffusion parameters obtained using sampling techniques (BEDPOSTX). Each tumour mask was transformed into participant-specific diffusion space and used as a ‘seed’ in GPU-accelerated probabilistic tracking with crossing fibres (ProbtrackX) in FSL with the following, default parameters: number of samples = 5,000; curvature threshold = 0.2; step length = 0.5 mm; subsidiary fibre volume fraction threshold = 0.01. b = 0 images were automatically segmented into CSF and grey matter with SynthSeg and were used as ‘exclusion’ and ‘waypoint’ masks, respectively in probabilistic tractography to obtain patient-specific whole-brain-to-tumour structural connectivity data⁶⁵. The streamline counts for each voxel were normalized by dividing each value by the total number of streamlines. Participant-level structural connectivity data were transformed to MINC1 60 MB and averaged across all individuals to generate population-level structural connectivity data for each tumour. Tumour-specific group-level structural connectivity maps were then analysed in PALM, in an identical manner to tumour network mapping described above for fMRI data, to examine the voxelwise association of tumour structural connectivity with overall survival²⁹. Convergence of control structural networks to the DMG structural network was tested through label and spin permutation analysis of spatial similarity, as described.

Reproducibility of the DMG network and the role of independent datasets

To ensure reproducibility of the DMG network and to examine its prognostic importance, we replicated these analyses across both independent, external validation cohorts, with identical methodology for tumour segmentation; brain MRI preprocessing; and tumour network mapping across functional and structural connectivity data. Convergence of both external networks to the DMG network was tested through label and spin permutation testing of spatial similarity, as previously described.

Evaluating the prognostic significance of the DMG network

Examining the impact of tumour functional connectivity on overall survival

We assessed the impact of mean network-to-tumour connectivity on patient overall survival. Mean network-to-tumour functional connectivity was derived by weighting the functional connections of each patient tumour to the DMG network with T values in the DMG network, restricted to only prognostically significant voxels, and then computing the mean of this weighted tumour connectivity. Hereafter, we refer to this as each patient’s DMG network-to-tumour functional connectivity. Such an approach has been used to test the out-of-sample robustness of connections associated with clinical improvements after effective deep brain stimulation for obsessive–compulsive disorder^25,66, Parkinson’s disease²⁵, Gilles de la Tourette’s disorder²⁵ and dystonia²⁵. Given systematic differences in absolute functional connectivity between pontine and thalamic tumours, for analyses in which both tumour locations were included, network-to-tumour connectivity values were normalized within each anatomic location to ensure that survival associations reflected relative network connectivity within a given anatomic context rather than gross, location-dependent differences in tumour connectivity. Patients were then stratified into two risk groups using a median cut-off: high risk (defined as median or above-median network-to-tumour connectivity) and low risk (defined as below-median network-to-tumour connectivity). Kaplan–Meier models were used to model overall survival for each connectivity-defined risk group, with the difference between groups evaluated using a two-tailed log-rank test. We also examined whether network-to-tumour connectivity acts as an independent predictor of patient overall survival through inclusion as a variable in both univariable and multivariable Cox proportional hazards models, as described for tumour location mapping analyses.

To confirm out-of-sample prediction, we fitted clinical-only and combined clinical-connectivity multivariable Cox proportional hazards models exclusively on the discovery cohort, incorporating patient age and sex and tumour location, DMG network-to-tumour connectivity, extent of resection, and adjuvant radiotherapy and chemotherapy. Model coefficients were subsequently frozen and the fully parameterized models were applied to generate overall survival estimates across the independent external validation cohort. Prognostic discrimination was quantified using the Harrell concordance index (C index) calculated from model-derived risk scores. The incremental predictive value of network-to-tumour connectivity was assessed as the difference in C index between combined and clinical-only models (ΔC index) across the external validation cohort. The ΔC index and the statistical uncertainty around the ΔC index were estimated using paired nonparametric bootstrap resampling of patients within the external validation cohort with 2,000 iterations. A 95% CI was calculated using the percentile bootstrap method. Two-tailed P values were derived from the bootstrap distribution of ΔC index relative to the null hypothesis ΔC index = 0.

To ensure that the prognostic impact of DMG network-to-tumour connectivity was robust to methodological choices, we repeated the above survival analysis using an orthogonal method to compute DMG network-to-tumour connectivity for each patient. DMG network peaks were used as weighted seeds in ABCD1000 to obtain a brain-wide distributed map representative of functional connectivity with brain regions of increased or decreased risk of shorter overall survival²⁹. Each patient tumour-specific functional connectivity map across the external validation cohort was spatially compared to this precomputed connectome discovery DMG network peak functional connectivity map with a whole-brain spatial Spearman’s correlation^26,29,31,67. The patients were stratified into high- and low-risk groups on the basis of a median cut-off in spatial similarity values, whereby high spatial similarities were indicative of high DMG network-to-tumour functional connectivity. We also repeated survival analyses using DMG network-to-tumour mean functional connectivity with maximum connectivity, minimum connectivity and the sum of connectivity values. Lastly, survival analyses were replicated across the DMG network as derived from Yeo1000 and GSP1000 adult connectomes, as well as tumour structural connectivity data obtained from the DMG structural network.

Testing the robustness of DMG network prognostic impact across DMG histone status and co-mutational burden

Given the potential for canonical and non-canonical oncohistone variants^41,93 and recent evidence for two prognostic, molecularly defined DMG subgroups containing co-mutations in either TP53 or BRAF/FGFR1 ^{(refs. 23,94}), we compared the mean network-to-tumour connectivity of each tumour in a subcohort of 87 children with paediatric solid tumour next-generation DNA panel sequencing using two-tailed Mann–Whitney U tests. The prognostic importance of pontine network-to-tumour connectivity was also confirmed across our external validation cohort restricting to children with biopsied pontine DMG only.

Specificity testing of DMG network prognostic impact

To ensure the prognostic specificity of the DMG network, the impact of DMG network-to-tumour connectivity on overall survival was assessed against across two independent cohorts of adult patients (aged >18 years) with GBM, IDH-wild type: (1) a previously reported, open-access 248-patient cohort from the University of California, San Francisco (UCSF-PDGM)⁴⁶ and (2) a newly generated institutional 272-patient cohort from the National Hospital for Neurology and Neurosurgery, Queen Square, London. Patients with midline tumours in the UCSF-PDGM cohort (n = 9) were excluded as tumours were not tested for histone alterations and a diagnosis of DMG could therefore not be wholly excluded. Previously generated tumour segmentation masks in the UCSF-PDGM data were checked against available T2-FLAIR sequences while tumour segmentation masks for patients at Queen Square were manually segmented from all available MRI sequences at diagnosis per patient, as previously described for children with DMG. GBM locations were used as seeds in the ABCD1000 connectome²⁶ to obtain brain-wide maps of tumour functional connections across both cohorts, enabling the computation of DMG network-to-tumour connectivity per patient. Survival analyses were replicated using identical methods and with the inclusion of clinical covariates which are reported to influence the overall survival of adult patients with GBM, IDH-wild type: patient demography (age at diagnosis and sex); O⁶-methylguanine-DNA methyltransferase (MGMT) promoter methylation; and clinical management (extent of resection, and completion of adjuvant radiotherapy and chemotherapy)⁷⁰. Furthermore, we investigated whether the prognostic effects of the DMG network reflected underlying common structure within the paediatric functional connectome⁴⁷. We first computed the row-summation vector of the group-average connectivity matrix of the ABCD1000 connectome (∑C), using previously reported methods⁴⁷. Next, we generated a ∑C-residual DMG network by regressing out ∑C from each discovery cohort patient T-map and reperformed tumour network mapping as previously described. ∑C-to-tumour and ∑C-residual DMG network-to-tumour functional connectivity was computed as previously described. To test the prognostic impact of ∑C alone and the DMG network with ∑C regressed out, respectively, connectivity scores were then used for Kaplan–Meier and univariable and multivariable Cox proportional hazards survival analyses. Spatial similarities of ∑C and the ∑C-residual DMG network to the primary result DMG network were tested using both spin and label permutation tests.

Evaluating the role of the DMG network in tumour anatomic progression

Tumour progression mapping

We assessed whether tumour anatomic progression mapped to DMG network trajectories. A subcohort of 71 individuals who underwent brain MRI within 2 months of death were identified; final follow-up tumours were segmented; and brain MRIs were preprocessed in an identical manner to that described above. Segmented tumour progression maps, not including baseline tumour volume, were combined across this subcohort to generate a frequency map of new tumour growth. Patterns of tumour progression were qualitatively reviewed by a central panel of three board-certified paediatric neuroradiologists (K.M., A.B. and S.S.) and compared to previously published post mortem series^15,17. The volume of new tumour growth within the DMG network and out of the DMG network was also computed and compared using a Wilcoxon signed-rank test. This subcohort was then further refined to 34 individuals whose tumour exhibited progression beyond the primary tumour site to other brain areas. To test the hypothesis that tumours infiltrate more functionally connected DMG network nodes, the degree of concordance between functional connectivity and tumour burden per node was assessed across pre-defined DMG network nodes for all 34 individuals using methods adapted from network diffusion modelling⁹⁵. DMG network nodes, or peak clusters, were defined as previously described using FSL cluster and brain regions representative of each node generated with paediatric MINCI 60 MB SynthSeg masks⁶⁵. The voxelwise degree of tumour overlap with each of the 18 DMG network nodes was computed for each tumour progression map, defined as the proportion (percentage overlap) of each DMG network node that intersected each tumour. Simultaneously, mean node-to-tumour functional connectivity values were computed from each patient-specific functional connectivity map. These mean functional connectivity values were then distance-corrected by division by the Euclidean distance (in mm) from the centroid of each DMG network node to each baseline tumour. Intraindividual Spearman’s correlation of distance-normalized functional connectivity values and tumour burden per node were computed and averaged for group-level inference. Group-level 95% CIs were obtained by bootstrapping with 2,000 iterations.

Spinal DMG cohort

An independent cohort of six children with intracranial parenchymal (non-leptomeningeal) metastases secondary to a biopsy-confirmed diagnosis of primary spinal DMG, H3K27-altered, were identified from a larger international cohort of 36 children with primary spinal DMG. Intracranial metastatic tumours were preprocessed, segmented and transformed to MNI space as previously described before anatomic comparison with the DMG network. The volume of metastatic disease within the DMG network and out of the DMG network was also computed and compared using a Wilcoxon signed-rank test.

Evaluating the spatiotemporal relationship between DMG incidence and brain metabolism

Modelling DMG incidence with [¹⁸F]FDG-PET data

DMG incidence was extracted from frequency of diagnosis at each age by rounding age at diagnosis across the whole study cohort down to the nearest year. Age-averaged [¹⁸F]FDG-PET data were obtained from a population-level paediatric [¹⁸F]FDG-PET atlas comprised of 795 children with inclusion criteria and acquisition parameters as previously reported⁵³. The voxelwise mean-squared difference between age-consecutive whole-brain standardized uptake values was computed, reflecting the mean difference in [¹⁸F]FDG-PET uptake between consecutive ages, and plotted against DMG incidence data across childhood (0–16 years) rounded to the nearest whole year. A linear regression model was fit between age-wise mean-squared difference and whole-cohort age incidence to assess the variance explained by model fit. 95% CI for the variance in age incidence data explained by [¹⁸F]FDG-PET signal changes was evaluated using 2,000 bootstrap iterations. Next, to assess the spatial specificity of the linear relationship between [¹⁸F]FDG-PET signal and age incidence, age-averaged [¹⁸F]FDG-PET data were parcellated according to paediatric MINCI 60 MB SynthSeg anatomic masks (n = 46). Brain region anatomic masks were defined as within the DMG network if ≥50% of voxels in each mask overlapped with those in the DMG network; brain regions not meeting this criterion were categorized as out of the DMG network. Region-wise and age-wise [18 F]FDG-PET signal changes within and out of DMG network were then compared in relation to their linear regression model fit with age incidence using a two-tailed Mann–Whitney U test with Benjamini–Hochberg correction for multiple comparisons.

Spatiotemporal specificity of peak childhood [¹⁸F]FDG-PET changes

[¹⁸F]FDG-PET maps were compared to identify voxelwise differences between consecutive ages. The age at which the maximal difference in normalized [¹⁸F]FDG-PET signal occurred across all ages was identified per voxel and assigned this age value. Voxelwise spatial maps of age-wise peak metabolic changes were thereby generated and compared for age-wise spatial similarity with the DMG network with label permutation and spin permutation tests, as described.

Neuropathological validation and chemoarchitectural characterization of the DMG network

Patient tumour tissue collection and dissociation for snRNA-seq

Tumour tissue samples (n = 7) collected at the time of primary biopsy or resection were snap-frozen and processed for snRNA-seq analysis under BRAIN UK approval⁹⁶. Tumour tissue was placed in 1.0 ml prechilled nucleus isolation medium containing 10 mM Tris-HCl pH 8.0, 250 mM sucrose, 25 mM KCl, 5 mM MgCl₂, 0.1% Triton X-100, 0.1 mM DTT, 1× protease inhibitor cocktail (Promega) and 1% Protector RNase Inhibitor (Merck)⁹⁷. Tumour tissue was then homogenized using a Wheaton Douce homogenizer (10 strokes loose pestle followed by 10 strokes fine pestle). Homogenized tissue was filtered through a 40-µm cell strainer and centrifuged (900g) for 10 min at 4 °C. The supernatant was aspirated to a residual 200 µl volume without disturbing pelleted nuclei and resuspended to total volume 1.0 ml using resuspension buffer (PBS with nuclease-free bovine serum albumin (1%) containing 0.5% Protector RNase Inhibitor). The samples were then centrifuged under the same conditions a second time and the supernatant aspirated to a residual 200 µl volume without disturbing pelleted nuclei which were then stained with DAPI (0.1 µg ml⁻¹) and resuspended to 500 µl total volume for sorting. All of the procedures were performed using aseptic technique on ice or at 4 °C.

Single-nucleus sorting was performed using a BD FACSymphony S6 cell sorter running FACSDiva v.9.0 using a 355 nm laser (band pass filter 450/50, long pass filter 410). Unstained and single-stained controls were included for all tumours. A standard gating strategy was applied to all of the samples (Extended Data Fig. 9a). Nuclei were identified by positive staining for DAPI. Doublets were excluded by stringent singlet-gating based on forward scatter height (FSC-H) versus width (FSC-W) and side scatter height (SSC-H) versus width (SSC-W). Singlet nuclei were sorted, using a 70-µm nozzle, into PCR tubes containing 0.2 µl prechilled resuspension buffer and centrifuged (200g) for 1 min at 4 °C twice before counting and processing for whole-transcriptome amplification, library preparation and sequencing.

snRNA-seq data generation

snRNA-seq libraries were generated using the 10x Chromium Controller and Chromium Next GEM Single Cell 3’ kit v.3.1 according to the manufacturer’s instructions (CG000315 Rev F; https://cdn.10xgenomics.com/image/upload/v1722285481/support-documents/CG000315_ChromiumNextGEMSingleCell3__GeneExpression_v3.1_DualIndex__RevF.pdf). Generated libraries were sequenced using NextSeq 2000 P3 XLEAP-SBS Reagent Kits (Illumina) at targeted 20,000 reads per nucleus.

snRNA-seq data analysis

Sequencing fastq files were processed with cellranger v.8.0.1, using the GRCh38 2024-A reference obtained from 10x Genomics (https://www.10xgenomics.com/support/software/cell-ranger/downloads#reference-downloads), and with the chemistry flag set to ‘threeprime’ and other arguments set as default. The resulting filtered matrix h5 files were analysed in R using Seurat v.5.3.1. Quality-control metrics for mitochondrial and ribosomal genes were calculated using the function PercentageFeatureSet with regular expression patterns ‘^MT-’ and ‘^PR[SL]’, respectively. Cells were filtered to have at least 1,000 and no more than 50,000 UMI counts and no more than 10% total counts attributed to mitochondrial genes.

Copy-number estimation was performed on filtered count matrices using infercnv v.1.26.0 and SCEVAN v.1.0.3. In both pipelines, a reference dataset was used comprising normal brain expression data from GSE168408, selecting only samples from individuals aged <21 years⁹⁸ such that 500 each of neurons, astrocytes, oligodendrocyte precursor cells, oligodendrocytes, endothelial, and vascular leptomeningeal cells, and 1,000 CD45⁺ immune cells were retained at random.

Cell type assignment was performed in Azimuth v.0.5.0 using the human motor cortex reference and filtering for class confidence score >0.5. Additional cell type assignment was performed using specific gene markers: BCAN and PDGFRA for malignant cells; PTPRC (also known as CD45) for immune cells; CD68 as a general macrophage and microglia marker; and P2RY12 and TMEM119 for microglia. Detection of lymphocytes was attempted using CD3E, CD4 and CD8A but no significant expression was observed. Oligodendrocytes were confirmed using MOBP and MOG expression; astrocytes were confirmed using expression of AQP4.

Normalization and variable feature detection were performed using the function SCTransform with variable features set to 2,000. Dataset integration was performed using Harmony v.1.2.4 with theta = 2, sigma = 0.1, dims.use = 1:30 and lambda estimation determined empirically by the algorithm. UMAP dimensionality reduction was performed on the first 20 dimensions of the Harmony output. Graphical clustering was obtained using FindNeighbours k.param = 20 and FindClusters resolution = 0.2.

Stemness score and differentiation score for triplot projections were calculated using methods and gene signatures from a previous study⁵². Stemness scoring was performed by obtaining the maximum of either astrocytic (AC) or oligodendrocytic (OC) state. If the higher score was AC, it was multiplied by −1. If either score was ≤0, the value was set to 0 with jitter. Differentiation scoring was obtained through subtraction of the maximum of either AC or OC from the OPC score. The expression of previously defined gene signatures in high- and low-connectivity DMG malignant cells was compared: neuronal¹², invasivity¹², DIPG synaptic score³, and NLGN3-upregulated⁸. The median difference in expression between high- and low-connectivity DMG per signature was computed and evaluated using a two-tailed Wilcoxon signed-rank test with Benjamini–Hochberg correction for multiple comparisons.

DNA methylation profiling

For a subset of individuals treated at GOSH or CHCO (n = 29), DNA was extracted from primary resection-derived 50-μm formalin-fixed paraffin-embedded (FFPE) tissue, taken either as 5 μm × 10 μm rolls or as macrodissected sections to enrich for tumour content, using the Maxwell 16 FFPE Tissue LEV DNA Purification Kit on a Maxwell 16 Research Instrument (Promega), consistent with the manufacturer’s instructions. Up to 500 ng of eluted DNA was subjected to bisulphite conversion using the Zymo EZ DNA Methylation-Gold Kit (ZymoResearch). Bisulphite-converted DNA was then further treated using the Infinium FFPE DNA Restore Kit before being assayed on the Illumina Infinium Human MethylationEPIC (850k) BeadChip platform, according to the Infinium HD FFPE Methylation Assay automated protocol (Illumina). DNA methylation data were also acquired for patient tumours included in the HERBY clinical trial (n = 7) with Illumina 450k DNA methylation arrays performed as previously reported⁹⁹. DNA methylation array results were submitted to the molecular neuropathology methylation classifier v.12.8 hosted by the German Cancer Research Centre (DKFZ) to verify tissue primary histopathological diagnosis¹⁰⁰. Individuals with a methylation class-specific calibrated score ≥0.90 were included, consistent with current clinical recommendations^100,101.

DNA methylation data preprocessing

IDAT files corresponding to raw DNA methylation array data were imported using the minfi package (v.1.48.0) in R and preprocessed using the preprocessIllumina function¹⁰². The samples in which >10% of probes failed to hybridize were designated suboptimal and excluded from subsequent analyses. Probes mapping to sex chromosomes, within 50 bp of a single-nucleotide polymorphism, with at least one cross-reactive target, or with a minor allele frequency >5%, were also excluded¹⁰³. Missing CpG beta values were then imputed using the inpute.knn function. CpG intensities were finally converted into methylation beta values.

Differential methylation analyses

Probes differentially methylated in neuroimaging-defined high- and low-connectivity groups using a median cut-off were identified using the dmpfinder function within minfi¹⁰². The resulting dataframe was filtered to satisfy the following cut-offs: Benjamini–Hochberg-corrected P < 0.05 and absolute beta value difference >0.25. Hypermethylated and hypomethylated probes in this dataset were annotated with CpG-derived gene annotations from the Illumina EPIC Manifest file (v.1.0 B5). All gene region feature categories were retained. We restricted differentially methylated CpG sites to previously defined gene sets implicated in high-grade glioma growth and assessed using our snRNA-seq dataset, namely ‘neuronal genes’, ‘DIPG synaptic genes’ and ‘invasivity genes’^3,12. The differential expression of these gene sets was statistically compared using a two-tailed Mann–Whitney U test.

Next-generation DNA panel sequencing

For a subset of children treated at GOSH (n = 26), tumour DNA was screened against a curated panel of genes defined as clinically actionable or reported as recurrently altered across paediatric solid tumours¹⁰⁴. Before sequencing, DNA was assayed using a Qubit 2.0 fluorometer (Thermo Fisher Scientific) and TapeStation 2200 (Agilent) to determine the quantity and degree of fragmentation. Library preparation, sequencing and variant calling were performed as described previously¹⁰⁴. DNA panel sequencing data were also collected for children included in the PNOC (n = 38), HERBY (n = 12)²⁴ and CHCO (n = 11)¹⁰⁵ datasets as previously reported.

Chemoarchitectural dominance analyses

The chemoarchitectural profile of the DMG network was investigated by using publicly available data representing the receptor and transporter densities of 18 unique neurotransmitter receptors and transporters derived from a PET atlas of nine human neurotransmitter systems (https://github.com/netneurolab/neuromaps)^55,56,106. These data were used to model the voxelwise brain-wide distribution of human neurotransmitter receptors and transporters and to identify which neurotransmitter receptor and transporter density data best characterized the DMG network^55,56. To this end, the DMG network and 18 neurotransmitter receptor and transporter density maps were resampled to 3 mm voxel space and masked to include only grey matter brain regions, ensuring consistent resolution and data availability across maps. A multiple linear regression model was then fit to the neurotransmitter receptor and transporter density data to predict connectivity strength across the DMG network. The dominance of each neurotransmitter density profile to this overall model fit was then assessed by fitting this same regression model to every possible combination of input variables. The total dominance for each neurotransmitter was reported, which is defined as the relative increase in adjusted R² when adding the same neurotransmitter data to each input variable combination, and represents the relative contribution of each neurotransmitter receptor and transporter distribution to the distribution of DMG network connectivity strength⁵⁵. To compare the connectivity profile of short-term and long-term survivors to brainstem neuromodulatory nuclei, the mean nucleus-to-tumour connectivity values were computed and compared using a one-way ANOVA with Tukey’s post-hoc test (factor neurotransmitter system with four levels: cholinergic, serotonergic, noradrenergic and dopaminergic).

Modelling the prognostic impact of extent of resection

A subcohort of 23 children with subtotal resections of primary thalamic DMG and brain MRI data available for analysis were identified across our study cohort. Two binary masks per patient were then manually segmented on high-resolution structural postoperative brain MRI acquired within 48 h of primary surgical resection: (1) residual tumour and (2) resection cavity. We hypothesized that surgical resections prioritizing cytoreduction and disconnection of infiltrated tumour networks would be associated with improved patient overall survival. To derive a biomarker of surgical disconnection (the disconnectome), delineated surgical paths were overlayed onto the DMG network and mean connectivity values for both the resected and non-resected masks were extracted as previously described^31,88. Patients were stratified by their resection ratio, defined as the mean network-to-tumour connectivity of the resected tumour divided by the mean network-to-tumour connectivity of the non-resected tumour⁸⁸, before inclusion in a Kaplan–Meier survival model. Univariable and multivariable Cox proportional hazards models in addition to statistical mediation analyses were used to examine whether the degree of connectomic resection was independent of tumour volume resected. The volume of tumour resected in each survival class was compared using a two-tailed Mann–Whitney U test.

Statistical analyses

Data are reported in line with the Strengthening Reporting of Observational Studies in Epidemiology (STROBE) statement¹⁰⁷. Modelling results are reported in line with the Transparent Reporting of a Multivariable Prediction Model for Individual Prognosis or Diagnosis (TRIPOD) statement¹⁰⁸. Data analyses and visualization were performed using the programming languages Python (v.3.12.1; Python Software Foundation); R (v.4.3.2; R Foundation for Statistical Computing); and MATLAB (v.R2024a; MathWorks). Statistical methods were independently verified by a consulting statistician (D.R.). The normality of continuous variables was tested using histogram visualization and the Shapiro–Wilk test. Normally (Gaussian) distributed continuous variables are reported as the mean and s.d. Non-normally distributed continuous variables are reported as the median and IQR. Categoric variables are reported descriptively as percentage frequencies and compared using the χ² or Fisher’s exact test, as appropriate. Pairwise differences in parametric data were examined using an unpaired two-tailed Student’s t test or one-way ANOVA with Tukey’s post-hoc test. Paired two-tailed Student’s t tests were performed for intraindividual pairwise comparisons. For nonparametric data, a two-tailed unpaired Mann–Whitney U test or one-tailed Wilcoxon signed-rank test was performed. Sample sizes are reported for all analyses. Unless otherwise stated, hypotheses were two-tailed and P < 0.05 was considered to be statistically significant with adjustment for multiple comparisons made where appropriate.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.