Delayed maturation of the structural brain connectome in neonates with congenital heart disease.

There is emerging evidence for delayed brain development in neonates with congenital heart disease. We hypothesize that the perioperative development of the structural brain connectome is a proxy to such delays. Therefore, we set out to quantify the alterations and longitudinal pre- to post-operative changes in the connectome in congenital heart disease neonates relative to healthy term newborns and assess factors contributing to disturbed perioperative network development. In this prospective cohort study, 114 term neonates with congenital heart disease underwent cardiac surgery at the University Children's Hospital Zurich. Forty-six healthy term newborns were included as controls. Pre- and post-operative structural connectomes were derived from mean fractional anisotropy values of fibre pathways traced using diffusion MR tractography. Graph theory parameters calculated across a proportional cost threshold range were compared between groups by multi-threshold permutation correction adjusting for confounders. Network-based statistic was calculated for edgewise network comparison. White-matter injury volume was quantified on 3D T1-weighted images. Random coefficient mixed models with interaction terms of (i) cardiac subtype and (ii) injury volume with post-menstrual age at MRI, respectively, were built to assess modifying effects on network development. Pre- and post-operatively, at the global level, efficiency, indicative of network integration, was lower in heart disease neonates than controls. In contrast, local efficiency and transitivity, indicative of network segregation, were higher compared to controls (all P < 0.025 for one-sided t-tests). Pre-operatively, these group differences were also found across multiple widespread nodes (all P < 0.025, accounting for multiple comparison), whereas post-operatively nodal differences were not evident. At the edge-level, the majority of weaker connections in heart disease neonates compared to controls involved inter-hemispheric connections (66.7% pre-operatively; 54.5% post-operatively). A trend showing a more rapid pre- to post-operative decrease in local efficiency was found in class I cardiac sub-type (biventricular defect without aortic arch obstruction) compared to controls. In congenital heart disease neonates, larger white-matter injury volume was associated with lower strength (P = 0.0026) and global efficiency (P = 0.0097). The maturation of the structural connectome is delayed in congenital heart disease neonates, with a pattern of lower structural integration and higher segregation compared to controls. Trend-level evidence indicated that normalized post-operative cardiac physiology in class I sub-types might improve structural network topology. In contrast, the burden of white-matter injury negatively impacts network strength and integration. Further research is needed to elucidate how aberrant structural network development in congenital heart disease represents neural correlates of later neurodevelopmental impairments.

Introduction

Congenital heart disease (CHD) occurs at a prevalence of about 8.2 per 1000 live births in Europe and accounts for one-third of all congenital anomalies (van der Linde ). More than one-third of all CHD infants require surgical correction or catheter interventions during early childhood, and are considered severe cases of CHD (Ferry, 1987). CHD survivors are at an increased risk for a spectrum of neurodevelopmental sequelae and cognitive dysfunction across the lifespan (Karsdorp ; Huisenga ).

Numerous studies have described aberrant brain development in CHD neonates that might underlie the risk for adverse neurodevelopment (Claessens ; Peyvandi ). With an onset in the foetal period (Limperopoulos ), delayed brain development in CHD neonates has been found on the macro-structural, micro-structural and metabolic scale (Miller ; Licht ; Ortinau ; von Rhein ; Claessens ), with the most severe CHD sub-types being most affected (Limperopoulos ; Peyvandi ). Furthermore, it has been suggested that brain dysmaturation is the substrate for white-matter injury (WMI), which is the most frequently found brain lesion in CHD neonates (Partridge ; Dimitropoulos ; Mulkey ; Guo ).

Although a link of aberrant brain development and injury with neurodevelopmental outcome has been demonstrated (Andropoulos ; Claessens ; Meuwly ), the complex phenotype of later neurodevelopmental impairments in CHD children remains insufficiently explained. Thus, more advanced neuroimaging analyses considering the complex network of the brain beyond specific structure and function mapping are needed to better capture adverse brain development in CHD neonates (Peyvandi ).

Structural connectomics is a rapidly emerging tool, which is promising in studying the emergent properties of the global brain network organization (Rubinov and Sporns, 2010) and revealing systems-level effects of disease in the developing neonatal brain (Tymofiyeva ; Song ; Keunen ; Jakab, 2019). Driven by neurobiological processes during foetal and neonatal development, the whole-brain network transitions from a highly segregated to a more integrated organisational architecture (Cao ), allowing for evolving parallel and efficient information processing across remote brain regions (Rubinov and Sporns, 2010). As previously shown, the connectome is subject to developmental disturbances in at risk populations, such as infants with neonatal ecephalopathy or after pre-term birth (Ziv ; Pandit ; Batalle ). As such, structural connectomics may detect aberrant brain development that appears normal on routine diagnostic imaging.

Recent studies demonstrated aberrations in the architecture of the neonatal whole-brain connectome in CHD, indicating possible vulnerability of brain network organization towards developmental disturbances in this population (De Asis-Cruz ; Schmithorst ). However, to date, the characteristics of the neonatal structural connectome and the development across the perioperative period in CHD neonates in comparison to healthy controls are not well studied. Furthermore, while there is evidence that WMI can disrupt the structural and functional connectome in pre-term born infants (Ceschin ; Cai ), the reason how the specific emergent properties of the whole-brain network are affected by WMI in neonates with CHD remains poorly understood.

