## Participants

This study was conducted as part of an ongoing project, investigating brain injury in critically born neonates. Participants were recruited between April 2021 and August 2022. Participants with HIE were recruited from the NICU at the Children’s Hospital of Southwestern Ontario, London, Canada. Eligibility criteria for the patients with HIE included: HIE diagnosis, gestational age ≥ 36 weeks, birth weight ≥ 2000 g. HIE diagnosis was based on a cord gas pH of ≤ 7.0 and/or base deficit of ≥ 16 mmol per L; if pH was between 7.01 and 7.15 or a base deficit was between 10 and 15.9 mmol per L, additional history of an acute perinatal event and an APGAR score at 10 minutes of ≤ 5, or need for assisted ventilation/resuscitation at birth, and the presence of seizures or evidence of moderate or severe encephalopathy based on a standardized neurological examination26. We excluded neonates if they had evidence of major anomalies of the brain or other organs, congenital infections (e.g., TORCH), intrauterine growth restriction (IUGR), identifiable metabolic disorder or ultrasound evidence of a large parenchymal hemorrhagic infarction.

Term-born healthy newborns with no reported brain injury were also recruited to the study from the Mother Baby Care Unit (MBCU). Inclusion criteria were admission to the MBCU, birth > 36 weeks’ GA, and inborn. Exclusion criteria were evidence of congenital malformation or syndrome, antenatal infections, antenatal exposure to illicit drugs, small for gestational age (SGA) and IUGR.

The study was approved by the Health Sciences Research Ethics Board at Western University. Informed consent was provided by the parents/caregivers of the newborns. The study was conducted in accordance with the Declaration of Helsinki.

## Clinical and Demographic Variables

Maternal and newborn health data were extracted from electronic medical records by a Paediatric Nurse, Paediatric Resident or NICU Fellow. For participants with HIE, the demographic data extracted included gestational age, birth weight, biological sex, HIE stage (based on Sarnat staging27), resuscitation details, Apgar scores and cord pH. We also collected the following data: 72 h treatment with TH, the presence of brain injury on MRI, and postnatal infections (clinical sepsis or positive culture infection, confirmed necrotizing enterocolitis).

## MRI Image Acquisition

