Structural-covariance networks identify topology-based cortical-thickness changes in children with persistent executive function impairments after traumatic brain injury.

Paediatric traumatic brain injury (pTBI) results in inconsistent changes to regional morphometry of the brain across studies. Structural-covariance networks represent the degree to which the morphology (typically cortical-thickness) of cortical-regions co-varies with other regions, driven by both biological and developmental factors. Understanding how heterogeneous regional changes may influence wider cortical network organization may more appropriately capture prognostic information in terms of long term outcome following a pTBI. The current study aimed to investigate the relationships between cortical organisation as measured by structural-covariance, and long-term cognitive impairment following pTBI. T1-weighted magnetic resonance imaging (MRI) from n = 83 pTBI patients and 33 typically developing controls underwent 3D-tissue segmentation using Freesurfer to estimate cortical-thickness across 68 cortical ROIs. Structural-covariance between regions was estimated using Pearson's correlations between cortical-thickness measures across 68 regions-of-interest (ROIs), generating a group-level 68 × 68 adjacency matrix for patients and controls. We grouped a subset of patients who underwent executive function testing at 2-years post-injury using a neuropsychological impairment (NPI) rule, defining impaired- and non-impaired subgroups. Despite finding no significant reductions in regional cortical-thickness between the control and pTBI groups, we found specific reductions in graph-level strength of the structural covariance graph only between controls and the pTBI group with executive function (EF) impairment. Node-level differences in strength for this group were primarily found in frontal regions. We also investigated whether the top n nodes in terms of effect-size of cortical-thickness reductions were nodes that had significantly greater strength in the typically developing brain than n randomly selected regions. We found that acute cortical-thickness reductions post-pTBI are loaded onto regions typically high in structural covariance. This association was found in those patients with persistent EF impairment at 2-years post-injury, but not in those for whom these abilities were spared. This study posits that the topography of post-injury cortical-thickness reductions in regions that are central to the typical structural-covariance topology of the brain, can explain which patients have poor EF at follow-up.

Introduction

Traumatic brain injury (TBI) in childhood and adolescence is a leading cause of disability (World Health Organization, 2006), and injuries occur in the individual-context of a still-developing brain (Wilde et al., 2012a). Paediatric TBI (pTBI) has a reported incidence between 1.10 and 1.85 cases per hundred for the 0–15 age range (McKinlay et al., 2008) and has specific adverse effects on neurodevelopment. TBI can result in pathology at a micro and macroscopic level, leading to both transient and permanent impairments (Bigler, 2007b, 2016; Maxwell, 2012). Damage may manifest as trauma-related, developmentally inappropriate atrophy (Bigler, 2013; Urban et al., 2017; Wilde et al., 2005) which can appear as relative changes in both brain volume (Bigler, 2016) and cortical-thickness (Urban et al., 2017) on structural magnetic resonance imaging (sMRI). In pTBI, injury-related brain abnormalities occur during a period of ongoing age- and development-dependent brain changes (Bigler, 2016; Maxwell, 2012)

Previous sMRI studies have shown that, from early to post-chronic timepoints, the morphometry of the injured brain differs from that of typically developing children (see King et al. (2019)). These cross-sectional differences are observed even up to 10 years post-injury (Beauchamp et al., 2011; Serra-Grabulosa et al., 2005) suggesting that alterations are non-transient, neither recovering nor being compensated for over time. These cross-sectional differences are evidence of a long-term effect of TBI on the morphometry of the brain.

Previous studies have suggested that, from the perspective of clinical characteristics, two traumatic brain injuries can superficially appear similar but result in vastly different outcomes (Bigler, 2007a; Schneider et al., 2014). The location and extent of focal lesions to the brain following a pTBI do not fully explain post-injury neuropsychological deficits (Bigler, 2001). There is also limited evidence of brain-behaviour relationships between brain morphometry differences and functional outcomes (King et al., 2019). The paucity of reliable relationships between brain structural measures and long-term outcomes of neuropsychological functions may be explained in part by the fact morphometric changes can be highly distributed across the cortex even within a single patient, and these changes vary across individuals (Bigler, 2007a; Bigler et al., 2013; Bigler and Maxwell, 2011). This spatial heterogeneity of damage and post-injury changes may limit the potential of univariate investigations of morphometry, which are therefore unable to tell us about the subtleties of the more diffuse effects of an injury. Therefore, looking more widely at the global effects of injury and how focal damage can change the wider ‘system’ of the brain may explain greater variance in functional outcomes post-injury. One way to explore this hypothesis is to investigate changes to the global neural-network following injury due to pTBI, in keeping with recent characterizations of TBI as a disorder of brain connectivity (Hannawi and Stevens, 2016; Hayes et al., 2016), utilising a graph-theory framework to quantitatively describe these networks. This will better capture the multifaceted nature by which the brain can experience pathological change post-injury.

Patterns of grey matter morphometry across the cortex can be interpreted as a biologically-meaningful brain network, capturing the meso‑scale organisation of brain structure across the cortex. This structural-covariance network models the degree to which the morphology (measured with cortical-thickness) of brain regions statistically co-varies across all possible pairs of ROIs (Alexander-Bloch et al., 2013a, 2013b; Evans, 2013; Lerch et al., 2006; Mechelli et al., 2005). These networks are sensitive to neurodevelopmental and age-related change (Alexander-Bloch et al., 2013; Fan et al., 2011; Khundrakpam et al., 2017; Khundrakpam et al., 2016, 2013; Raznahan et al., 2011; Váša et al., 2017), with regions showing similar/shared developmental trajectories being more similar in morphometry (Alexander-Bloch et al., 2013), likely driven by the gene-controlled patterning of cortical-thickness and structural-covariance across the cortex (Romero-Garcia et al., 2018; Yee et al., 2017).