Therefore, the primary aim of this study is to characterize the neonatal structural brain connectome in pre- and post-operative neonates with CHD in comparison to healthy controls. Second, we aim to assess risk factors for adverse network development, such as CHD severity sub-type and WMI. We hypothesize that (i) structural whole-brain connectomics will reveal systems-level aberrant brain development beyond what is evident on routine diagnostic imaging; (ii) there is an effect of CHD subgroups on pre- to post-operative network development and (iii) WMI is a risk factor for the perturbations of the neonatal connectome.

Materials and methods

Study population

In this prospective cohort study, neonates with CHD who required neonatal corrective or palliative cardiac surgery during the first weeks of life at the University Children’s Hospital Zurich between December 2009 and March 2019 were eligible. Neonatal corrective or palliative cardiac surgery was considered irrespective of the use of cardiopulmonary bypass surgery. Furthermore, cases with single ventricle physiology who underwent a hybrid approach instead of Norwood I procedure were also included. Exclusion criteria were a suspected or confirmed genetic syndrome or a gestational age below 36 weeks. A comparison group of healthy term born neonates was recruited between January 2011 and April 2019 from the post-natal ward at the University Hospital Zurich. Inclusion criteria for healthy controls were birth >36 weeks of gestation and unremarkable post-natal adaptation. None of the healthy controls were admitted to the neonatal intensive care unit.

Neonatal, perioperative, surgical and demographic characteristics were prospectively extracted from patients’ charts. CHD diagnoses were classified into anatomical sub-classes according to Clancy et al.: class I, biventricular CHD without aortic arch obstruction; class II, biventricular CHD with aortic arch obstruction; class III, univentricular CHD without aortic arch obstruction and class IV, univentricular CHD with aortic arch obstruction (Clancy ). Parental written informed consent was obtained, and the study was approved by the ethical committee of the Kanton Zurich, Switzerland (KEK StV-23/619/04). The study was carried out in accordance with the principles enunciated in the Declaration of Helsinki and the guidelines of Good Clinical Practice. In this sample, all participants enrolled in the recruitment period, who met the inclusion criteria and had good-quality neonatal diffusion tensor-imaging data were included (for more information, see flow chart in Supplementary Fig. 1). Neuroimaging and neurodevelopmental outcomes in subsamples of this cohort were published previously (Bertholdt ; Claessens ; Jakab ; Meuwly ).

Brain MRI

Both CHD neonates and healthy controls underwent brain MRI on a 3.0T clinical MRI scanner using an 8-channel head coil (GE Signa MR750). CHD neonates were scanned pre- and post-operatively, with serial MRI at both time points if possible. Healthy controls received one post-natal MRI. MRI was acquired during natural sleep with a feed and wrap technique; hearing protection was provided. The anatomical scanning protocol included an axial T1-weighted spin-echo sequence [TR (repetition time)/TE (echo time) 680/21 ms] with a slice thickness of 2.5 mm; a T2-weighted fast-spin-echo sequence in axial, coronal and sagittal planes with a slice thickness of 2.5 mm and a slice gap of 0.2 mm (TR/TE 3910/114 ms); either a 3D T1-weighted fast-spoiled gradient echo sequence (TR = 11 ms, TE at ‘minimum setting’, ranging from 4.8 to 5.3 ms, automatically determined by the scanner based on the signal absorption model used) or a 3DT1 rapid gradient-echo sequence with an identical protocol (except TE: 2.9–3.1 ms), inversion time = 450 ms, flip angle = 12.1 mm3 isotropic, acquired in sagittal plane. Diffusion tensor imaging (DTI) was axially acquired using a pulsed gradient spin-echo echo-planar imaging sequence (TR/TE 3950/90.5 ms, field of view =18 cm, matrix = 128 × 128, slice thickness = 3 mm) with 35 diffusion encoding gradient directions at a b-value of 700 s/mm2 and four b = 0 images. During the study period, an upgrade was performed from the HD.xt to the MR750 platform. This was accounted for by including MRI cohort as covariate in all multivariate regression models. All scans were reviewed for abnormal MRI findings by a neuroradiologist (R.K.).

WMI volume quantification

