Abstract
Therapeutic hypothermia following neonatal encephalopathy due to birth asphyxia reduces death and cerebral palsy. However, school-age children without cerebral palsy treated with therapeutic hypothermia for neonatal encephalopathy still have reduced performance on cognitive and motor tests, attention difficulties, slower reaction times and reduced visuo-spatial processing abilities compared to typically developing controls. We acquired diffusion-weighted imaging data from school-age children without cerebral palsy treated with therapeutic hypothermia for neonatal encephalopathy at birth, and a matched control group. Voxelwise comparison of fractional anisotropy (33 cases, 36 controls) confirmed microstructural alterations in widespread areas of white matter in cases, particularly in the fornix, corpus callosum, anterior and posterior limbs of the internal capsule bilaterally and cingulum bilaterally. In structural brain networks constructed using probabilistic tractography (22 cases, 32 controls), graph-theoretic measures of strength, local and global efficiency, clustering coefficient and characteristic path length were found to correlate with IQ in cases but not controls. Network-based statistic analysis implicated brain regions involved in visuo-spatial processing and attention, aligning with previous behavioural findings. These included the precuneus, thalamus, left superior parietal gyrus and left inferior temporal gyrus. Our findings demonstrate that, despite the manifest successes of therapeutic hypothermia, brain development is impaired in these children.
1 Introduction
Neonatal encephalopathy (NE), which often results from perinatal asphyxia, leads to a high risk of death or disability, including cerebral palsy (CP) (Azzopardi et al., 2014; Marlow, 2005; O’Connor et al., 2017; Robertson et al., 1989). In the UK, approximately 2.6 per 1000 live births in 2015 were affected by NE secondary to perinatal asphyxia (Gale et al., 2018). The recommended treatment for NE (National Institute for Clinical Excellence (NICE), 2010) is therapeutic hypothermia (TH), which consists of reducing the infant’s core temperature to 33.5°C for three days, commencing as soon as possible after the asphyxia (Azzopardi et al., 2009; Rutherford et al., 2010). TH reduces the chance of death and disability at 18 months (Jacobs et al., 2013), reduces likelihood and severity of CP (Jary et al., 2015) and increases the incidence of survival with an IQ > 85 (Azzopardi et al., 2014). However, recent studies have shown that children aged 6-8 years, who underwent TH at birth for NE and did not develop CP, perform worse in motor and cognitive tests than controls (Jary et al., 2019; Lee-Kelland et al., 2020) and have attention difficulties, slower reaction times and reduced visuo-spatial processing abilities (Tonks et al., 2019). These motor deficits are not predicted by 18-month developmental scores (Jary et al., 2019). Thus, despite the reduced occurrence of severe disabilities following TH, aspects of brain development remain affected by NE.
Studies on children born with NE, prior to widespread use of TH (Gao et al., 2012; Ly et al., 2015; Martinez-Biarge et al., 2012), and on animal models (Chakkarapani et al., 2010; Kyng et al., 2015; Yue et al., 1997) indicate damage to white matter and subcortical structures, caused by hypoxic-ischaemic brain injury. Some studies have shown an association between hypothermia/rewarming and subcortical white matter apoptosis independent of hypoxic-ischemic brain injury (O’Brien et al., 2019; Wang et al., 2016), whereas other findings suggest no impact of hypothermia on the subcortical white matter (Gressens et al., 2008). It is unknown how the interplay between the damage mechanisms of NE and the effects of TH impact brain development.
Diffusion-weighted imaging (DWI) provides a non-invasive tool for investigating white matter microstructure. Measurement of diffusion of water molecules through brain tissue allows calculation of diffusion metrics such as fractional anisotropy (FA), which is related to its microstructural properties. FA is affected by properties such as myelination and fibre density (Le Bihan and Johansen-Berg, 2012) and has clinical relevance in patient cohorts (Assaf and Pasternak, 2008; Dennis and Thompson, 2013a; Assaf et al., 2019). We used tract-based spatial statistics (TBSS) (Smith et al., 2006) to perform voxel-wise comparison of FA across the brain’s white matter, whilst controlling for multiple comparisons. We further investigated white matter connectivity by constructing structural brain networks, or connectomes (Sporns et al., 2005), in which nodes represent brain regions and edges were determined by probabilistic tractography. We characterised structural networks by drawing on techniques from graph theory (Bullmore and Sporns, 2009; Hagmann et al., 2010a; Fornito et al., 2013; Bassett and Sporns, 2017), allowing comparison of quantitative differences in whole-brain network structure. Such techniques have previously been used to characterise the developing human connectome (Hagmann et al., 2010b; Dennis and Thompson, 2013b; Morgan et al., 2018), as well as in the study of specific neurodevelopmental complications such as CP (Arrigoni et al., 2016; Pannek et al., 2014) and neurodevelopmental impairments following preterm birth (Brown et al., 2014; Muñoz-Moreno et al., 2016). We then used the network-based statistic (NBS) (Zalesky et al., 2012, 2010) to look for subsets of connections (subnetworks) which were weakened in cases, and subnetwork which related to measures of cognition.
Our findings demonstrate that, although TH reduces severe disabilities after NE, underlying structural deficits are present which are associated with the cognitive differences found between cases and controls at school-age. These differences are often overlooked as most children given TH for NE do not demonstrate significant deficits in performance at 18 months.
2 Materials and Methods
2.1 Participants
Informed and written consent was obtained from the parents of participants, in accordance with the Declaration of Helsinki. Ethical approval was obtained from the North Bristol Research Ethics Committee and the Health Research Authority (REC ID: 15/SW/0148).
2.1.1 Cases
Eligibility criteria were as follows: gestation at birth ≥36 weeks and treatment with TH as standard clinical care based on TOBY trial eligibility criteria including signs of perinatal asphyxia and moderate to severe encephalopathy, confirmed by amplitude integrated electroencephalogram (Azzopardi et al., 2009). Children were excluded if they had started cooling later than six hours after birth, were cooled for less than three days, had received Xenon as part of a neuroprotective feasibility study, had been found to have a metabolic or genetic disorder, or if any major intracranial haemorrhage or structural brain abnormality could be seen on the neonatal MRI scan. Cases were sequentially selected from the cohort of children who received TH between 2008 and 2011. These data are maintained by the Bristol Neonatal Neurosciences group at St Michael’s Hospital, Bristol, UK, under previous ethics approval (REC ID: 09/H0106/3). A diagnosis of CP was ruled out at 2 years and reconfirmed at 6-8 years. Children were native English speakers and had no additional medical diagnosis other than NE.
2.1.2 Controls
The control group consisted of children matched for age, sex and socio-economic status (Lee-Kelland et al., 2020). Children were excluded if they were born before 36 weeks gestation, had any history of NE or other medical issues of a neurological nature, or were not native English speakers.
2.2 Cognitive Assessment
Cognitive performance was assessed using the Wechsler Intelligence Scale for Children 4th Edition (WISC-IV) (Kaufman et al., 2006), which summarises raw score performance from 10 subsets into 10 scaled scores. These 10 scores are summed in four domains – verbal comprehension, perceptual reasoning, processing speed and working memory – which are combined to form a full-scale intelligence quotient (FSIQ) score. Cognitive testing was administered by assessors who were not previously involved with the patients’ care and were blinded to case-control status.
2.3 Image Acquisition
T1-weighted images and DWI data were acquired with a Siemens 3 tesla Magnetom Skyra MRI scanner at the Clinical Research and Imaging Centre (CRiCBristol), Bristol, UK. An experienced radiographer placed children supine within the 32-channel receive only head-coil, and head movement was minimised with memory-foam padding. Children wore earplugs and were able to watch a film of their choice. A volumetric T1-weighted anatomical scan was acquired for tissue segmentation and parcellation, with the magnetisation-prepared rapid acquisition gradient echo (MPRAGE) sequence using the following parameters: echo time (TE) = 2.19 ms; inversion time (TI) = 800 ms; repetition time (TR) = 1500 ms; flip angle = 9°; field of view (FoV) 234 × 250 mm; 240 slices; 1.0 mm isotropic voxels. DWI data were acquired for tractography and microstructural analysis, with a multiband echo-planar imaging (EPI) sequence, using the following parameters: TE = 70 ms; TR = 3150 ms; FoV 192 × 192 mm; 60 slices; 2.0 mm isotropic voxels, flip angle 90°, phase encoding in the anterior-posterior direction, in-plane acceleration factor = 2 (GRAPPA (Griswold et al., 2002)), through-plane multi-band factor = 2 (Moeller et al., 2010; Setsompop et al., 2012a, b). For the purpose of data averaging and eddy-current distortion correction, two sets of diffusion-weighted images were acquired with b = 1,000 s mm-2 in 60 diffusion directions, equally distributed according to an electrostatic repulsion model, as well as 8 interspersed b = 0 images, with one data set acquired with positive phase encoding steps, then repeated with negative steps (so-called, “blip-up, blip-down”), giving a total of 136 images.
2.4 Quality Control
The quality of the DWI data was assessed using the EddyQC tool (Bastiani et al., 2019) from the FMRIB Software Library (FSL, http://fsl.fmrib.ox.ac.uk) (Smith et al., 2004). Scans were rejected if the root-mean-square of all movement and eddy current metrics from EddyQC was greater than one standard deviation above the mean for all participants.
T1-weighted anatomical images were assessed visually; any scans with severe movement artefacts were rejected. The remaining scans were processed with the structural pipeline described below, followed by further visual inspection of the parcellation and tissue segmentation. Scans were further rejected at this stage if any moderate artefacts had caused errors in the parcellation or segmentation.
Figure 1 shows the process of recruitment and scan quality control. We recruited 51 cases and 43 controls for this study. Of these, 7 cases and 4 controls did not want to undergo scanning. A further 4 cases had incomplete data due to movement during the scan. DWI quality control led to the rejection of a further 6 cases and 2 controls. One further case and one control were rejected due to incorrect image volume placement. This left 33 case and 36 control scans which passed the DWI quality control, which were used in the TBSS analysis. Of these remaining 69 datasets, the anatomical scan for 11 cases and 4 controls was not of sufficient quality to allow segmentation and parcellation, leaving 22 cases and 32 controls for network analysis. Participant demographics are shown in Table 1.
Demographics are shown for each of the TBSS and network analysis groups. Median (range) is shown for age, socio-economic status and FSIQ. Socio-economic status is defined as follows: A = upper middle class, B = middle class, C1 = lower middle class, C2 = skilled working class, D = working class, E = casual worker or unemployed.
Flowchart of participants at each stage of quality control. FoV = field of view, indicating the scans which were rejected due to incorrect image volume placement.
Anatomical images were visually assessed for focal lesions and abnormal signal intensities. In the TBSS datasets, lesions were present in 1 case and 2 controls. In the network analysis datasets, lesions were present in 1 control. These lesions were judged by the blinded assessor (FC) to be non-severe, consequently these subjects were not excluded.
Note that previous findings from the same cohort demonstrate reduced performance in cases in all WISC-IV domains (Lee-Kelland et al., 2020), whereas in the smaller group which passed quality control in this study cases exhibit significantly reduced performance in perceptual reasoning, verbal comprehension, working memory and FSIQ. Though processing speed was reduced in cases in this study, the difference was not significant (see Table 1).
2.5 TBSS
Voxelwise statistical analysis of the FA data was carried out using TBSS (Smith et al., 2006), part of FSL. FA images were generated by fitting a tensor model to the diffusion data using FSL’s FDT software. All images were then nonlinearly registered to one subject, chosen automatically by finding the most representative subject, which was then affine registered to MNI152 standard space. This is the recommended procedure when testing data from children, which may not register well to an adult template (Smith et al., 2006). The mean FA image was then eroded to create a skeletonised representation of the white matter tracts. Each subject’s registered FA image was then projected onto this skeleton to allow voxelwise statistics.
2.6 Structural Network Construction
A weighted connectome was constructed for each subject, with nodes defined by parcellation of the anatomical scan and edges determined by probabilistic tractography using the DWI data. The processing pipeline, described in more detail below, is summarised in Figure 2.
Method for constructing structural brain networks from T1 and DWI data. Cortical and sub-cortical nodes were defined by segmentation of the T1-weighted structural scan. Edges were determined by seeding streamlines from the cortical grey/white matter interface and performing tractography using the fibre orientation distribution obtained by spherical deconvolution of the measured diffusion signal. Edges were weighted by the mean FA along all streamlines passing between the corresponding pair of nodes, and the resulting network was represented by a connectivity matrix.
2.6.1 Node Definition
T1-weighted anatomical images were denoised with the Advanced Normalization Tools DenoiseImage tool (http://github.com/ANTsX/ANTs) (Manjón et al., 2010). Brain extraction was performed with either SPM8-VBM (http://fil.ion.ucl.ac.uk/spm) (Ashburner and Friston, 2005) or CAT12 (http://www.neuro.uni-jena.de/cat) (Gaser and Dahnke, 2016) depending on which gave better delineation of the brain surface for each subject. Each subject’s T1-weighted image was parcellated, using FreeSurfer (http://surfer.nmr.mgh.harvard.edu) (Fischl, 2012), according to the Desikan-Killiany atlas (Desikan et al., 2006) (a total of 84 regions; 34 cortical, 7 subcortical and 1 cerebellar per hemisphere). The FIRST (Patenaude et al., 2011) subcortical segmentation tool from FSL was found to give better segmentation of subcortical structures (including the hippocampus and amygdala) than FreeSurfer, so this was combined with the cortical parcellation from FreeSurfer using the labelsgmfix tool from MRtrix (www.mrtrix.org) (Tournier et al., 2019).
2.6.2 Edge Definition
DWI data were corrected for eddy current induced distortions and subject movements using EDDY (Andersson and Sotiropoulos, 2016) and TOPUP (Andersson et al., 2003), from FSL. Subsequent DWI processing and tractography steps were performed using MRtrix. The response function (the DWI signal for a typical fibre population) was estimated from the data (Tournier et al., 2013) in order to calculate the fibre orientation distribution (FOD) by performing constrained-spherical deconvolution of the response function from the measured DWI signal (Tournier et al., 2007). The normalised FOD image and the five-tissue-type segmentation of the T1-weighted anatomical image were used to perform anatomically-constrained tractography (Smith et al., 2012) using second-order integration over FODs (Tournier et al., 2010), with the following parameters: step size = 1 mm, minimum length = 50 mm, cutoff FOD magnitude = 0.1, maximum angle between steps = 30°. Streamlines were seeded in the interface between grey and white matter and only accepted if they terminated in subcortical or cortical grey matter. Terminated streamlines which were not accepted were allowed to backtrack to a valid point to be resampled (Smith et al., 2012). This method was used to generate 10 million streamlines which were subsequently filtered to 1 million using spherical-convolution informed filtering of tractograms (Smith et al., 2013) in order to improve biological plausibility and remove length bias. FA images were then used to assign a weight to each streamline according to the mean FA along its path. In order to construct a weighted graph for each subject, edges were defined between any pair of nodes connected by at least one streamline, with the connection strength defined by the mean FA along all streamlines connecting the nodes.
2.7 Network Metrics
We selected the following metrics to quantify properties of the FA-weighted structural connectivity networks: average strength, characteristic path length, global efficiency, local efficiency, clustering coefficient, modularity and small-worldness. These are defined below (for an in-depth description see Rubinov and Sporns, 2010).
The strength of a node is defined as the sum of the weights of all edges connected to the node. The average weight for the entire graph is equal to the average node strength across all nodes. The characteristic path length of the graph is the average of the shortest path from each node to every other node, where the edge distances used to calculate path lengths are defined inversely to edge weights (making stronger connections equivalent to shorter paths). Note that this does not reflect physical distance between regions in the brain. A shorter characteristic path length indicates stronger connectivity across brain regions, thus implying stronger potential for integration (Rubinov and Sporns, 2010). Global efficiency is the average of the inverse of the shortest path length. This has a roughly inverse relationship with characteristic path length, and therefore indicates integration (Bullmore and Sporns, 2009). However, the two metrics differ in the edges they are influenced by; the calculation of characteristic path length is more dependent on longer paths, whereas global efficiency is more dependent on shorter paths.
Local efficiency of a given node is the average of the inverse of the shortest path length between the immediate neighbours of that node. This is then averaged across all nodes to give a single measure for the whole graph. The clustering coefficient gives the number of connections between the nearest neighbours of a node as a fraction of the maximum number of possible connections. Modularity indicates how well the network can be split up into relatively separate communities (i.e. modules) of nodes by measuring a normalised ratio of the number of within-module connections to the number of between-module connections. Local efficiency, clustering coefficient and modularity indicate the efficiency of local information transfer, thus indicating the potential for segregated functional processing (Bullmore and Sporns, 2009; Rubinov and Sporns, 2010).
Both integration and segregation are required for brain networks to carry out localised and distributed processing simultaneously (Tononi et al., 1994). The degree to which a network exhibits both segregation and integration is measured by the small-worldness of the network (Muldoon et al., 2016; Rubinov and Sporns, 2010). A high degree of small-worldness is characterised by a high clustering coefficient and low characteristic path length compared to random graphs. We measured small-worldness with small-world propensity (Muldoon et al., 2016). All other metrics were calculated with the Brain Connectivity Toolbox (http://www.brain-connectivity-toolbox.net) (Rubinov and Sporns, 2010).
2.8 Statistical Analysis
Group differences between case and control network metrics were tested using two-tailed, unpaired t-tests. Correlation of network metrics with cognitive score was then tested by calculating the partial Pearson correlation coefficient, including age and sex as covariates. In order to reduce the effect of multiple comparisons and increase statistical power, each network metric was tested for correlation with FSIQ, not with every WISC-IV domain. Bonferroni correction was applied to correct for multiple comparisons. Statistical analysis of the network metrics was performed in MATLAB (R2018b, Mathworks). For TBSS, significance was tested using FSL’s non-parametric permutation testing software, RANDOMISE (Winkler et al., 2014). We used 10,000 permutations and applied threshold-free cluster enhancement to correct for multiple comparisons. Significant results have corrected p < 0.05.
2.8.1 Network-Based Statistic (NBS)
We used NBS to test the hypothesis that cases exhibit reduced connectivity (i.e. reduced FA) compared to controls, based on previously reported findings of reduced FA in white matter in neonates treated with TH for NE (Lally et al., 2019; Tusor et al., 2012). We also explored group differences in the relationship between cognitive scores and connectivity.
NBS (Zalesky et al., 2012, 2010) is a nonparametric, permutation-based approach for controlling family-wise error rate (FWER) on the level of subnetworks. NBS identifies connected subnetworks in which each edge satisfies the given contrast (e.g. group differences in connectivity). The t-statistic is calculated for each edge in the network, then thresholded at a chosen value. Of the remaining suprathreshold edges, the size of each connected subnetwork (given by the number of edges) is stored. This process is repeated for random permutations of the data to estimate the null distribution. The FWER-corrected p-value for each subnetwork is given by the number of permutations for which the largest connected subnetwork in the permuted data is the same size or larger than the given subnetwork, normalised by the number of permutations.
We tested for reduced connectivity in cases compared to controls (one-tailed) and for group differences in the dependence of cognitive scores on edge weights (two-tailed). We tested all four cognitive domains for correlation (perceptual reasoning, processing speed, verbal comprehension, working memory) in addition to FSIQ. We used 10,000 permutations to calculate the p-value. In order to only test robust edges, only connections present in >50% of cases and >50% of controls were assessed. Age and sex were included as covariates in a general linear model in all tests (design matrices are shown in Supplementary Tables 2 and 3). As recommended in the literature (Zalesky et al., 2012, 2010), a range of t-statistic thresholds were tested (2.5-3.5) to find the value which gave robust results (Supplementary Figure 2). This procedure allows identification of large subnetworks with subtle effects (at low primary thresholds) as well as smaller subnetworks with strong effects (at high primary thresholds). Significant results have p < 0.05 (FWER-corrected).
2.9 Visualisation
Subnetworks were visualised with the BrainNet Viewer (https://www.nitrc.org/projects/bnv/) (Xia et al., 2013) and as Circos connectograms (http://www.circos.ca) (Krzywinski et al., 2009).
2.10 Data Availability
The data that support the findings in this article are available upon reasonable request to the corresponding author.
3 Results
3.1 TBSS
Figure 3 shows the results of voxelwise comparison of FA using TBSS, demonstrating widespread reduction in FA in cases compared to controls. The effect is most prominent in the fornix, the corpus callosum, anterior and posterior limbs of the internal capsule bilaterally, and the cingulum bilaterally, but can also be seen in other distributed areas of white matter. These results demonstrate extensive alterations to white matter microstructure in cases. This analysis was repeated with age and sex included as covariates in a general linear model; the results were largely unchanged (see Supplementary Figure 1).
Results of voxelwise comparison of FA on the white matter skeleton (green). Significant results are indicated by the colour bar (p < 0.05, TFCE corrected). These are overlaid on the MNI standard template. Labels indicate some major white matter tracts and regions. Abbreviations are as follows: anterior thalamic radiation (ATR), corpus callosum (CC), external capsule (EC), internal capsule (IC), inferior fronto-occipital fasciculus (IFOF), inferior longitudinal fasciculus (ILF), superior longitudinal fasciculus (SLF), uncinate fasciculus (UF). Differences are most notable in the IC, CC, fornix and cingulum.
3.2 Network Metrics
3.2.1 Group Differences
No significant group differences were found in network metrics (see Supplementary Table 1). Notably, small-world characteristics were expressed robustly across the entire cohort with all subjects expressing a small-world propensity greater than 0.82 (networks with small-world propensity > 0.6 are considered small-world (Muldoon et al., 2016)).
3.2.2 Cognitive Correlations
Figure 4 shows the correlation of network metrics with FSIQ. In cases, FSIQ was significantly correlated with average node strength (r = 0.6858, p = 0.0059), local efficiency (r = 0.6320, p = 0.0196), global efficiency (r = 0.6672, p = 0.0092), clustering coefficient (r =0.6817, p = 0.0065) and characteristic path length (r = −0.6704, p = 0.0085). In controls, network metrics exhibited the same general trends as in cases, however none of the correlations were significant, despite there being a comparable spread in the residuals.
Network metrics and FSIQ were controlled for age and sex, with residuals plotted for cases (blue circles) and controls (red triangles). These are fitted with a blue solid line and red dashed line for cases and controls respectively. Where p > 0.05, plots are labelled as not significant (n.s.). p-values are Bonferroni corrected for the number of correlations performed.
3.3 NBS
Figures 5 and 6 show the significant subnetworks identified by NBS. To reiterate; in each of the subnetworks, the tested contrast is expressed significantly at the level of each individual connection, with FWER controlled for the whole subnetwork. Significant results were found for: reduced connectivity (equating to reduced FA) in cases compared to controls; stronger relationship between connectivity and FSIQ in cases than in controls; and stronger relationship between connectivity and processing speed in cases than in controls. No results were found for group differences in the relationship between connectivity and perceptual reasoning, verbal comprehension or working memory.
Subnetworks are shown for case-control comparison (top), correlation with FSIQ (middle) and correlation with processing speed (bottom). In the case-control status subnetwork, all connections shown are significantly weakened in cases compared to controls (reflecting lower FA in cases). In the FSIQ and processing speed subnetworks, the dependence of cognitive score on connection strength is significantly higher in cases than in controls. The dorsal (axial) view shows all connections, while the lateral (sagittal) views of the left and right cortices show the intrahemispheric connections.
Subnetworks given by NBS analysis of case-control comparison (A), correlation with FSIQ (B) and correlation with processing speed (C). Connectograms are shown with interhemispheric connections in green and intrahemispheric connections in red (left) and blue (right). Panel A also shows box plots of the mean FA across all connections in the case-control subnetwork for both the true FA values (left) and residual values after controlling for age and sex (right). In the box plots, the circle is the median, the solid box represents the 25th to 75th percentiles, and the lines extend to the minimum and maximum data points. Panels B and C also show scatter plots of the mean FA across all connections in the FSIQ and processing speed subnetwork, respectively, for both cases (blue circles, blue solid line) and controls (red triangles, red dashed line). The complete list of node label abbreviations is shown in the Supplementary Table 7.
Connectivity was significantly reduced in cases compared to controls (t = 2.8, p = 0.0304) in a subnetwork comprising 19 nodes (10 left, 9 right) and 20 edges (14 interhemispheric, 6 intrahemispheric). In this subnetwork, the five most well-connected nodes were the right precuneus cortex, left superior parietal gyrus, left precuneus cortex, left thalamus and left inferior temporal gyrus.
The relationship between connectivity and FSIQ was significantly stronger in cases than controls (t = 3.5, p = 0.0132) in a subnetwork comprising 23 nodes (10 left, 13 right) and 22 edges (all interhemispheric). The subnetwork was entirely composed of interhemispheric connections connecting parietal, limbic, temporal and occipital areas. The five most well-connected nodes were the left precuneus cortex, left and right supramarginal gyrus, left superior parietal gyrus and right parahippocampal gyrus.
The relationship between connectivity and processing speed was significantly stronger in cases than controls (t = 3.3, p = 0.0122) in a subnetwork comprising 28 nodes (6 left, 22 right) and 30 edges (7 interhemispheric, 23 intrahemispheric). The three most well-connected nodes were the right lingual gyrus, left superior parietal gyrus and right cuneus cortex. See Supplementary Tables 4-6 for the complete list of nodes in each subnetwork.
To provide graphical demonstration of each effect, the average FA in each of these subnetworks was calculated for each subject and plotted as a box plot for case-control differences (Fig. 6A) and plotted against FSIQ (Fig. 6B) and processing speed (Fig. 6C). This figure clearly demonstrates the effect captured by each subnetwork. The median of the subnetwork-averaged FA in cases was 6% lower in cases than in controls (p < 0.0001) in the case-control status subnetwork. The dependence of FSIQ on connectivity was much stronger in cases than controls in the FSIQ subnetwork (ANCOVA with age and sex as covariates; p < 0.0001). Similarly, the dependence of processing speed on connectivity was stronger in cases than controls in the processing speed subnetwork (ANCOVA with age and sex as covariates; p < 0.0001).
Removing the outlier in Figure 6B strengthened the correlation (r = 0.7860, p < 0.0001).
Discussion
This study assessed white matter microstructure and connectivity properties in children aged 6-8 years who underwent TH for NE at birth and who did not develop CP, compared to a matched group of control children with no history of neurological issues. TBSS was used to compare white matter microstructural properties, derived from diffusion weighted imaging data at the voxel level, between cases and controls. Network analysis was used to further investigate the relationship between brain connectivity and cognitive measures in cases and controls, using graph theory to interpret connectome data. NBS was used to determine the specific connections associated with case-control status and those associated with cognitive performance.
Children who were treated with TH for NE at birth exhibited widespread reduction in FA compared with controls. Correlations with FSIQ were found in strength, local efficiency, global efficiency, clustering coefficient and characteristic path length of the structural networks, in cases only. NBS revealed subnetworks associated with case-control status, FSIQ and processing speed.
4.1 Cases Exhibit Widespread Alterations to White Matter Microstructure
Several factors can cause a reduction in FA, including reduced fibre density, cross-sectional area or myelination. Previous studies of neonates treated with TH for NE have investigated the relationship between white matter diffusion properties, measured in the first weeks following birth, and neurodevelopmental outcome at 2 years of age; these studies found a significant reduction in FA in infants with adverse outcomes, compared to those with favourable outcomes, in widespread areas of white matter including the centrum semiovale, corpus callosum, anterior and posterior limbs of the internal capsule, external capsules, fornix, cingulum, cerebral peduncles, optic radiations and inferior longitudinal fasciculus (Lally et al., 2019; Tusor et al., 2012). In addition, FA in many of these regions was found to correlate with developmental scores of children with and without CP (Tusor et al., 2012). Our findings in a select group of school-age children who did not develop CP, who had developmental scores in the normal range at 18 months and who were attending mainstream school, demonstrated reduced FA in many of the same areas of white matter as those highlighted in neonates, providing evidence that these microstructural differences persist to an older age group, even in the absence of CP. This suggests that children cooled for NE have an altered neurodevelopmental trajectory.
The question remains whether these alterations are caused by the cooling treatment, or if there is residual damage from the initial injury resulting from NE. There is conflicting evidence regarding the impact of TH on subcortical white matter. While one experimental study reported no adverse effect of hypothermia on subcortical white matter, brain maturation or neuronal death markers (Gressens et al., 2008), other studies have suggested that TH causes cell death in subcortical white matter (O’Brien et al., 2019; Wang et al., 2016). However, the damage resulting from NE without TH (Gao et al., 2012; Gosar et al., 2020; Ly et al., 2015; van Kooij et al., 2010, 2008), and the reduction in white matter lesions with TH compared to standard care following NE (Cheong et al., 2012; Rutherford et al., 2010) suggest that these microstructural alterations are likely attributable to the hypoxic-ischemic insult that preceded NE.
4.2 Structural Connectivity Correlates with Cognitive Outcome in Cases Only
We found no significant group mean differences in the network metrics. However, when considering cognitive performance, a close relationship was revealed between structural connectivity and functional outcome in cases only. In controls, though each metric exhibited the same general trend as in cases, none of the correlations with cognitive performance were significant, indicating that individual differences in structural connectivity play a bigger role in determining FSIQ in cases than in controls. The fact that this trend emerged in relation to cognitive performance, despite finding no significant group differences in network metrics, suggests that cases exhibit a broad spectrum of connectivity impairments, ranging from mild to severe, which relate to functional outcome in these children.
In cases, the positive correlation of global efficiency with FSIQ and negative correlation of characteristic path length with FSIQ indicate a relationship between cognitive performance and network integration, which reflects the brain’s ability to carry out distributed processing (Bullmore and Sporns, 2009; Rubinov and Sporns, 2010). Also in cases, the positive correlation of local efficiency and clustering coefficient with FSIQ demonstrate a relationship between cognitive performance and network segregation, which reflects localised processing capabilities (Rubinov and Sporns, 2010). These relationships provide further evidence for the link between the severity of connectivity impairment and cognitive outcome following the brain injury.
During development, increasing network segregation is thought to be associated with pruning, while increasing strength and integration are thought to be associated with myelination (Dennis and Thompson, 2013b; Tymofiyeva et al., 2014). We found an association between reduced cognitive performance and measures of segregation and integration, reinforcing the hypothesis that the developmental trajectory of the TH children is altered, potentially impacting the processes of myelination and pruning and resulting in a ceiling effect on functional outcome at school age.
Despite the association between FSIQ and network strength, efficiency, clustering and characteristic path length, no relationship was found with small-worldness or modularity. This suggests that brain reorganisation during development prioritises small-world, modular characteristics, such that no relationship emerges between these properties and the level of cognitive impairment resulting from NE. Similar findings have been reported in school-age children born extremely preterm or with intrauterine growth restriction (Fischi-Gomez et al., 2016).
4.3 Regions Involved in Attention and Visuo-spatial Processing have Impaired Connectivity in Cases
Connectivity, measured by FA, was significantly reduced in cases compared to controls in a subnetwork comprising several sensorimotor areas including the thalamus, putamen, precentral gyrus, postcentral gyrus, paracentral gyrus and the superior parietal gyrus. The superior parietal gyrus is concerned with aspects of attention and visuo-spatial perception, including the representation and manipulation of objects. The precuneus, which appears bilaterally as two of the three most well-connected nodes in the subnetwork, is associated with numerous highly integrated tasks including visuo-spatial imagery and episodic memory retrieval (Cavanna, 2007; Cavanna and Trimble, 2006). Other nodes in the subnetwork include the insula (sensorimotor as well as higher-level cognitive function (Uddin et al., 2017)), isthmus of the cingulate cortex (which has a role in memory), inferior temporal gyrus (visual processing and visual object recognition), superior temporal gyrus (visual information integration (Karnath, 2001; Shen et al., 2017)), fusiform gyrus (object and face recognition (Kleinhans et al., 2008; Pelphrey et al., 2007)), amygdala (emotional behaviour) and the posterior cingulate cortex (internally directed thought (Leech et al., 2011) and task management (Pearson et al., 2011)). The posterior cingulate cortex is also involved in controlling attention via interaction with the cognitive control network and has been linked to attentional impairments in brain injury, autism, attention deficit hyperactivity disorder and schizophrenia (Leech et al., 2011; Leech and Sharp, 2014). Both the precuneus and the posterior cingulate cortex feature in the default mode network (Raichle et al., 2001), suggesting a role in the neural correlates of consciousness (Cavanna, 2007).
The reduced connectivity to numerous regions involved in visuo-spatial processing and attention aligns with behavioural findings from a study by Tonks et al., demonstrating reduced visuo-spatial processing, attention difficulties and slower reaction times in this group of children (Tonks et al., 2019). Similarly, the sensorimotor regions included in the network (in particular the numerous thalamocortical connections) may account for the reduced motor performance in the absence of CP (Jary et al., 2019; Lee-Kelland et al., 2020) while the impaired connectivity to the amygdala may be linked to the increased likelihood of emotional behavioural difficulties (Lee-Kelland et al., 2020).
4.4 Connectivity to Regions Involved in Visuo-spatial Processing Correlates with Cognitive Outcome
Subnetworks were found in which there is a stronger dependence of aspects of cognitive outcome (FSIQ and processing speed) on connectivity in cases than in controls. Processing speed aims to measure the mental speed and cognitive flexibility of the child; however, the score is also affected by other cognitive factors such as visuo-motor coordination, visual discrimination, attention, short-term visual memory and concentration. FSIQ is a measure of the overall cognitive ability of an individual based on performance on all WISC-IV subtests (Kaufman et al., 2006). There were no edges common to the two subnetworks, indicating that the correlation with FSIQ was not driven by correlation with processing speed.
The most well-connected nodes in the FSIQ subnetwork are involved in visuo-spatial processing, memory and attention, but there are also connections to several association cortices and visual processing areas. All connections in the FSIQ subnetwork are interhemispheric, suggesting involvement of the corpus callosum. The processing speed subnetwork consists of predominantly visual processing regions, as well as areas involved in visuo-spatial function and attention, and sensorimotor areas. Importantly, the relationship between connectivity and outcome is significantly stronger in cases than in controls, as demonstrated in Figures 6B and C. This provides an extension to the idea of ceiling effects being imposed on the cognitive processing abilities of cases, whereby the connections in the subnetwork restrict cognitive outcome in cases, whereas the cognitive processing abilities of controls are less dependent on the strength of these particular connections.
Though perceptual reasoning, verbal comprehension and working memory were reduced in cases (see Table 1), group differences in the dependence of these domains on connectivity was not found. This could be due to the dependence of these domains on connectivity being equal across subjects regardless of case-control status, or due to these domains being dependent on different connections in each subject rather than on any distinct subnetwork. This may also be dependent on how well each WISC-IV domain reflects fundamental cognitive processes versus higher-level thinking.
4.5 Major Hubs in the Human Connectome are Among Those Affected in Cases
Several studies have investigated structural brain network properties to determine key, densely connected hub nodes which constitute a structural core, or “rich club”, of the human connectome (Gong et al., 2009; Hagmann et al., 2008; van den Heuvel and Sporns, 2011). These hub nodes are thought to play a central role in information integration.
These studies consistently identified the precuneus cortex as a key node in the rich club, as well as highlighting the posterior cingulate cortex, superior parietal cortex, paracentral lobule, isthmus of the cingulate cortex, superior temporal cortex and thalamus. Additionally, sensorimotor areas were among those found to be hubs during the neonatal period (Fransson et al., 2011; van den Heuvel et al., 2015) and have been shown to be affected in dyskinetic cerebral palsy (Ballester-Plané et al., 2017), which can also result from hypoxia at birth. Many of these rich club nodes were implicated in the relationship between connectivity and case-control status, FSIQ and processing speed.
It has been suggested that, due to their topological centrality and high biological cost, rich club nodes are particularly vulnerable to a wide range of pathogenic factors (Crossley et al., 2014; van den Heuvel and Sporns, 2013). The high metabolic rates of the precuneus cortex (Cavanna and Trimble, 2006) and posterior cingulate cortex (Leech and Sharp, 2014) support this suggestion of vulnerability. Increased vulnerability may be a reason for these nodes being implicated in NE children; these nodes are affected the most by the lack of oxygen during birth therefore they sustain lasting developmental alterations.
4.6 Strengths and Limitations
To our knowledge, this is the first study to investigate whole-brain structural connectivity in school-age children treated with TH for NE, who did not develop CP. We used a robust methodology of high angular resolution DWI combined with an anatomically-constrained tractography method capable of resolving crossing fibres. Movement can be a common issue when scanning children, therefore we applied a robust quality control pipeline. The rejection of scans due to movement artefact, as well as the incomplete or unobtained scans, resulted in a relatively small sample size. However, there were no significant differences between the cognitive scores of the rejected subjects and those included in the analysis. In order to increase the robustness of the NBS results, connections were only included in the analysis if expressed in >50% of cases and >50% of controls.
5 Conclusions
We demonstrate structural connectivity deficits relating to white matter microstructure and network connectivity properties in school-age children treated with TH for NE, who did not develop CP, compared to typically developing controls. We provide evidence for a relationship between structural connectivity and cognitive outcome and further demonstrate specific brain regions and connections which are associated with case-control status and with cognitive outcome. This work demonstrates that, despite the successes of TH in reducing cases of severe disability and death, there are still aspects of brain structure which are impacted by the process of NE despite treatment with TH. Further study involving neonatal scans and longitudinal investigation of the developmental aspects of these impairments could guide follow-up care and inform future therapeutic intervention strategies.
Data Availability
The data that support the findings in this article are available upon reasonable request to the corresponding author.
Funding
This work was supported by the Baily Thomas Charitable Fund (TRUST/VC/AC/SG4681-7596), David Telling Charitable Trust, as well as Sparks (05/BTL/01 and 14/BTL/01) and the Moulton Foundation. AS is supported by the Wellcome Trust (WT220070/Z/20/Z). JB is supported by the UK Medical Research Council (MR/N026969/1). MG is supported by the EPSRC (EP/N014391/1) and by a Wellcome Trust Institutional Strategic Support Award (WT105618MA).
Acknowledgements
We thank the children and their families for participating, Ngoc Jade Thai for her assistance with MR sequences, Aileen Wilson for her radiographical expertise, and Charlotte Whitfield and Emily Broadbridge for their assistance with neuropsychological assessment..
References