Given this highly coordinated, genetically programmed developmental ‘blueprint’ of brain maturation, neurological disruption to the structure of the brain during this period can have a significant impact on subsequent brain development, detectable as an abnormality in structural-covariance across the cortex. In other forms of paediatric brain insult, including malformations of cortical development in neonates, structural-covariance metrics change as a function of the specific gestational-timing of disruption (Hong et al., 2017). Therefore, the structural-covariance approach may be sensitive to the effects of pTBI on the developmental trajectory of the brain, respecting the complex organisation of the GM across the whole cortex, which cannot be captured with techniques that use univariate ROIs.

Using the correlational structure of regional-level morphometry following pTBI compared to controls, Spanos et al. (2007) investigated volumetric correlations across cerebro-cerebellar regions and found a significant positive relationship between DLPFC/cerebellum in the typically developing, but not in the TBI group. Drijkoningen et al. (2017) estimated the correlational structure of ‘atrophy’ scores between regions, finding moderate to very strong positive correlations. These regional morphometric relationships lend support to the hypothesis of a diffuse pattern of pathology following TBI. The current study expands on these previous findings by investigating these relationships across the whole brain utilising structural-covariance which is novel to the field of TBI.

The topology of brain networks ordinarily makes them inherently robust to insult (Hillary and Grafman, 2017). However, targeted damage to topologically central regions may have a disproportionate impact on the network and is more likely to result in behavioural symptomatology (Crossley et al., 2014; Hillary and Grafman, 2017). In the clinical setting, neurological and psychiatric disorders emerging in childhood may also be linked to abnormal development of these ‘hub’ regions (Morgan et al., 2018). During development, there is early formation of hubs in the structural network, providing a stable scaffold to build upon during subsequent development. Later, maturational change is most prominent in hub regions, with hub locations becoming more adult-like across childhood (Csermely et al., 2013; Morgan et al., 2018; Oldham and Fornito, 2019). The protracted development of structural-covariance hubs over childhood and adolescence, with those responsible for higher integrative functions developing most slowly (Khundrakpam et al., 2013; Whitaker et al., 2016), may put them at greatest risk to pathology which may result in delayed or disrupted development (Morgan et al., 2018). Therefore, damage to regions central to the network during this period, may result in behaviourally relevant changes to the developmental trajectory of the brain.

Current research echoes this sentiment, suggesting that brain damage is specifically linked to network structure of the brain. Across multiple neurological disorders, the probability of a region showing case-control differences in grey matter morphology is significantly related to the nodal-degree (the summation of the number of connections) of the region (Crossley et al., 2014). Similarly, voxels with significant case-control differences in grey matter volume/density, specifically reductions, belonged to regions with a greater median degree than ‘undamaged’ voxels (Crossley et al., 2014). In adult TBI specifically, reductions in ‘hubness’ of nodes (betweenness and eigenvector centrality) derived from a tractography network were related to greater cognitive impairment, including executive functions (Fagerholm et al., 2015). These results highlight the fact that the behavioural consequences of brain damage are highly dependent on the topological position of the damage within, and the organization of, neural networks (Aerts et al., 2016; Hillary and Grafman, 2017). It is therefore reasonable to hypothesise this may be the case in pTBI.

Aims and hypotheses

The aims of the current study were twofold; (1) to capture global, diffuse nature of the effects of TBI on the still-developing brain using multivariate-network methodologies, and (2) to investigate whether the cortical topography of post-pTBI cortical-thickness reductions in relation to the typical topology of the brain is related to poor cognitive functioning at two-year follow-up.

We examined three hypotheses:

Patient groups (including those subgroups with good/poor EF outcome) would show significant cortical-thickness reductions in comparison to heathy controls.

Patient groups would show differences in structural-covariance compared to controls.

We predicted that for pTBI, cortical-thickness reductions would have occurred in regions that have higher structural-covariance than randomly selected regions.

Methods

Ethics statement

Data from the TBI cohort in the current study was obtained under a material transfer agreement between the Murdoch Children's Research Institute and Aston University. This cohort were originally acquired for a study that had previously received ethical approval via the Human Research and Ethics Committee of Royal Children's Hospital, Melbourne, Australia. We also acquired additional control data through the public Autism Brain Imaging Data Exchange (ABIDE) database, as shared by the Preprocessed Connectome Project (PCP). A favourable opinion was granted by Aston University ethics panel for the secondary analysis of both the TBI and ABIDE datasets.

Participants

TBI cohort

The data used in the current experiment are a subset of an existing dataset of children who have experienced a TBI between the ages of five and 16 years of age. 157 children (patients n = 114) were recruited between 2007 and 2010 into a study on ‘Prevention and Treatment of Social Problems Following TBI in Children and Adolescents’. Further details have recently been published elsewhere (Anderson et al., 2013, 2017; Catroppa et al., 2017). In brief, children with TBI were recruited on presentation to the emergency department at the Royal Childrens’ Hospital, Melbourne, Australia. Eligibility for the study was determined if they: (i) were aged between five and 16 years at the time of injury, (ii) had recorded evidence of both a closed-head injury and also two post-concussive symptoms (such as headaches, dizziness, nausea, irritability, poor concentration), (iii) had sufficient detail within medical records (Glasgow Coma Scale (GCS; Teasdale and Jennett (1974)), neurological and radiological findings) with which to determine the severity of the injury, (iv) had no prior history of neurological or neurodevelopmental disorder, non-accidental injuries or previous TBI, and (v) were English-speaking. TD controls were also recruited and were required to meet criteria (i), (iv) and (v).