WMI identified as areas of T1 hyper-intensity were manually delineated on 3D T1-weighted images using Display software (http://www.bic.mni.mcgill.ca/software/Display/Display.html, 5 January 2020, date last accessed) as previously described in detail (Guo ). As applied by others (Peyvandi ), cumulative WMI volume in mm3 was defined as the maximum lesion volume on either pre- or post-operative MRI (i.e. the larger lesion volume in cases where both MRI were affected).

Diffusion tensor image processing

Details of the diffusion tensor image processing can be found in the Supplementary Methods.

Anatomical parcellation of regions of interest

Anatomical regions of interest (ROIs) were parcellated in the anatomical space of the UNC infant brain atlas for neonates (Shi ), in which a labelling system corresponding to the automated anatomical labelling (Tzourio-Mazoyer ) was used. In this study, we used 90 cortical and subcortical areas. Next, a custom FA template was created in the UNC space based on the FA images of the control subjects to transfer the UNC-automated anatomical labelling ROIs to the DTI space of each subject. The template generation and non-linear co-registration process has been described in detail previously (Jakab ).

Structural connectivity and network construction

Seeds for probabilistic diffusion tractography were defined in native DTI space by keeping the voxels with higher than 50th percentile FA value, determined from the histogram of FA values of the brain parenchyma. This outlined the white matter and the most anisotropic voxels belonging to the deep grey-matter and brainstem structures. Tractography was carried out using a probabilistic Runge–Kutta 4th-order tractography method implemented in Camino (Parker ), which used the dyadic vector output of the BedpostX algorithm of FSL. Per voxel, 100 probabilistic samples were tracked, probabilistic nearest-neighbour interpolation was used. A termination criterion of FA < 0.1 was used for stopping tractography. Tracts <15 mm were ignored to reduce the effects of spurious connections between neighbouring cortical ROIs. Undirected structural connectivity networks were generated, where each edge represented the number of probabilistic tracts passing through the white matter and connecting any two ROIs (conmat command in Camino). Weighted, undirected connectivity networks were formed. The connectivity or edge weight was defined as the mean fractional anisotropy (FAmean) value of all tractography streamlines connecting the regions at the two end-points, as these networks have been shown to be most sensitive to microstructural dysmaturation of the white matter in neonates with CHD (Schmithorst ). Nodes corresponded to the automated anatomical labelling ROIs in subject space.

Graph theoretical analysis

Network generation and analysis was performed in R Version 3.6.0 (R Core Team, 2019) with the packages igraph Version 1.2.4.2 (Csardi and Nepusz, 2006) and brainGraph Version 2.7.3 (Watson, 2019). Raw FAmean networks were thresholded by means of an iterative proportional cost thresholding procedure to filter spurious connections and eliminate the bias of cost for the comparison of network topology (van den Heuvel ). To ensure that thresholding resulted in the same absolute number of connections across the majority of subjects, we determined the 5th centile of the distribution of network densities. This T95 served as upper level of the proportional thresholding procedure, thus ensuring that 95% of the participants had the same amount of connections at each threshold. Hence, the bias of cost was reduced, while allowing for the inclusion of very sparse networks. Networks were consecutively thresholded over the cost range from 1% to T95. Graph theoretical parameters were calculated for all networks at each threshold.

Graph theory parameters

To quantify the topological organization of brain networks, and specifically depict the main characteristics of network reorganization during rapid neonatal brain development, parameters of structural integration and segregation were measured. Calculations were performed as implemented in the igraph and brainGraph package in R and as described by Rubinov and Sporns (2010). As parameter of structural integration, global/nodal efficiency was calculated. As measure of structural segregation, local efficiency and transitivity were assessed. Furthermore, network strength was quantified. A detailed description of the calculated graph parameters is given in the Supplementary Methods.

Multi-threshold permutation correction

Multi-threshold permutation correction (MTPC) is a comprehensive and sensitive approach for the statistical comparison of network topologies. It takes the trajectory of graph metrics over the whole network cost threshold range into account and thus overcomes the bias of arbitrary threshold selection for network comparison (Drakesmith ). Pairwise comparison of global and nodal graph parameters between pre- and post-operative CHD neonates and controls was performed by MTPC as implemented in the brainGraph package. To test group differences with MTPC, general linear models were calculated at each individual cost threshold and test statistics were noted. By permuting the group assignment in the general linear model across n iterations, the test-statistics under the null hypothesis were sampled for each threshold and the maximum of the null test-statistics across all thresholds was determined and noted for each permutation. The critical test-statistics threshold above which the results were significant was identified as the top αth percentile of the null-test statistics. For each cluster of ≥3 consecutive thresholds with test-statistics above the critical test-statistics threshold, the area under the curve (AUC), denoted as was calculated. Similarly, the mean of the super-critical AUCs of the permuted statistics under the null hypothesis was calculated (). The null-hypothesis was rejected if . The number of permutations was set to 5000 for nodal and 10 000 for global graph metric comparisons for each threshold as suggested by Watson .The α-level for statistically significant group differences was set to 0.025 as two one-sided hypothesis tests were performed to determine the direction of group differences. To control the inflation of the false-discovery rate arising from comparing nodal graph theoretical parameters across 90 ROIs, P-values were adjusted with the Benjamini–Hochberg procedure (Benjamini and Hochberg, 1995).

Network-based statistic

Network-based statistic (NBS) was performed to compare the strength of each network connection, i.e. at the edge-level, between CHD neonates and controls. For NBS, networks were thresholded at a single cost threshold at T95. In brief, in NBS a general linear model was specified at each individual network connection and all edges with test statistics above a pre-defined threshold were identified (here P = 0.0005 for two one-sided hypothesis tests to determine the direction of potential group differences). Among those supra-threshold connections, any connected components were identified. By means of permutation testing (here 5000 permutations) the P-value controlled for the family wise error rate was determined for each connected component based on its size (Zalesky ). We used the implementation of NBS in the brainGraph package in R.

Statistical analysis

Statistical analyses were performed using R Version 3.6.0 (R Core Team, 2019). Descriptive statistics for continuous variables were reported as mean and standard deviation or median and interquartile range as appropriate for sample distribution. Categorical variables were reported as proportions.

All statistical models specified within MTPC or NBS analyses included sex, post-menstrual age (PMA) at MRI and scanning cohort as nuisance variables to correct for confounding. The brain networks and nodal MTPC differences were visualized with BrainNet Viewer for Matlab R2014 (Xia ).

Random coefficient linear mixed effect models from the nlme package were built to assess longitudinal growth trajectories of network parameters among CHD sub-classes and the effects of WMI on network development. Therefore, a single-network threshold at the T95 maximum cost threshold was chosen. To account for the longitudinal structure of the data and allow for individual trajectories (Schielzeth and Forstmeier, 2009), the subject identifier and time at scan were included as random effects. All models included sex and scanning cohort as fixed effects. PMA at scan was mean centred and scaled to have a mean of 0 and SD of 1. Two main mixed model analyses were conducted. First, to assess network parameter development among CHD sub-classes and controls, an interaction term between sub-class and PMA at scan was specified. The control group served as the reference group. To test for further interactions between cardiac subclasses, the slopes of the subtypes were compared with each other using the glht command from the multcomp package in R. Second, to assess the effect of WMI on network development, an interaction term between cumulative WMI volume and PMA was specified. Non-significant interaction terms were split into non-interacting fixed effects in the final models. As each mixed model was calculated for four global graph theoretical parameters, the family-wise error rate was controlled with the Bonferroni correction. The Bonferroni-adjusted P-values are given as Padj.

To assess the robustness of global group differences between CHD neonates and healthy controls, we excluded all CHD neonates with WMI or arterial ischaemic stroke in a post-hoc sensitivity analysis and repeated the MTPC comparison of pre- and post-operative global graph theory parameters.

Data availability statement

The de-identified data that support the findings of this study will be made available upon reasonable request from the corresponding author.

Results

In total, 122 CHD neonates were included in the study. Among those five had to be excluded due to clinical reasons (severe brain injury after meningitis n = 1; no MRI obtained due to clinical instability n = 1; post-ponement of surgery beyond first year of life n = 1; later diagnosis of genetic syndrome n = 2). Among those remaining 117 CHD neonates, 114 had DTI sequences with sufficient quality at least at one time point, which allowed for the reconstruction of the structural connectome and the inclusion in further analyses. For further details on participants excluded for clinical reasons or poor imaging quality, we can refer to the flow chart shown in Supplementary Fig. 1. Ultimately, 79 pre-operative and 91 post-operative connectomes (56 cases with both pre- and post-operative connectomes) were included in this analysis. Of 48 included healthy controls, structural connectomes could be generated in 46 (excluded for motion n = 2). Cardiac diagnoses were classified as class I in 66/114 (57.9%), class II in 24/114 (21.1%), class III in 7/114 (6.14%) and class IV in 17/114 (14.9%). Demographic and clinical characteristics of the study cohort are summarized in Table 1.

Table 1

Demographic and clinical characteristics of study population

	CHD	Controls	P-value
N	114	46
Gestational age (weeks) mean (SD)	39.4 (1.3)	39.5 (1.3)	0.75
Birth weight (g) median [IQR]	3300.0 [3000.0, 3670.0]	3340.0 [3050.0, 3650.0]	0.67
Head circumference (cm) median [IQR]	34.5 [34.0, 35.2]	35.0 [34.0, 36.0]	0.057
Male, n (%)	83 (72.8)	21 (45.7)	0.002
Apgar score at 5 min, median [IQR]	9.0 [8.0, 9.0]	9.0 [9.0, 9.0]	0.22
Mechanical ventilation (days), median [IQR]	3.0 [2.0, 4.0]	NA
Intensive care unit stay (days), median [IQR]	6.0 [4.0, 8.0]	NA
Univentricular CHD, n (%)	24 (21.1)	NA
Cyanotic CHD, n (%)	14 (12.3)	NA
CPB surgery, n (%)	91 (79.8)	NA
Age at surgery (days), median [IQR]	10.0 [8.0, 13.8]	NA
Age at MRI (days), median [IQR]		21.0 [16.0, 27.2]
Pre-operatively	7.0 [6.0, 10.0]
Post-operatively	26.0 [21.0, 34.8]
Post-menstrual age at MRI (weeks), median [IQR]		42.3 [41.2, 43.9]
Pre-operatively	40.3 [39.4, 41.6]
Post-operatively	43.4 [42.0, 44.4]
Days pre- to post-operative MRI, median [IQR]	19 [14.0, 24.0]	NA
Days surgery/intervention to MRI, median [IQR]	14 [10.75, 19.25]	NA

Abbreviations: CHD, congenital heart disease; PMA, post-menstrual age; CPB, cardiopulmonary bypass.

Brain injury and white-matter injury burden

Arterial ischaemic stroke was observed in 4/79 (5.1%) CHD neonates pre-operatively and 4/91 (4.4%) post-operatively. In patients with serial pre- and post-operative MRI, none of the strokes was new on post-operative MRI. All were focal cortical strokes except one case with involvement of the main branch of the middle cerebral artery. WMI was detected in 13/79 (16.5%) CHD neonates on pre-operative MRI. Median [IQR] pre-operative WMI volume was 57.72 mm3 [34.41, 518.59] (missing quantification n = 1 due to missing 3D T1). On post-operative MRI WMI was evident in 12/91 (13.2%) CHD neonates with a median volume of 26.30 mm3 [16.41, 48.79]. In 1/7 (14.3%) neonates with post-operative WMI and MRI at both scanning time points, WMI was new. No brain abnormalities on MRI were found in healthy controls.

Less integrated and more segregated structural brain networks in CHD neonates compared to controls

The mean network density (cost) of the non-thresholded connectivity matrices was 0.29 ± 0.06 (range, 0.13–0.41) and the 5th centile of the cost distribution was 0.19. Thus, 0.19 was defined as T95 and network analyses were performed across the range from 1 to 19% of network cost.

Global network topology

Pre-operatively, the MTPC-based comparison of global network topology between CHD neonates and controls revealed significantly higher global efficiency in controls compared to CHD neonates. The opposite was found for local efficiency and transitivity which were higher in CHD neonates compared to controls. Post-operatively, the same differences were found: global efficiency was higher, whereas local efficiency and transitivity were lower in controls compared to post-operative CHD neonates (Table 2 and Fig. 1).

Table 2

MTPC results of global network comparison between pre- and post-operative CHD neonates and controls

Parameter	Contrast	Threshold*	β*	SE*	95% CI*	P-values*	Amtpc	Acrit
Pre-operatively
Global efficiency	Controls > CHD	0.14	0.0099	0.003	[0.0032; 0.017]	<0.001	0.26	0.16
Local efficiency	CHD > controls	0.19	0.26	0.077	[0.084; 0.44]	<0.001	0.176	0.175
Transitivity	CHD > controls	0.15	0.044	0.0095	[0.022; 0.065]	<0.001	0.42	0.18
Post-operatively
Global efficiency	Control > CHD	0.19	0.0087	0.0027	[0.0025; 0.015]	<0.001	0.25	0.14
Local efficiency	CHD > control	0.19	0.19	0.065	[0.044; 0.34]	0.0019	0.44	0.19
Transitivity	CHD > control	0.16	0.034	0.01	[0.011; 0.057]	<0.001	0.23	0.17

Significant global level differences between pre- and post-operative CHD neonates and healthy controls as revealed by MTPC. As two one-sided tests were performed to determine the direction of the effects, results are grouped by the contrast ‘controls > CHD’ or ‘CHD > controls’ and are significant at the α-level of 0.025. Amtpc and Acrit denote the results of the overall MTPC comparison across the whole range of thresholds. *Threshold indicates the cost threshold at which the strongest β* coefficient was observed. Statistical parameters are given for that threshold.

Figure 1

(A) Global graph theory parameters in pre- and post-operative CHD neonates and controls across network cost threshold range. Individual observations at each threshold are plotted as scattered points in background. (B) Results from MTPC comparison of pre-operative CHD neonates and (C) of post-operative CHD neonates versus healthy controls. Comparison corrected for age at scan, MRI cohort and sex. Green dots correspond to maximum permuted test statistics across threshold (Acrit), dashed horizontal black line depicts significant test statistics threshold defined as the top αth percentile of the null statistics distribution (α = 0.025 for two one-sided hypothesis tests). Area in magenta indicates threshold cluster at which test statistics of observed group differences were above the critical test statistics threshold (Amtpc). Null hypotheses were rejected if Amtpc > Acrit. Here for global efficiency, transitivity and local efficiency. CHD, congenital heart disease; MTPC, multi-threshold permutation correction.

Nodal network topology

Pre-operatively, similar differences in nodal level network parameters were found in the comparison of CHD neonates with controls. MTPC revealed higher nodal efficiency and strength in controls compared to pre-operative CHD neonates, whereas local efficiency and transitivity were higher in CHD neonates compared to controls (Fig. 2). Differences in nodal efficiency and local efficiency involved the frontal, limbic, parietal, occipital and insular lobe and nodes in the sub-cortical grey matter bilaterally. Higher nodal transitivity in pre-operative CHD neonates compared to controls was found predominantly in the limbic lobe in both hemispheres, but also affected parietal, frontal and sub-cortical regions. Higher strength in controls compared to pre-operative CHD neonates was found in multiple nodes in the frontal lobe, in the right pre-central gyrus and in the right inferior occipital gyrus. In contrast, significantly higher strength was found in the thalamus bilaterally in pre-operative CHD neonates compared to controls. Supplementary Table 1 provides a complete list of nodal graph theory parameter differences. Post-operatively, when comparing CHD neonates with controls no statistically significant differences in network topology were found at the nodal level.

Figure 2

Significant nodal network parameter differences among pre-operative CHD neonates and controls tested with MTPC. Comparison corrected for age at scan, MRI cohort and sex. P-values were adjusted for multiple comparison across all 90 anatomical regions of interest with the Benjamini–Hochberg procedure. The left hemisphere is displayed on the left-hand side of the image. Nodal size corresponds to the Amtpc value. Nodes are coloured according to lobe membership (pink: frontal; orange: insula; dark blue: limbic; light blue: occipital; green: sub-cortical grey matter; yellow: parietal). Supplementary Table 2 provides a list of node label abbreviations. CHD, congenital heart disease; MPTC, multi-threshold permutation correction.

Edge-level network comparison

NBS carried out at a network cost of 0.19 revealed multiple connected components with significantly different connectivity strength in CHD neonates compared to controls. Pre-operatively, two connected sub-networks with 18 nodes sharing 17 edges (P = 0.001), and 8 nodes sharing 7 edges (P = 0.005) with significantly stronger connections in controls compared to pre-operative CHD neonates were found. Of these edges, 66.7% were inter-hemispheric connections. Furthermore, three connected sub-networks with significantly stronger edges in CHD neonates compared to controls were found (16 nodes, 18 edges, P = 0.001; 11 nodes, 11 edges, P = 0.006; 8 nodes, 7 edges, P = 0.009) (Fig. 3A and B). Of these edges, 8.33% were inter-hemispheric connections. Post-operatively, two connected sub-networks were found with significantly higher connectivity strength in controls compared to CHD neonates (9 nodes, 8 edges, P = 0.003; 4 nodes, 3 edges, P = 0.033) (Fig. 3C). Of these edges, 54.5% were inter-hemispheric. Similarly, two sub-networks with stronger connections in CHD compared to controls were found (4 nodes, 3 edges, P = 0.012; 3 nodes, 2 edges, P = 0.04) (Fig. 3D) of which no edges were inter-hemispheric.

Figure 3

Results of network-based statistic showing edgewise network differences between pre- (A,B) and post-operative CHD (C,D) neonates and healthy controls. (A) Connected components with higher connectivity strength in controls compared to pre-operative CHD neonates. (B) Components with higher connectivity strength in pre-operative CHD neonates compared to controls. (C) Network components with higher connectivity strength in controls compared to post-operative CHD neonates. (D) Components with higher connectivity strength in post-operative CHD neonates compared to controls. Network-based statistic was carried out at the cost threshold of 0.19. Comparison corrected for age at scan, MRI cohort and sex. Axial and coronal views of network components are shown. The left hemisphere is displayed on the left-hand side of the image. Nodal size corresponds to the nodal degree (number of connections), node colour corresponds to lobe membership (pink: frontal; orange: insula; dark blue: limbic; light blue: occipital; green: sub-cortical grey matter; yellow: parietal). Edge size corresponds to test statistics value. CHD, congenital heart disease.

Trend-level evidence for the improvement of structural network topology in class I CHD subgroup

In random coefficient linear mixed effects models, the development of network parameters across the first weeks of life was investigated. Global efficiency and strength were positively associated with PMA at scan (β = 0.015, 95% CI = 0.013–0.017, P < 0,0001, Padj < 0.0001; β = 0.44, 95% CI = 0.38–0.50, P < 0.0001, Padj < 0.0001), whereas the opposite relationship was found for local efficiency and transitivity (β = −0.26, 95% CI = −0.31 to −0.21, P < 0.0001, Padj < 0.0001; β = −0.018, 95% CI = −0.025 to −0.012, P < 0.0001, Padj < 0.0001) (Fig. 4A). Sex was not significantly associated with network parameters or developmental slopes. An interaction term of PMA at scan and subgroup (levels: control, cardiac class I, class II, class III and class IV) was introduced to test whether the slope of the network parameter development was modified by CHD severity sub-type. For local efficiency, a significant interaction term was found, suggesting a faster decrease in local efficiency in CHD neonates with class I sub-type compared to controls. However, this interaction term did not survive Bonferroni correction for multiple comparison (β = −0.17, 95% CI = −0.32 to −0.013, P = 0.034, Padj = 0.14) (Fig. 4B). Furthermore, there was no difference in slopes among other sub-classes compared to controls or compared to each other. No significant interaction was found for the remaining network parameters. Neither did the grouping of CHD sub-types by arch obstruction or bi- versus uni-ventricular cardiac physiology show evidence for differences between CHD subtypes.

Figure 4

Results from network development analysis with random coefficient mixed models. (A) Association of global graph theory parameters with PMA at scan including all CHD and control connectomes. (B) Trajectory of local efficiency development among CHD severity sub-types and healthy controls. Trajectories are overlaid on individual data points plotted as grey dots and connected by a thin grey line, if they represent longitudinal measurements. Differences in slopes were not statistically significant after correction for multiple comparison. Cardiac sub-types were grouped according to cardiac severity classes by Clancy . PMA, post-menstrual age; CHD, congenital heart disease.

White-matter injury is associated with lower strength and integration of global networks across CHD sub-types

In a second random coefficient model including only CHD neonates, an interaction term of cumulative WMI volume and age at scan was introduced to test the effect of WMI burden on network development. Higher cumulative WMI volume was negatively associated with global efficiency (β = −0.0036, 95% CI = −0.0062 to −0.00092, P = 0.0086, Padj = 0.034) and strength (β = −0.15, 95% CI = −0.24 to −0.055, P = 0.0021, Padj = 0.0084). However, WMI did not modify the effect of age at scan on network parameter development.

Post-hoc sensitivity analysis

To test the robustness of global graph theory differences between CHD neonates and controls, we repeated the global network comparison after excluding all participants with WMI or arterial ischaemic stroke (n = 18 pre-operative, n = 18 post-operative CHD neonates). Pre-operatively, global efficiency remained significantly lower, and transitivity remained significantly higher in CHD neonates without brain lesions compared to healthy controls. Post-operatively, no significant differences in global network topology were found between CHD neonates without brain lesions and controls.

Discussion

In this prospective cohort study, we found less mature structural whole-brain connectomes in neonates with CHD compared to healthy controls on all levels of network analysis. The maturational delay characterized by lower structural integration and higher segregation was more prominent in the pre-operative but persisted to the post-operative period. Trend-level evidence suggested that cardiac physiology might modify the trajectory of perioperative network development. Importantly, larger WMI volume was associated with reduced global efficiency and strength of the whole-brain network topology and was thus identified as a risk factor for perturbation of the neonatal network development.

Reduced integration and higher segregation of global and nodal network topology in CHD neonates

We found significantly altered structural whole-brain network topology in CHD neonates compared to healthy controls. The differences were evident pre-operatively and persisted, albeit less prominent, post-operatively. They can be summarized as a pattern of reduced integration and increased segregation. Similar findings of a perturbated network architecture in CHD neonates compared to controls were reported by Schmithorst . In their study, CHD neonates underwent either pre- or post-operative MRI. Pre- and post-operatively reduced global and nodal efficiency, i.e. structural integration were found in FAmean networks, however, when controlling for network cost, these differences did not remain significant. In contrast, our proportional cost thresholding paired with MTPC revealed network topological differences between CHD neonates and controls that were independent of network cost. Furthermore, we were able to reveal perturbations in network segregation. In contrast to our analysis, CHD neonates with any type of brain injury were excluded by Schmithorst . We assume that our results are more representative of network disturbances found in CHD neonates as perioperative brain injury, such as punctate WMI and focal stroke, are common in CHD (Peyvandi ). However, the different findings in our study and the analysis by Schmithorst cannot be entirely explained by the cases with brain injury as revealed by our post-hoc sensitivity analysis. We found that after exclusion of cases with injury pre-operative global network differences remained, whereas post-operatively significant differences did not. Beside the effects of injury, the loss of post-operative differences might have also resulted from reduced sample size and statistical power.

Opposed to our structural connectomics findings, a previous analysis of the functional connectome demonstrated a preserved functional global network architecture in CHD neonates (De Asis-Cruz ). This contrast in global topological findings could be partly attributed to the exclusion of CHD neonates with parenchymal brain lesions, but could as well arise from temporal differences in functional and structural network maturation. At the neonatal age, structural precedes functional network development, whereas the coupling of both networks gradually increases from the last trimester of gestation to 20 years of age (Cao ). Therefore, structural might precede functional connectomic disturbances, with the latter only becoming apparent at later ages.

Our network edge-level analysis revealed multiple sub-networks that differed between CHD neonates and controls pre- and post-operatively, with inter-hemispheric connections appearing most vulnerable to perturbations in CHD. Impaired inter-hemispheric connectivity evident as volumetric and microstructural dysmaturation in the corpus callosum has been reported in a small sub-sample of our CHD cohort (Hagmann ). Our findings also accord with a study, showing that ∼60% of functional connections that were weaker in CHD neonates than in controls were inter-hemispheric (De Asis-Cruz ).

Connectome topology reflects structural network dysmaturation in CHD neonates

The longitudinal trajectory of network development in our cohort reflects that global efficiency and strength increased significantly across the investigated period of neonatal development, whereas parameters of segregation such as local efficiency and transitivity significantly decreased. This is in line with the suggested normative trajectory of network development during foetal and neonatal brain development, characterized by an increase in network integration (i.e. global efficiency) and decrease in network segregation (i.e. local efficiency and transitivity) (Takahashi ; Tymofiyeva ; Jakab ; Song ; Keunen ; Turk ). Axonal growth, synaptogenesis and myelination starting during the late second and early third trimester of gestation lead to the emergence of inter-hemispheric connections and intra-hemispheric long-range association fibres (Keunen ) and contribute to an increasing structural integration of the whole-brain network (Tymofiyeva ; Song ). This increasing efficiency is paralleled by decreasing segregation as result of synaptic pruning and the refinement of the specialized sub-networks (Damaraju ). Given this pattern of structural connectome development, our finding of reduced global efficiency, and higher local efficiency and transitivity in CHD neonates indicates a maturational delay in network architecture. These observations in the early pre-operative period suggest that network dysmaturation in CHD neonates might be of foetal origin with a persistence into the post-operative period. Furthermore, our findings demonstrate an impairment beyond the previously described multi-faceted brain developmental delay in CHD neonates that involves structural and metabolic aberrations and originates in the last trimester of gestation (Limperopoulos ; Berman ; Clouchoux ).

Cardiac physiology might modify perioperative network development in CHD neonates

The trajectories of brain development in CHD neonates can be modified by the type of cardiac defect and post-operative cardiac physiology as shown previously (Peyvandi ). When comparing trajectories between CHD severity sub-groups and controls, we found trend-level evidence for a difference in local efficiency development between class I CHD neonates and controls. More rapid decrease of local efficiency from the pre- to the post-operative period could suggest a catch-up in network development in this CHD sub-class on the mild end of the severity spectrum. The majority of class I CHD neonates had a transposition of the great arteries, and thus received surgical restoration of normal cardiac physiology. This improvement in network topology from the pre- to the post-operative period evident in the largest CHD subgroup could drive our observation that the difference in network topology between CHD neonates and controls was less prominent post-operatively. However, when conservatively correcting for multiple comparison, this difference did not remain significant. While this is consistent with similar white-matter microstructural change rates in neonates with different sub-types of CHD (Claessens ), further research with even larger and more balanced subgroups is needed to explore the modifying effect of cardiac physiology on network development.

White-matter injury is associated with widespread network perturbations

We found evidence for a negative effect of WMI burden on global efficiency and strength of the structural whole brain network in CHD neonates, identifying WMI as potential risk factor for network perturbations. This disruptive effect of WMI on connectivity is in line with functional connectomic studies in pre-term born infants, revealing a negative association between WMI and thalamocortical (Duerden ) as well as inter-hemispheric connectivity (Smyser ).

Given the differences in spatial distribution of WMI in CHD compared to pre-term infants, with a paucity of central WMI in CHD (Guo ), it seems striking that the mainly small punctate lesions have such widespread effects on global network integration and strength in CHD neonates. Although the impact of lesion location on network dysconnectivity warrants further investigation, this could highlight that the overt WMI depicted on diagnostic MRI is accompanied by covert and widespread microstructural white-matter dysmaturation (Dimitropoulos ). The latter is only evident on advanced neuroimaging and could explain the systems-level effects of WMI observed in our study. There is evidence that WMI furthermore impacts longitudinal micro- and macro-structural brain development in CHD infants with moderate to severe WMI or stroke when compared to infants without ischaemic injury (Claessens ). Moreover, brain growth was found to be impaired in CHD neonates with moderate to severe WMI (Peyvandi ). In our cohort, covering a narrow range of age at scan, we did not find a modifying effect of WMI on the slope of network development.

Similar findings in other age groups and potential functional relevance

In line with our findings, a previous study found that trend level reduced global efficiency and significantly higher modularity, indicative of reduced structural integration and increased segregation, in adolescents with transposition of the great arteries (Panigrahy ). Together, our findings suggest that neonatal network perturbations might persist into adolescence; however, this has yet to be tested in a longitudinal data set. Furthermore, Panigrahy et demonstrated that disrupted network organization in CHD adolescents is of functional relevance as network topology was found to mediate cognitive performance across multiple domains (Panigrahy ). However, the reason how the altered neonatal network architecture in CHD infants impacts later neurodevelopment remains to be elucidated.

Limitations

The following limitations of our study need to be mentioned. Our analysis is based on DTI tractography, which is inevitably an imperfect reconstruction of complex WMI fibre tract architecture (Tymofiyeva ). The neonatal brain is still undergoing myelination and many of the white-matter structures underlying cortico-cortical connectivity are unmyelinated, which leads to lower diffusion anisotropy than in children or adults. Lower anisotropy results in higher uncertainty when estimating fibre orientations, which is particularly challenging in crossing-fibre regions. While our study protocol represents a realistic compromise between image quality and scan time in this sensitive patient population, studies using higher angular resolution data sets will likely provide better estimates of whole-brain structural connectivity.

As only two imaging time points were acquired in our cohort, the exploration and modelling of more complex growth trajectories were hindered. However, other aspects of brain development such as structural brain growth and functional network development do not follow linear but rather asymptotic growth trajectories (Dosenbach ; Holland ). This might be particularly true during the first months of life, as structural brain growth shows a tapering off of the daily growth rate at about 3 months (Holland ). Thus, we expect that the development of network parameters during this time might follow more complex growth trajectories too. More densely sampled longitudinal neuroimaging data across a longer timespan would be necessary to establish and investigate more complex growth trajectories of structural network architecture in healthy and aberrant neurodevelopment. However, in a typical clinical setting, this remains a challenging goal particularly for neonates with critical CHD.

In a small fraction of included CHD neonates, in our study, arterial ischaemic strokes were observed. These strokes were mainly focal and restricted to the cortex. As the group of neonates with stroke was very small, no further testing of the effect of stroke on the structural connectome was performed. This remains to be addressed in future research.

A further limitation of the interpretation of our findings is that results might depend on the connectomics approach deployed. For instance, the parcellation scheme, the tractography algorithm or the way cost-thresholding is applied. It remains a challenge of this rapidly emerging field to establish methodological standards and ensure replicability and comparability of results (Tymofiyeva ). Further reproducibility studies are needed to help depict and control technical bias.

Conclusion

In our cohort of neonates with CHD undergoing cardiopulmonary bypass surgery, we found evidence for a maturational delay in the structural brain connectome characterized by lower integration and higher segregation. This network dysmaturation was most prominent in the pre-operative but persisted to the post-operative period. We found trend-level evidence for a modifying effect of CHD severity sub-type on network development, suggesting that post-operative normalization of cardiac physiology may improve structural network topology. WMI burden was associated with perturbation of the structural connectome, impacting global network strength and efficiency. Our findings advance our understanding of aberrant brain development in CHD neonates. We revealed potential risk factors for aberrant network development, and mechanisms how these findings could represent neural correlates of later neurodevelopmental sequelae.

Supplemental material

Supplementary material is available at Brain Communications online.

Click here for additional data file.