Newborns with HIE underwent at least one MRI scan post cooling after rewarming. Healthy newborns did not undergo MRI scanning. Newborns with HIE were scanned on a 1.5 T GE MRI scanner. A T1-weighted structural image was acquired (TR = 8.4–11.5 ms [depending on clinical requirements], TE = 4.2 ms, flip angle = 12/25°, matrix size 512 × 512, 99–268 slices, voxel size typically 0.39 × 0.39 × 0.5 mm (0.31 × 31 × 5 to 0.43 × 0.43 × 0.6 for some neonates), as well as a T2-weighted structural image (TR = 3517–9832 ms, TE = 7.3–8.4 ms, flip angle = 90/160°, matrix size 256 × 256, 19–60 slices, 0.7 × 0.7 × 2–5 mm voxel resolution). Additionally, an echo planar imaging sequence to measure Blood Oxygen Level–Dependent (BOLD) fMRI data was also acquired to examine RSFC (TR = 3000ms, TE = 50 ms, flip angle = 70°, matrix size 64 × 64, 39 slices, voxel size 3 x 3 x 3 mm, total volumes 35).

## fMRI preprocessing and analysis

The fMRI data were preprocessed using FSL (FMRIB Software Library28, https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FSL). The preprocessing pipeline included brain extraction, motion correction, spatial smoothing (full width at half maximum = 5 mm), band-pass filtering (0.01–0.2 Hz) and registration to a neonatal atlas29. From the frontal, parietal, temporal and occipital lobes of both hemispheres, we extracted the average BOLD sequences. Then, we calculated the Pearson correlation among them to build an 8-by-8 RSFC map for each neonate (Fig. 1B).

## T1-weighted image processing and segmentation

The T1 weighted scans were automatically segmented using infant FreeSurfer30. Automatic processing steps included intensity normalization, skull stripping, and segmentation of the cortex, white matter and subcortical structures30. A multi-atlas approach was employed for the segmentation. Multiple brain atlases of newborns were first registered to native space, and structure labels were transferred. The atlases were developed from neonatal MRI scans31. To initially create the atlases, manually segmented labels were developed using MRI scans from a representative sample of neonates (0–2 years of age). In the current study, developmentally appropriate atlases for newborns were employed. Labels were fused into a single segmentation result, providing higher accuracy than single-atlas approaches32. Volumetric measurements for anatomical features could then be extracted30. Brain subcortical brain structure volumes were extracted, as well as cerebral white matter and cerebral cortex volumes, to compute total cerebral volumes (TCV). Each segmented T1-weighted image was visually inspected using Freeview software, available within FreeSurfer. Manual segmentation was employed to correct any segmentation errors (i.e., partial volume effects) using ITK-SNAP (http://www.itksnap.org/).

## Brain injury characterization

A pediatric neuroradiology fellow (LT) scored the T1-weighted anatomical images for brain injury severity. White matter injury (WMI) was defined as foci exhibiting T1 hyperintensity without T2 hypointensity or by low-intensity T1 foci and was scored on a 3-point scale (none = 0, minimal = 1, moderate-severe = 2–3 combined) using the methods of de Vries33. Intraventricular haemorrhage (IVH) was graded (none = 0, mild = 1–2, and moderate-severe = 3–4) using Papile’s method34. Only supratentorial injuries were scored.

## fNIRS acquisition and analysis

fNIRS was acquired in newborns with HIE and healthy newborns using the same experimental setup. Participants were scanned using a NIRSport2 (NIRx, Berlin, Germany) system with 8 light sources and 8 detectors covering the whole brain. The system operates at two wavelengths of light (760 nm and 850 nm), and the sampling rate was 10.1725 Hz. For either hemisphere, 4 channels were located on the temporal, 2 on the parietal, frontal and occipital lobes (Fig. 1A). fNIRS scans lasted for a minimum of 6 minutes during rest or natural sleep with the neonates lying in the incubator, cot, or caregiver’s arms.

After data acquisition, Homer3 software35 was used for preprocessing. The pipeline included spline interpolation for motion correction36, band-pass filtering of 0.01–0.1 Hz and conversion to HbO and Hbr fluctuations with modified Beer-Lambert law37. For each neonate, two connectivity maps, sized differently, were calculated using Pearson correlation. Channels corresponding to one lobe were averaged and then correlated with other cortices to build lobe-wise connectivity maps (8 by 8). 20-by-20 channel-wise RSFC maps were also calculated with individual channels corresponding with each other. HbO and Hbr data were both used for building connectivity maps separately (Fig. 1A).

## Statistical Analysis

Statistical analyses were performed using Matlab (R2020b, Natick, Massachusetts: The MathWorks Inc).

For the first aim, which was comparing the RSFC maps yielded from the fNIRS and fMRI, Euclidean distances and Jaccard distances were calculated between the two 8-by-8 maps at various levels of sparsity38. It is considered common practice to filter out negative or low-weighted connections when analyzing RSFC maps, where the percentage of connections remaining is the level of sparsity of a map. For weighted maps, we kept the original weights of connections, while for binarized maps, the remaining connections after filtering were rounded to 1. Given adjacency matrices of two maps \({A}_{1}=\left[{a}_{ij}^{1}\right]\) and \({A}_{2}=\left[{a}_{ij}^{2}\right]\), where \({a}_{ij}\) was the weight of connection between node \(i\) and \(j\), Euclidean distance between the two was defined as

$${d}_{E}=\sqrt{{\sum }_{{ }_{i,j}^{ }}{\left({a}_{ij}^{1}-{a}_{ij}^{2}\right)}^{2}}$$

1

,

while Jaccard distance was defined as

$${d}_{J}=1-J, \text{w}\text{h}\text{e}\text{r}\text{e} J=\left\{\begin{array}{c}\frac{{\sum }_{i,j}\text{min}\left({a}_{ij}^{1},{a}_{ij}^{2}\right)}{{\sum }_{i,j}\text{max}\left({a}_{ij}^{1},{a}_{ij}^{2}\right)}, if {\sum }_{i,j}\text{max}\left({a}_{ij}^{1},{a}_{ij}^{2}\right)>0\\ 1, otherwise\end{array}\right..$$

2

We also calculated similarity maps for weighted and binarized RSFC maps between fNIRS and fMRI. Note that similarity maps had the same dimension as connectivity maps. For weighted cases, entries of similarity maps were the average Euclidean distance among subjects for the corresponding connections. In contrast, for binarized cases, they were the percentage of subjects sharing the connection. Volume counts of cortical grey matter obtained from T1-weighted images were associated with overall fNIRS and fMRI connectivity strength, respectively, using linear regression, adjusted for GA at birth, postmenstrual age (PMA) at scan and sex.

To address the second aim, which was identifying altered connectivity patterns between HIE neonates and healthy controls, we used graph theory-based measurements to quantify the patterns of 20-by-20 RSFC maps. On HbO and Hbr maps, we calculated the clustering coefficient of a map as

$${C}_{p}=\frac{1}{N}{\sum }_{i}\frac{{\sum }_{j,k}{a}_{ij}{a}_{ik}{a}_{jk}}{\left({\sum }_{j}{a}_{ij}-1\right){\sum }_{j}{a}_{ij}}$$

3

,

where \(N\) was the number of nodes in the map. Network efficiency was also calculated locally as

$$NE=\frac{1}{N}{\sum }_{i}\left[\frac{1}{\left|{G}_{i}\right|\left(\left|{G}_{i}\right|-1\right)}{\sum }_{k\ne j\in {G}_{i}}\frac{1}{{d}_{jk}}\right]$$

4

,

where \({d}_{ij}\) was the shortest distance \({N}_{i}\)between \(i\) and \(j\), \({G}_{i}\) the node set containing \(i\) and its direct neighbors only, and \(\left|{G}_{i}\right|\) the number of nodes of \({G}_{i}\). Modularity was defined as

$$M={\sum }_{p}\left[{e}_{pp}-{\left({\sum }_{q}{e}_{pq}\right)}^{2}\right]$$

5

,

where \(C\) was the number of communities identified by a Louvain-like algorithm39 and \({e}_{pq}\) denoted the fraction of connections between community \(p\) and \(q\). Based on clustering, the between-subject variation was calculated using mutual information40. Given a confusion matrix \(C=\left[{c}_{wv}\right]\), where \({c}_{wv}\) was the number of nodes in both community \(w\) and community \(v\) of two separate clusterings, respectively, mutual information was defined as

$$MI=\frac{-2{\sum }_{w,v}({c}_{wv}\text{ln}\frac{{c}_{wv}N}{{\sum }_{w}{c}_{wv}{\sum }_{v}{c}_{wv}}) }{{\sum }_{w}\left({\sum }_{v}{c}_{wv}\text{ln}\frac{{\sum }_{v}{c}_{wv}}{N}\right)+{\sum }_{v}\left({\sum }_{w}{c}_{wv}\text{ln}\frac{{\sum }_{w}{c}_{wv}}{N}\right)}$$

6

.

T-tests were used to determine if the two groups were significantly different over a measurement.

To further demonstrate the distinctive power of connectivity measurements separating HIE neonates from healthy controls, we trained machine learning models of support vector machine (SVM) with several local metrics, including connectivity, clustering coefficient, nodal efficiency, degree centrality and closeness centrality, with leave-one-out validation. The accuracies of classification were calculated, and receiver operating characteristic (ROC) curves were plotted. Areas under ROC curves (AUC) were also obtained.