Injury severities were defined clinically using a combination of imaging findings, pathology and GCS score. These were defined as: (i) Mild TBI; those with GCS of 13–15, no neurological deficit and either no clinical CT/MRI scan was conducted or no evidence of mass lesion on clinical CT/MRI findings, (ii) Mild-complicated TBI; GCS of 13–15 but abnormal clinical CT/MRI findings indicating mass lesion, (iii) Moderate TBI; GCS of 9–12 and/or evidence of mass-lesion or evidence of specific injury on clinical CT/MRI and/or neurological deficit and, (iv) Severe TBI; GCS of 3–8 and/or evidence of mass-lesion or evidence of specific injury on clinical CT/MRI and/or neurological deficit.

We applied a number of inclusion criteria to the dataset. We only included subjects who; (a) met strict quality control criteria of Freesurfer outputs, and (b) had MRI data available and were scanned <90 days post-injury. This resulted in a subset of n = 116 subjects (TBI patients (n = 83); healthy controls (n = 33)) who underwent MRI acutely after injury (range = 1–88 days). To note, one pTBI subject included had a diagnosis of ADHD. Table 1 shows the group demographics, whilst Fig. 1 outlines how participant samples were derived.

Table 1

Demographics for patients and controls.

Group	pTBI	Controls	Comparison
N	83	33	–
M/F	54/29	20/13	OR= 0.83, p=.67a
Age at Scanning (median, yrs.)	10.92	9.99	F(1114)= 0.262, p=.61b
(Range, yrs.)	6.09–14.82	6.53–15.47	–
Age at Injury (median, yrs.)	10.92	–	–
(Range, yrs.)	5.92–14.67	–	–
Injury-Scan Interval (median, days)	34	–	–

Fisher's exact test (OR = odds-ratio).

One-Way ANOVA.

Fig. 1

Flowchart indicating how samples were derived for this secondary analysis of data from the initial recruitment into the ‘Prevention and Treatment of Social Problems Following TBI in Children and Adolescents’ study. Adapted from Consort Diagram (http://www.consort-statement.org/consort-statement/flow-diagram).

MRI acquisition

MRI images were acquired at 3T on a Siemens Trio scanner (Siemens Medical Systems, Erlangen, Germany) using a 32-channel matrix head coil. The acquisition included a sagittal three-dimensional (3D) MPRAGE [TR = 1900 ms; TE = 2.15 ms; IR prep = 900 ms; parallel imaging factor (GRAPPA) 2; flip angle 9°; BW 200 Hz/Px; 176 slices; resolution 1 × 1 × 1 mm] and sagittal 3D T2-FLAIR non-selective inversion preparation SPACE (Sampling Perfection with Application-optimised Contrast using different flip-angle Evolution) [TR = 6000 ms; TE = 405 ms; inversion time (TI) = 2100 ms; water excitation; GRAPPA Pat2; 176 slices; 1 × 1 × 1 mm resolution matched in alignment to the 3D T1-weighted sequence].

ABIDE dataset

To provide a second healthy reference group for validation of findings, we employed the open-access data from the Autism Brain Imaging Data Exchange (ABIDE, Di Martino et al. (2014)), specifically the pre-processed version of the dataset made available by the Preprocessed Connectome Project (PCP, Bellec et al. (2013); see http://preprocessed-connectomes-project.org/). The ABIDE dataset consists of a large sample of 532 individuals with autism spectrum disorders and 573 typical controls, composed of MRI (functional and structural) and phenotypic information for each subject, accumulated across 17 independent sites. The scan procedures and parameters are described in more detail elsewhere (http://fcon_1000.projects.nitrc.org/indi/abide/).

We applied similar inclusion criteria to this dataset, only including subjects who; (a) passed a strict MRI quality control criteria of raw structural MRI (see supplementary materials), (b) were recorded as controls within the ABIDE database, (c) at time of scan were aged < 17 years and (d) had pre-processed Freesurfer data available as part of the PCP release. This resulted in a final reference group of n = 327 (M/F = 259/68, median age (years.) = 12.49, age range (years.) = 6.47 – 16.93). The list of IDs for ABIDE subjects included in these analyses can be found in supplementary materials, as per ABIDE's recommendations.

Both controls in the experimental cohort and the ABIDE cohort had qualitatively similar mean IQ (M = 105.4 and M = 109.8) as measured across multiple age-appropriate IQ tests (in the experimental cohort IQ was assessed by WASI 2 scale IQ whereas the measures used by the ABIDE dataset were varied, see ABIDE documentation for details).

MRI processing

3D tissue segmentation and estimation of cortical-thickness and estimated total intracranial volume (eTIV) from T1-weighted (T1w) MR images were conducted using an established pipeline (Freesurfer version 6.0; see Fischl (2012) for review). The steps involved are documented elsewhere (Fischl et al., 2004) but briefly, T1w images were stripped of non-brain tissues (Segonne et al., 2004), GM/WM boundaries were tessellated and topology was automatically corrected (Fischl et al., 2001; Segonne et al., 2007). Finally, deformation of this surface was performed, to optimally define the pial (Cerebro-spinal fluid/GM) and white (GM/WM) surfaces using intensity gradients to estimate where intensity maximally shifts to define boundaries of these tissue classes (Dale et al., 1999; Dale and Sereno, 1993; Fischl and Dale, 2000). All patients and controls included in the analyses here (n = 83 & n = 33 respectively, including patients with and without lesions) had FLAIR images acquired which were used to supplement the Freesurfer segmentation algorithm as per the standard Freesurfer pipeline (using the –FLAIRpial command). All cases without these T2 FLAIR MRI failed Freesurfer QA and were not included in the analyses presented here.

In this study, Freesurfer was used to estimate the cortical volume/thickness for 34 regions-of-interest per hemisphere, based upon the cortical parcellation of the Desikan–Killiany atlas (Desikan et al., 2006). This parcellation was chosen over a more fine-grained parcellation scheme due to concerns over statistical power if a greater number of ROIs were analysed.

The quality of Freesurfer outputs was assessed using Qoala-T (Klapwijk et al., 2019) as a decision support tool to guide the systematic and replicable selection of which cases required manual editing. Visual inspection of surfaces after manual editing resulted in the exclusion of n = 15 pTBI cases (including those without T2 FLAIR MRI as described above) due to continued poor quality of surface reconstruction.

Multiple cases within the original TBI cohort also had frank parenchymal lesions. For these cases, Freesurfer has limited applicability with its standard processing pipeline and thus an adjusted pipeline was utilised and is described in Supplementary Materials. Eight lesion cases were retained for analysis using this pipeline.

Processing using the Freesurfer pipeline had already been done for the ABIDE dataset within the PCP, using the standard pipeline as described above (however using an older version of Freesurfer (version 5.1). Details of quality assurance of the anatomical processing using Freesurfer for the ABIDE data, and steps to control for ABIDE site effects, can be found in Supplementary materials.

Executive functions (EF)

EF was assessed for patients in the TBI cohort (patients and controls) at 24-months post injury/recruitment using performance-based neuropsychological testing. Several standard, age-appropriate neuropsychological tests were administered to participants to index EF skills, and these were from three typical, age-appropriate test batteries; (i) Tests of Everyday Attention – Children (TEA-Ch; (Manly et al., 1999)), (ii) Delis-Kaplan Executive Function System (D-KEFS, (Delis et al., 2001)), and (iii) Wechsler Intelligence Scale for Children (WISC-IV, (Wechsler, 2003)). These measures were selected from a wider battery of administered neuropsychological tests as part of the wider study. Specific subtests used in the current study were selected to represent components of a three-factor EF model (Miyake et al., 2000) and can be found in Table 2 As per Diamond (2013), we also included a measure of selective-attention in the domain of inhibition, given evidence highlighting the high correlation of this skill with other EF domain (Downing, 2015; Santa-Cruz and Rosas, 2017). It is important to note, whilst we have assigned these subtests to specific subdomains of EF based upon our theoretical model of EF, these assignments are not used in the designation of the EF impairment/non-impairment groupings.

Table 2

Neuropsychological tests and subtests used to group patients on executive functioning outcome 2 years post-injury.

EF Domain	Battery	Subtest	Measure
Set Shifting	TEA-Ch	Creature counting	Accuracy (no. correct)
	TEA-Ch	Creature counting	Time taken
Inhibition	D-KEFS	Colour-word interference – condition 3	Time Taken
	D-KEFS	Colour-word interference – condition 4	Time Taken
	TEA-Ch	Walk-don't-walk	Score
	TEA-Ch	Skysearch	Attention Score
Working Memory	WISC-IV	Digit span backwards	Score

An approach to define those individuals exhibiting clinically relevant cognitive impairment was selected (a-priori) to group patients in terms of executive (dys)function at 2 years post-injury. The current study adopted the neuropsychological impairment (NPI) rule proposed by Beauchamp et al. (2015). This rule has previously been shown to be sensitive to TBI with an increase in identification of impaired individuals (Beauchamp et al., 2015), and has been used to detect behavioural impairment (Donders and DeWit, 2017), and cognitive inefficiency (Beauchamp et al., 2018) following paediatric TBI and concussion respectively.

Briefly, performance scores for the neuropsychological test batteries were converted to age-scaled scores (M = 10, SD=3). To identify those with a clinically relevant impairment in executive functioning a cut-off of 1 SD outside ‘average’ functioning in the direction of worse performance. To be assigned to the group who were experiencing clinically relevant cognitive impairment (poor EF outcome (EFPoor)), participants had to have shown impaired functioning on two or more individual EF measures, whereas those who were impaired less than two measures were assigned to the without cognitive impairment group (good EF outcome (EFGood)). A minimum of two cases of impairment identifies a pattern of deficit, unlikely to be due to typical variability due to individual differences or measurement error for instance. We only calculated the NPI rule for those cases that had the full battery of EF tests. The demographics of these two subgroups (EFPoor and EFGood) are shown in Table 3. These groups did not differ on injury variables, such as injury severity, as can be seen in Table 3.

Table 3

Demographics for patient subgroups.

Group	Control	EFGood	EFPoor	Statistical comparison
N	33	42	17	–
M/F	20/13	27/15	12/5	p=.78a
Age at Scanning (median, yrs.)	9.99	10.95	11.13	F(2,89)= 0.366, p=.70b
(Range, yrs.)	6.53–15.47	6.69–14.82	6.09–14.17	–
Age at Injury (median, yrs.)	–	10.75	11.00	F(1,57)= 0.027, p=.87b
(Range, yrs.)	–	6.58–14.67	5.92–14.00	–
Injury-Scan Interval (median, days)	–	35.5	30.0	F(1,57)= 1.971, p=.17b
Injury Severity
Mild	–	23	10	p= .58a
Mild-Complicated	–	4	3
Moderate	–	11	4
Severe	–	4	0
Note. a Fisher's exact test, b One-Way ANOVA

It is important to note that of the control group, n = 4 (∼∼12%) had EF scores that satisfied the NPI rule. Given the NPI rule used a SD of 1.5 below mean, standard distribution would suggest that the 4 (12%) controls with poor EF is representative of expected normal variation. Thus, these controls were retained.

Statistical analysis

The analysis plan of the current study was inspired by that of Wannan et al. (2019). All analyses were conducted with a series of packages in R (R Core Team, 2016), with network analyses being specifically conducted using ‘brainGraph’ version 2.2 (Watson, 2016), which is an expansion of the iGraph package (Csardi and Nepusz, 2006). All analyses were conducted over three group-contrasts; (i) pTBI patients vs. controls, (ii) pTBI EF intact vs. controls and iii) pTBI EF impaired vs. controls. Only case-control comparisons were conducted, and thus did not include case-case contrasts (i.e. EF Impaired vs EF Intact). This is because we specifically wanted to investigate pathological deviations to the typical development of the brain.

Differences in cortical-thickness between pTBI and controls

Firstly, we investigated cross-sectional differences in cortical-thickness between patients and our experimental controls. For each ROI (n = 68) a general linear model (GLM) was generated to test the effect of group (patient vs control) on cortical-thickness, whilst controlling for the effects of age at scanning, sex, and eTIV. A t-test was used to test the directional hypothesis of cortical-thickness reductions in the patient group compared to controls. When calculating p-values, the false discovery rate was maintained at αfdr = 0.05 using the Benjamini and Hochberg (1995) correction to control for multiple comparisons across all ROIs. The effect size was reported as Hedges’ g (Hedges and Olkin, 2014) corrected for unequal sample sizes as per Rosnow et al. (2000). This was repeated for the three pairwise contrasts.

Differences in structural-covariance between pTBI and controls

Structural-covariance networks were generated using the Freesurfer-derived structural parcellation as the nodes (n = 68) and the edges of the network the similarity of cortical-thickness between as pairs of ROIs. As is common in the structural-covariance literature, cortical-thickness was used as the dependant variable for general linear models run across all ROIs with covariates of age at scanning, sex, and estimated total intracranial volume. This is to control for the fact that cortical-thickness has been shown to decrease with age (Magnotta, 1999), and increase with total intracranial volume (Im et al., 2008) and to differ across genders (Sowell et al., 2007). The studentised residuals were then retained for analysis and used to generate graphs of structural covariance. Pearson's correlations between residuals of each ROI generated a single 68 × 68 adjacency matrix data. This will represent an undirected, unthresholded, weighted network, with ROIs as the nodes and correlation coefficients as the edge-weights between nodes. This network represents age-invariant structural covariance (Váša et al., 2017) with age at scanning controlled for in the model.

For each graph/network, the ‘magnitude’ of structural covariance for each node was measured as node strength. For node i, this is the sum of the connectivity weights of all edges connected to node I (Fornito et al., 2016). We did not normalize these measures based on number of edges as we utilised the fully-connected, unthresholded networks and thus the number of edges connected to each node was equal across all nodes. To calculate an estimate of graph-level strength, we calculated the average nodal strength over all nodes. To generate confidence intervals for each group, these measures were bootstrapped over 5000 resamplings. In order to assess significant differences in structural covariance, permutation testing (5000 permutations) generated a null distribution of differences (t-values) in graph metrics between two groups with a two-tailed α-level of 0.05. These comparisons were conducted for each of the three pair-wise contrasts and were conducted at the graph-level (mean graph strength) and at the nodal level. p-values for nodal-level comparisons were also FDR-corrected over the 68 nodes, whilst the graph level comparisons were FDR-corrected over the three comparisons.

Structural-covariance between regions with cortical-thickness reductions in pTBI

To assess whether structural-covariance was significantly greater between regions with cortical-thickness reductions in pTBI compared to randomly selected regions, we conducted permutation testing. Briefly, for each contrast, ROIs were ranked in terms of the effect size of cortical-thickness reductions in the patient group compared to controls. For the top n-regions in terms of effect size, mean nodal strength was calculated (where n = 2,3,4…68) based on the structural-covariance graph calculated for the control group only. A null distribution of this mean nodal strength was generated by calculating mean structural-covariance for 5000 sets of randomly selected sets of n-nodes (without replacement). For each value of n, a one-tailed p-value was calculated as the proportion of permutation cases where the mean nodal strength of randomly selected nodes exceeded that of the observed mean nodal strength. p-values were corrected across values of n using the FDR-correction. A significant result suggests that structural-covariance of regions where cortical-thickness reductions exist is significantly greater than expected for randomly selected regions. We also repeated this analysis using the larger ABIDE cohort with which to provide an estimate of age-invariant structural covariance across a larger, more representative dataset compared to the experimental controls.

Results

Differences in cortical-thickness between pTBI and controls

Analyses showed that 46/68 ROIs had a cross-sectional reduction in cortical-thickness (adjusted for age at scanning, sex, and eTIV) in the patient group compared to experimental controls, with a small mean effect size (). However, no ROI showed significantly thinner cortex in patients compared to controls after FDR correction. When comparing the cognitively spared group (EFGood) to experimental controls, similar results were found, 42/68 regions showed a reduced cortical-thickness in patients compared to controls () yet no differences survived FDR correction.

There appeared a greater amount and severity of a reduction in cortical-thickness for the cognitively impaired group (EFPoor) compared to controls, 62/68 regions had an effect size in the direction of reduced cortical-thickness for patients, with the mean effect size being bigger than that of the other contrasts ().

The effect sizes of these contrasts can be seen in Fig. 2, where positive effect size indicates a reduction in cortical-thickness in the patient group relative to controls. However, it is important to remember that, across all contrasts, no regional cortical-thickness reductions in the TBI group were significant (after FDR correction).

Fig. 2

(A) Effect size (Hedge's g) maps for three planned contrasts; (i) Patients < Controls, (ii) EFGood < Controls and (iii) EFPoor < Controls, each contrast representing thinner cortical-thickness in pTBI compared to controls. (B) Permuted difference in structural-covariance for three planned contrast (i) Patients vs Controls, (ii) EFGood vs Controls and (iii) EFPoor vs Controls with greater permuted difference representing greater structural-covariance in the patient group.

Differences in structural-covariance between pTBI and controls

Mean graph strength for each of the groups and subgroups can be found in Table 4. No significant difference (after FDR correction between the three contrasts) in mean graph strength was found between patients and our experimental controls (permuted difference (PermDiff) = −0.483, pfdr = 0.062). When investigating subgroups, significant differences were found between experimental controls and EFpoor (PermDiff = 0.947, pfdr = 0.008) but not EFGood (PermDiff = −0.307, pfdr = 0.062, respectively). However, it is important to note that, whilst the observed between-group difference between EFPoor and experimental controls was significant in comparison to the permuted-distribution, the confidence intervals of the differences all crossed zero.

Table 4

Mean graph strength and bootstrappeda 95% confidence intervals.

Group	Graph Strength	CILow	CIHigh	PermDiffb	pfdrb
Controls	17.1	11.8	22.6	NA	NA
pTBI	28.0	21.9	34.8	−0.483	.062
EFGood	28.0	19.5	37.5	.307	.062
EFPoor	37.1	27.9	48.1	.947	.008

5000 resamplings.

compared to experimental controls, greater permuted difference representing greater structural-covariance in the patient group.

After FDR correction across ROIs (n = 68), no nodal differences remained significant between control and the whole pTBI group or EFGood subgroups. However, when comparing the EFPoor group to controls, multiple regions (44/68) showed significantly greater nodal strength in the patient group. These regions can be seen in Fig. 2. These regions were widely distributed across the cortex, yet a high proportion of these significant regions were found in the frontal lobe (41% frontal lobe, 25% temporal lobe, 20% parietal lobe, 9% cingulate, 5% occipital lobe).

Structural-covariance between regions with cortical-thickness reductions in pTBI

We conducted permutation testing to estimate whether, for either the whole group or either of the two subgroups, regions which showed cortical-thickness reductions compared to controls were those regions which, in the typically developing population (i.e. using the experimental controls or ABIDE data controls), show higher levels of structural-covariance. When considering the whole group of pTBI patients, for no value of n number of regions with greatestcortical-thickness reductions was the mean strength of regions in the experimental control group significantly greater than that of n randomly selected regions (see Fig. 3). This was also true of the EFGood subgroup.

Fig. 3

Observed mean strength in the experimental control group across n nodes (n = 2–68) with greatestcortical-thickness reductions in the whole pTBI group and both the EFGood and EFPoor subgroups, grey region represents the mean nodal strength for 5000 permutations of n randomly selected nodes.

However, for the EFPoor group, the mean strength in the experimental controls of the n nodes with greatest cortical-thickness reductions was significantly greater than the mean strength of n randomly selected regions for 59/67 values of n (n = 8, 10 – 67, pfdr all < 0.05).

We validated these results by repeating the analysis using the ABIDE dataset to estimate mean strength of nodes in the typically developing brain, as seen in Fig. 4. The results using our experimental controls were largely replicated; mean node strength was significantly greater than that of n randomly selected regions for multiple values of n in the EFPoor group (n = 19 – 65, 67, pfdr all < 0.05), but neither the whole pTBI sample or the EFGood subgroup (as seen in Fig. 5).

Fig. 4

(A) Nodal strength as estimated from the typically-developing controls from the ABIDE cohort across all ROIs in the Desikan-Killany atlas; (B) Distribution of nodal strength values across the ROIs.

Fig. 5

Observed mean strength in the ABIDE control group across n nodes (n = 2–68) with greatest cortical-thickness reductions in the whole pTBI group and both the EFGood and EFPoor subgroups, grey shaded region represents the mean nodal strength for 5000 permutations of n randomly selected nodes, with the grey line representing the mean nodal strength across all permutations.

Lesion cases

To examine whether the results were driven by a bias towards cases with cortical grey matter lesions who were processed using our custom Freesurfer pipeline, we repeated all analyses excluding the lesion cases (n = 8). The results can be seen in supplementary material but, briefly, these were qualitatively the same as the results presented above, with significant effects seen in the EFPoor group but not the whole group or the EFGood group.

Discussion

The findings presented here demonstrate a potential network-based mechanism for the association between (sub)acute cortical thinning and long-term executive dysfunction. For those patients who exhibited poor, long-term EF outcomes, cortical-thickness reductions were more likely to be localised to regions that typically have higher structural-covariance than randomly selected regions in normal development. In other words, we found that, at the group-level, for cases where long term executive function outcome is poor, cortical damage (measured as cortical-thickness reductions) is seemingly preferentially loaded onto regions high in structural covariance in the typically developing population, but not when EF is spared. A key strength of this study was that we were able to replicate these findings across two different control groups.

In traditional ROI-based analyses of cortical thickness, we found that no regions showed significant thinning, regardless of chronic executive functioning. However, our findings of network-based loading of the subtle cortical differences in the executive dysfunction group, suggests that rather than the topography (the physical distribution across cortex) of specific alterations being important to functional outcomes, it is the topology (the connectivity of a node) of these regions in the wider cortical network which is important. There is spatial inconsistency in alterations to brain morphometry associated both between and within clinical manifestations of neurological disorders (Cauda et al., 2019), including pTBI (King et al., 2019), and the location and extent of focal brain damage is seemingly insufficient to fully explain the neuropsychological deficits that persist post-injury (Bigler, 2001; King et al., 2019).

Our findings suggest that, this may be because spatially-disparate lesions and/or morphometric changes may occupy similar topological positions in the network and in these cases result in a similar neuropsychological profile. Damage to topologically central regions likely has a disproportionate impact on the broader network, especially as they are particularly relevant to the development of the brain (Csermely et al., 2013; Morgan et al., 2018; Oldham and Fornito, 2019), and this in turn renders this ‘damage’ more likely to be behaviourally symptomatic (Crossley et al., 2014; Hillary and Grafman, 2017).

This was apparent in these results with damage loading preferentially onto these regions seemingly more likely to result in EF impairment. The neuropsychological nature of executive functions are also more likely to amplify this effect, given that they are supported by widely distributed neural networks (Beauchamp et al., 2011; Collette et al., 2006; Nowrangi et al., 2014; Slomine et al., 2002) and are therefore particularly vulnerable to the distributed, multifocal mechanisms of TBI (Treble-Barna et al., 2017).

The findings here, if replicable, could have significant clinical implications for pTBI. This association identified by these results, that the topological positioning of (sub)acute cortical thinning being related to long-term (2 years post-injury) cognitive sequelae suggests that the structural-basis of poorer outcomes is set in motion from the early damage caused by primary injury, specifically the topology of this damage. Whilst the effects of secondary injury are treatable and potentially reversible with prompt and proper disease-modifying treatment (Ghajar, 2000; McKee and Daneshvar, 2015), primary injury is considered more permanent.

Therefore, these results suggest that some variance in chronic poor cognitive outcomes can be explained by the unchangeable (sub)acute, primary injury. This would therefore mean that there is potential for predicting, at this early timepoint, those young people at greatest risk of these long term outcomes potentially allowing for early identification of these young people, and thus early neurorehabilitative treatment, ensuring these young people can maximise their long-term outcomes.

It is important to note the acute/sub-acute timing of the MRI in these cases. MRI changes to both cortical thickness and volume have been reported in both the acute and sub-acute stages of pTBI previous findings at this acute time-point post injury (King et al., 2019; McCauley et al., 2010; Urban et al., 2017; Wilde et al., 2012), however, we did not replicate this in our current sample. Our findings also challenge previous analyses of post-injury cortical-thickness being specifically related to executive functioning (Wilde et al., 2012). Observed differences were generally in the expected direction, with patients having thinner cortices, but effect sizes were very small across regions.

Cortical thickness reductions, as measured in-vivo with MRI, aim to assess the potential atrophic effects of the cascade of mechanisms that occur post-injury (Bigler, 2013). At this early stage post-injury, cortical thinning could be due to primary injury mechanisms, such as impaired perfusion (van der Kleij et al., 2020; Xu et al., 2010) or trauma induced cell loss (Bigler, 2013; Cullen et al., 2011) for instance. Cortical thinning is likely to be associated with loss of cortical neurons (Maxwell et al., 2010), although this continues over months and years post-injury, likely due to the ongoing effects of secondary injury (Bigler, 2013) and potentially explaining neuroimaging findings consistent with chronic neurodegeneration over time (King et al., 2019). The lack of significant findings in terms of cortical thinning may suggest that, without the chronic effects of neurodegeneration, cortical thinning is relatively subtle and undetectable (especially in smaller samples). It may be the case that the changes seen this acutely post-injury are too subtle to detect at these sample sizes, especially over multiple ROIs. However, the only two studies to report cortical-thickness reductions in patients compared to controls at a similarly acute time point post-injury conducted vertex-wise analyses (McCauley et al., 2010; Urban et al., 2017). It may in fact be that cortical thickness reductions are local to the sites of primary injury (coup and contra-coup) and the process of parcellating and averaging these changes over a region for ROI-based analyses makes these harder to detect.

Despite the non-significance of these univariate tests, we hypothesised that the cumulative effect of these subtle differences in cortical-thickness still has a functionally meaningful effect on the developing brain, post-injury. Hence, in a novel set of analyses, we investigated the structural covariance network as a method to investigate the multivariate relationships between cortical-thicknesses across the cortex. We found that differences in the mean graph strength, the average magnitude of structural-covariance across all nodes, was not different when comparing pTBI patients to our experimental controls, however, when stratifying based on EF outcome, significant differences from controls were found for the EFPoor but not the EFGood group. This suggests that the structural covariance network is only ‘abnormal’ in the impaired group, with the non-impaired group showing a network structure similar to controls. This pattern was repeated for the nodal-level findings and, somewhat unsurprisingly, the significant differences in the EFPoor group were overly represented by nodes in the frontal and temporal lobes, regions commonly implicated in brain morphometry differences post-pTBI (King et al., 2019), with frontal regions being key regions in the widely distributed neural networks supporting executive functions (Beauchamp et al., 2011; Collette et al., 2006; Nowrangi et al., 2014; Slomine et al., 2002).

At both the nodal- and graph-level, the magnitude of structural-covariance was greater in the impaired patient group than our experimental groups. This would suggest that, in these patients, the morphometry of regions across the cortex was less differentiated. Whilst it remains unclear how this may translate into changes to the underlying neuroanatomy of the brain, this represents a marked change from the gene-controlled patterning of structural-covariance across the cortex (Romero-Garcia et al., 2018; Yee et al., 2017). This group-level finding aligns with our previous individual-level analysis of deviation from typical structural covariance network being related to EF abilities (King et al., 2020).

Given the acute timing of the MRI in this study, it is unclear how the ongoing neuropathophysiology of brain injury may alter the ongoing development of the morphometry of the brain after the injury. This will be of particular interest to the study of structural covariance after pTBI, given the developmental drivers of this higher-order organisation of brain morphometry (Alexander-Bloch et al., 2013; Fan et al., 2011; Khundrakpam et al., 2017, 2016, 2013; Raznahan et al., 2011; Váša et al., 2017), future longitudinal research would be well positioned to answer this.

Diffusion MRI (DWI) is the current standard for estimating structural connectivity of white matter fibre-bundles between brain regions (Batalle et al., 2018). White-matter connectivity of the brain is specifically susceptible to the effects of pTBI, due to diffuse axonal injury, reflected in topological differences after injury, compared to controls acutely after mild pTBI (Yuan et al., 2015) but also across wider injury severities chronically after injury (Caeyenberghs et al., 2012; Konigs et al., 2017; Yuan et al., 2016). However, high quality DWI sequences have long acquisition lengths and thus may not be suitable for paediatric populations (Batalle et al., 2013). Thus, structural-covariance may be a more feasible alternative to studying the structural effects of pTBI, by investigating the meaningful meso‑scale organisation of brain structure across the cortex. Previous research has also highlighted the potential role of WM connectivity as a driver of structural covariance between regions (Gong et al., 2012; Reid et al., 2016) as regions which are similar in cytoarchitectural organisation are more likely to be anatomically connected (Goulas et al., 2017; Wei et al., 2019). Therefore, the structural-covariance approach may in fact capture not only pathological grey-matter alterations, but also the effects of DAI. However, it is vital to emphasise that DWI and associated tractography methods estimate ‘actual’ connectivity, rather than the network-level organisation of brain structure measured by structural covariance. These methods are thus likely to provide overlapping and complementary information to DWI, and future studies should combine these methodologies in multimodal studies of the cortex post-injury to better understand how they capture injury mechanisms.

Limitations and future considerations/directions

It is important to note that in this cohort, there were a limited number who met the criteria we set for executive dysfunction, compared to those with a favourable outcome at two years post-injury. This means that, the sub-analysis had much smaller sample sizes, especially the EF-Poor group, although we conducted a number of steps to best support the findings reported here (i.e. validating in a second control cohort, correcting for unequal sample sizes etc.). Future research should firstly validate the claims made here before any work can be done to investigate the clinical significance of this work.

Structural-covariance is limited to the application to population-level covariance in neuroanatomy (Alexander-Bloch et al., 2013). We specifically investigate the age-invariant structural covariance network (Váša et al., 2017), since the analysis combines data across childhood and adolescence, modelling the common network structure across development. The limited numbers in our experimental control group could result in limited accuracy of the estimation of this age-invariant structural-covariance network. Utilising the ABIDE reference group allowed us to replicate the results using a more reliable estimate of the age-invariant structural-covariance network due to the much larger sample size. However, it is important to note that hubs of the structural-covariance network have distinct developmental trajectories over the time course of childhood and adolescence (Khundrakpam et al., 2013; Whitaker et al., 2016). Therefore, future research should try to resolve the relationship between cortical-thickness reductions post-TBI and structural-covariance across age-matched structural-covariance networks.

It is important to consider the potential methodological reasons for detecting differences in structural covariance but not cortical thickness. It may in fact be the case that the SCN approach, due to the fact it considers higher order interactions and organisation of regional morphometry, these methods are more sensitive to the subtleties of injury. This may be because (sub-significant) thinning to region X may not be detected in a standard univariate analysis, but in the SCN approach, this thinning impacts on nodal metrics for all nodes – thus amplifying these subtle differences. This amplification and thus increased sensitivity to structural differences may be beneficial in studies of mild traumatic brain injury, where damage can be difficult to detect even when there is significant functional impairment caused by injury. Future work should further investigate the role of SCN analyses in TBI, especially in more chronic populations where even univariate cortical thickness differences are detectable, to identify the additional benefit of these methods. Methodologically, greater consideration also needs to be given to the apparent power of SCN analyses to detect abnormalities in comparison to traditional univariate analyses.

We focussed on nodes with high topological strength in the network, namely regions that have high summed structural covariance with all other regions of the brain. This metric is a relatively simplistic, although intuitive, measure of ‘hubness’ but may not capture the more nuanced aspects of the centrality of a node in a network (Oldham and Fornito, 2019; van den Heuvel and Sporns, 2013). However, due to the fact that the structural-covariance networks do not adhere to typical assumptions of networks (edges representing definitive real connections) we utilised strength as a simpler metric which makes fewer assumptions about the underlying neurophysiology of the network. Once a more complete understanding of communication dynamics throughout the structural covariance network has been understood, future studies may investigate other, more nuanced measures of nodal centrality, which may capture greater information about their role in the wider network.

We specifically try to elucidate effects specific to pTBI, and where possible have accounted for variables which may also impact upon cortical development such as gender and age. However, other important factors may have influenced both cortical development and/or post-injury executive functioning such as pre-morbid IQ. Accurate estimates of pre-morbid functioning where not collected in the current dataset and reflects a common limitation of research in this field.

Conclusion

There is strong theoretical support for future studies of brain insults to focus on generating hypotheses about underlying pathophysiology to the neural network biology (Aerts et al., 2016). Here we provide a novel methodological advance in support of this goal, offering an analysis of structural-covariance post-TBI in which we specifically proposed that the topology of nodes which were damaged would be important for understanding which children go on to experience functionally relevant impairment post-injury. Given disparity in outcomes after a pTBI, these results are a key first step in utilising knowledge of topology in typical development for predicting which children go on to experience significant impairment post-injury as opposed to those who will recover well. Future research needs to expand these findings to investigate the causal nature of these relationships, and to see whether these patterns expand beyond the structural-covariance network of the brain (i.e. DTI).

CRediT authorship contribution statement

Daniel J. King: Conceptualization, Methodology, Formal analysis, Writing – original draft. Stefano Seri: Writing – review & editing, Supervision. Cathy Catroppa: Investigation, Resources, Writing – review & editing. Vicki A. Anderson: Investigation, Resources, Writing – review & editing. Amanda G. Wood: Conceptualization, Writing – review & editing, Supervision, Funding acquisition.