Developmental divergence of structural brain networks as an indicator of future cognitive impairments in childhood brain injury: Executive functions.

Brain insults during childhood can perturb the already non-linear trajectory of typical brain maturation. The diffuse effects of injury can be modelled using structural covariance networks (SCN), which change as a function of neurodevelopment. However, SCNs are estimated at the group-level, limiting applicability to predicting individual-subject outcomes. This study aimed to measure the divergence of the brain networks in paediatric traumatic brain injury (pTBI) patients and controls, and investigate relationships with executive functioning (EF) at 24 months post-injury. T1-weighted MRI acquired acutely in 78 child survivors of pTBI and 33 controls underwent 3D-tissue segmentation to estimate cortical thickness (CT) across 68 atlas-based regions-of-interest (ROIs). Using an 'add-one-patient' approach, we estimate a developmental divergence index (DDI). Our approach adopts a novel analytic framework in which age-appropriate reference networks to calculate the DDI were generated from control participants from the ABIDE dataset using a sliding-window approach. Divergence from the age-appropriate SCN was related to reduced EF performance and an increase in behaviours related to executive dysfunctions. The DDI measure showed predictive value with regard to executive functions, highlighting that early imaging can assist in prognosis for cognition.

The pathological effects of neurological conditions occurring during childhood, necessarily interact with the highly-programmed maturation of the brain, perturbing the trajectory of normal brain development, which is in itself non-linear (Gogtay et al., 2004; Shaw et al., 2008). Previous research has suggested that deviations from the developmental trajectory of the brain may act as a marker of brain health, neurological disorders and cognitive functioning (Bigler, 2013; Cole and Franke, 2017; Erus et al., 2015). Thus, the degree to which the injury alters normal development may be an important factor to consider when trying to understand subsequent cognitive sequalae post-insult including impairments to intellectual and executive functioning, as well as attention and processing speed (Crowe et al., 2015). The current study investigates this idea using a measure of divergence of the structural network to investigate levels of post-insult cognitive impairment, with a focus on executive functioning.

Specifically, the current study focuses on traumatic brain injury (TBI) in childhood and adolescence, a leading cause of disability (World Health Organization, 2006). Many injuries occur in the context of a still-developing brain (Wilde et al., 2012a,b), with an incidence between 1.10–1.85 cases per hundred for the 0–15 age range (McKinlay et al., 2008). Paediatric TBI (pTBI) has specific adverse effects on neurodevelopment. The traumatic, external force to the brain can result in pathology at both a cellular and tissue level, leading to transient or even permanent impairment (Bigler, 2007, 2016; Maxwell, 2012). Some damage is realised as trauma-related, developmentally inappropriate atrophy (Bigler, 2013; Urban et al., 2017; Wilde et al., 2005) which, when imaged using techniques such as structural magnetic resonance imaging (sMRI), can appear as relative decreases to both brain volume (Bigler, 2016) and cortical thickness (CT) measures (Urban et al., 2017). However, in pTBI, these negative consequences of injury occur against a backdrop of ongoing age- and development-dependent changes to the brain (Bigler, 2016; Maxwell, 2012) leading to differential vulnerability to TBI depending on the developmental stage at which injury occurs (Anderson et al., 2011; Goldstrohm and Arffa, 2005; McCrory et al., 2004). For example, the state of development of myelinated axons at the time of injury influences the magnitude of degeneration of nerve fibres (Adelson and Kochanek, 1998; Kochanek et al., 2000; Maxwell, 2012; Staal and Vickers, 2011). Thus, disruption at different ‘critical’ periods of the developmental trajectory could result in very different functional outcomes long term (Anderson et al., 2011; Resch et al., 2019).

Previous sMRI studies have shown that, from early to post-chronic timepoints post-injury, the morphometry of the injured brain differs from that of typically developing children (see King et al. (2019) for a systematic review of findings). These cross-sectional differences are found even up to 10 years post-injury (Beauchamp et al., 2011; Serra-Grabulosa et al., 2005) suggesting alterations which 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.

Whilst these cross-sectional studies can provide evidence that differences exist, longitudinal studies are needed to provide explanation of the basis of these changes (ie whether pathologic-injury related change or developmental change) and if they resolve over time. Longitudinal morphometric studies of paediatric cohorts have investigated changes between patients and controls across multiple timepoints post-injury (Dennis et al., 2017a,b; Dennis et al., 2016; Mayer et al., 2015; Wilde et al., 2012a,b; Wu et al., 2010, 2018). The majority of these studies show a reduction in volume or cortical thinning over time in the TBI group, as well as cross-sectional differences from controls. Interestingly however, they also show an interaction between group (patient vs. controls) and time post-injury on cortical thickness (CT) measures (Mayer et al., 2015; Wilde et al., 2012a,b) and corpus callosum volumes (Wu et al., 2010), with greater atrophy over time seen post pTBI. Dennis et al. (2017a,b) also found differences in longitudinal morphometric change between a TBI patient group who experienced slowed inter-hemispheric communication, compared to those with a normal inter-hemispheric transfer time, suggesting that these structural differences are not only a result of an injury, but also relevant to post-injury functioning of the brain. Due to the highly programmed trajectories of white matter (WM) and grey matter (GM) development during childhood and adolescence (Batalle et al., 2018; Mills et al., 2016; Raznahan et al., 2011a,b; Shaw et al., 2008), these group differences in longitudinal change (between TBI patients and controls) suggest that the developmental trajectory of the brain is in fact altered to some degree by a TBI. However, previous research has not investigated the magnitude to which a pTBI interferes with the developmental trajectory at an individual level or how this may change as a function of age at which the injury occurs. Overall, these studies suggest that pTBI has potentially lifelong consequences, owing to the persistent and ongoing differences to the structural development of the cortex post-injury.

The effects of a pTBI on the brain are highly diffuse with morphometric differences found across widespread brain regions (in both cortex and subcortically) even within a single individual (Bigler, 2007; Bigler et al., 2013; Bigler and Maxwell, 2011). This diffuse (rather than focal) injury can also vary in location across individuals, samples and studies, although commonly fronto-temporal regions are affected (King et al., 2019). Thus, previous studies investigating regions of interest (ROIs) with a univariate approach, may not capture the multivariate and heterogeneous nature of injury. One way to interrogate the multivariate structure of the brain is the structural covariance network (SCN) approach (Bigler, 2016; Lerch et al., 2017), modelling the degree to which the morphology of brain regions statistically co-varies across all possible pairs of ROIs (Alexander-Bloch et al., 2013a, b; Evans, 2013; Mechelli et al., 2005).

These whole-brain, network approaches to morphometric data, within a graph theoretic framework (Bullmore and Sporns, 2009), will allow us to investigate additional information beyond that which is offered by univariate, local approaches (Bullmore and Sporns, 2009; Pagani et al., 2016).

SCNs are both biologically meaningful and sensitive to changes to the developing-brain. The topological organization of these networks are quantifiably non-random and complex (Alexander-Bloch, Giedd, et al., 2013; Evans, 2013), similarly to brain connectivity networks estimated from both resting-state fMRI and diffusion-weighted imaging (DWI). Structural covariance across the cortex changes as a function of neurodevelopment, age and maturational change (Alexander-Bloch, Raznahan, et al., 2013; Fan et al., 2011; Khundrakpam et al., 2017, 2016; Khundrakpam et al., 2013; Raznahan et al., 2011a,b; Váša et al., 2017) and may be related to shared expression of genes related to controlling cortical development between ROIs (Romero-Garcia et al., 2018). Age-related change in the SCN captures variation in changes to the brain beyond that of neurodevelopmental processes such as thinning and myelination (Váša et al., 2017). These networks are also sensitive to differences due to other types of paediatric brain insult, including malformations of cortical development and cortico-genesis disruption in neonates, with degree of differences changing as a function of the specific gestational-timing of disruption (Hong et al., 2017), suggesting that the SCN can index divergence of the typical maturational trajectory of the cortex. Thus, these approaches may allow us to capture the developmental-divergence of morphology after a pTBI and investigate its relationship to functional outcomes post-injury.

However, as multiple participants are required to sample enough cortical measurements to generate a correlation between all possible regional-pairs (each participant can only contribute a single measurement per region), this SCN approach can only generate group-level brain networks, expressing population-level covariance in neuroanatomy (Alexander-Bloch, Raznahan, et al., 2013). Thus, studies have tried to develop methods that can translate this information to the individual-subject level (for example Seidlitz et al. (2018) or Tijms et al. (2012). We specifically utilise the individual contribution metric (proposed by Saggar et al. (2015)), as a potential solution to this problem, which allows us to estimate distance of a patient from a group-level, reference SCN.

Aims and hypotheses

We aimed to measure subject-level divergence of the structural covariance network following brain insult and potential perturbation of brain development. Specifically, we investigate deviation in a cohort of paediatric TBI patients from a reference network of typically developing control participants by leveraging a large-scale, open-access MRI database. Our approach adopts a novel analytic framework of a sliding-window approach to calculate these developmentally-appropriate reference networks. We predict that there will be greater divergence of structural networks for cases with a pTBI compared to control cases. We also aimed to use these divergence metrics as a proxy of perturbations in brain development and as a predictor of long-term functional outcome, specifically hypothesising that greater structural divergence will be associated with greater executive dysfunction. The current study focused upon executive functioning because cognitive-skills are more vulnerable to damage occurring during the period of skill-maturation (Ewing-Cobbs et al., 2004; Krasny-Pacini et al., 2017), thus the protracted period of EF development (Diamond, 2013; Friedman et al., 2016; Perone et al., 2018) means EF is likely to have an extended window of vulnerability (Krasny-Pacini et al., 2017).

We also hypothesised that stronger associations would be found between structural divergence and executive dysfunction when investigating a sub-graph of the whole-brain SCN, which consists of regions known to subserve core executive function skills.

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 for a study which had previously received ethical approval via the Human Research and Ethics Committee of Royal Children’s Hospital, Melbourne, Australia. The reference data used in this research was acquired through the public Autism Brain Imaging Data Exchange (ABIDE) database, as shared by the Preprocessed Connectome Project (PCP). The database has de-identified all the patient health information associated with the data. A favourable ethical opinion was granted by Aston University 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, controls n = 43) 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).

TBI severity was categorized as follows: mild TBI: GCS 13–15 on hospital presentation, no evidence of mass lesion on CT or clinical MRI and no neurologic deficits (if there was evidence of intra-cranial pathology, these were classified as mild complicated); moderate TBI: GCS 9–12 on hospital presentation, and/or mass lesion or other evidence of specific injury on CT/MRI, and/or neurological impairment; and, severe TBI: GCS 3–8 on hospital presentation, and/or mass lesion or other evidence of specific injury on CT/MRI, and/or neurological impairment.. Due to small group sizes in relevant analyses, the mild-complicated, moderate and severe groups were collapsed for analyses.

MR-Images were acquired for the patient group acutely after injury (<90 days post-injury, range = 1–88 days). MRI images were acquired at 3 T as a part of an existing research protocol on a Siemens Trio scanner (Siemens Medical Systems, Erlangen, Germany) using a 32-channel matrix head coil. The standard 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 degrees; BW 200 Hz/Px; 176 slices; resolution 1 × .5 × .5 mm] and sagittal 3D T2-w 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 × .5 × .5 mm resolution matched in alignment to the 3D T1-weighted sequence].

We applied a number of inclusion criteria to the dataset, only including subjects who; a) met strict quality control criteria of Freesurfer outputs, b) had no gross/frank pathology/lesions identified within the grey matter ribbon (as this may bias image processing with Freesurfer (King et al., 2020), c) had available MRI data and were scanned <90 days post-injury. This resulted in a subset of n = 108 subjects (TBI patients (n = 75) and healthy controls (n = 33)). Group demographics can be seen in Table 1.

Table 1

Demographics for each cohort by group.

Cohort/Group	TBI Cohort -Patients	TBI Cohort - Controls	ABIDE
N	n = 75	n = 33	n = 327
M/F	51/24	20/13	259/68
Age at MRI (median, yrs)	10.81	9.99	12.49
(range, yrs)	6.18-14.91	6.53-15.47	6.47-16.93
Age at Injury (median, yrs)	10.58	–	–
(range, yrs)	6.08–14.67
Injury-MRI interval (median, days)	34
(range, days)	1-88
Injury Severity
Mild	47	–	–
Moderate/Severea	28	–	–

Mild Complicated TBI + Moderate TBI + Severe TBI.

ABIDE dataset

In order to provide a healthy reference group for the calculation of our divergence metric, 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), for full details see Pre-processed Connectome Project website 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 on the ABIDE website (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 sMRI (see supplementary materials for further details), 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. As per ABIDE’s recommendations to share the data ID list used for primary analyses, this can be found in supplementary materials. Group demographics can be seen in Table 1.

Both controls in the experimental cohort and the ABIDE cohort had 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 (Wechler, 1999) whereas the measures used by the ABIDE dataset were varied, see ABIDE documentation for details).

MRI processing

3D tissue segmentation and estimation of CT from T1-weighted (T1w) MR images was 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). Where available, 3D T2-weighted (T2w) FLAIR MRI were used to refine the boundary between the pial surface and dura. 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). The quality of Freesurfer outputs was assessed using Qoala-T (Klapwijk et al., 2019) as a decision support tool to guide systematic and replicable selection of which cases required manual editing. 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 Freeurfer for the ABIDE data, and steps to control for ABIDE site and cohort effects (TBI cohort vs ABIDE), can be found in Supplementary materials.

Graphs of structural covariance

All network analysis were conducted with a series of packages in R version 3.5.0 (R Core Team, 2016), specifically brainGraph version 2.2.0 (Watson, 2016), which is an expansion of the iGraph package (Csardi and Nepusz, 2006).

As is common in the SCN literature, CT 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 CT 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 for the ABIDE reference data. This will represent an undirected, unthresholded, weighted network, with ROIs as the nodes and correlation coefficients as the edge-weights between nodes.

Divergence metrics

Since structural covariance networks are derived from correlations between regions within participants, graphs are compiled at a group level. Our hypotheses suggest that the individual deviation from ‘typical’ maturation will be an important variable in the prediction of executive function. Therefore, it is important to identify methods by which to extract estimates of this deviation at the individual subject level as a proxy of perturbations in brain development. A developmental-divergence index (DDI) is therefore generated for each patient using the ‘Add-One-Patient’ (AOP) approach (Saggar et al., 2015). This measure is further outlined below. Saggar et al. (2015) term this ‘individual contribution’ to the group-level network. From our perspective, those that are most different from the group/reference network will be those whose development is furthest from typical, expected trajectories. Hence, we refer to this ‘individual contribution’ metric as a (developmental) divergence index.

The AOP approach allows the direct comparison of the weighted SCN by assessing the matrix of CT residuals. The approach compares the structural network of a reference group and a second matrix comprising of the reference group, plus a single patient (hence AOP). This means that the existing correlation matrix for a reference control group, denoted R, will be combined with each patient individually, to generate a new matrix, denoted R + P (where i is the individual patients, i = 1, 2,. .n). Subsequently, a normalized Mantel test (Mantel, 1967) is conducted to assess similarity of these matrices calculated as;Where X and Y represent and respectively, and are elements of these matrices, and are the standard deviations for the matrices and n is the number of nodes (in the case of this study, 68 ROIs) for each correlation matrix (Saggar et al., 2015). This metric of similarity (whereas r increases this represents two matrices with higher similarity) is subtracted from one to compute the divergence from the reference group matrix where;

These divergence metrics will provide individual-level distance (at the level of the whole graph) from the reference group. If and are similar (Mantel test trending toward 1), subj(i) has not altered the group-level network, and therefore Subj(i) does not show divergent morphology, thus DDI will be low. If highly dissimilar (Mantel test trending toward 0), addition of subj(i) has significantly altered the group-level network, thus subj(i) is different from typically developing peers (and DDI is greater). Essentially, if the patient exhibits developmentally-appropriate morphometry, the reference-plus-patient network will be similar to that of the reference group alone. Therefore, the less similar the networks, the more developmentally divergent the patient’s morphometry. Thus, for each patient the analysis will output a single DDI to estimate divergence across the whole cortex.

Reference networks

In order to calculate developmental divergence for both the control and TBI cases from the TBI cohort, we used the ABIDE dataset as a reference group to generate the R SCN. We calculated developmental divergence from the typically-developing SCN using two approaches, an age-invariant SCN and an age-matched SCN. For the age-invariant network, the all participants selected for our ABIDE cohort were utilised as a reference group. Similarly to Váša et al. (2017), we termed this the age-invariant SCN since the analysis combines data across childhood and adolescence, with participants of all ages in the ABIDE sample included. Age-invariant DDI (DDIinv.) was therefore calculated for each subject in the TBI cohort (both patients and controls) from this whole-group reference. Previous studies have adopted this approach of using a single reference group to calculate Saggar et al.’s (2015) ‘individual contribution’ metric. A single reference group combines a wide range of ages and thus ignores known variations in developmental changes of grey and white matter across childhood and into adolescence (Gogtay et al., 2004; Sowell et al., 2004). Thus, we also adopted a novel analytic framework in which developmentally-appropriate reference networks to calculate an age-matched DDI (DDIage) were generated from control participants from the ABIDE dataset using a sliding-window approach (outlined below).

Sliding-Window

Similarly to Váša et al. (2017), we used a sliding window approach in order to calculate developmentally appropriate, age-matched reference SCNs. In brief, subjects form the ABIDE dataset were ordered by age at scanning. Subject-level CT residuals were then correlated within equal-sized windows of participants, with the window being ‘slid’ across the age-range of participants (Váša et al., 2017). A window-size of 26 participants and a step-size of 15 participants was selected, subsequently 21 half-overlapping windows across the ABIDE cohort were selected, resulting in a single reference SCN per window.

Window-size was selected against a number of criteria: a) based on recommendations by Saggar et al. (2015) in relation to stability of their AOP metric, b) maximised the difference-statistic for control vs TBI differences in DDIage measure, and c) which resulted in an nth window (where number of windows is 1: n) which was as close to the defined window size as possible (due to the remainder from the calculation of:the final window was not guaranteed to have the full number of subjects). Details of this window-size selection process can be found in supplementary materials.

Once reference SCNs for each window were generated, the median age of participants within the window were calculated. Each participant within the TBI cohort (patients and controls) was individually-matched to the reference window which minimised the difference between their age at scanning and the median-age of the reference-window. This matched reference window was then used to calculate the DDIage for that individual.

Executive functions (EF)

We investigate EFs as they are commonly impaired, both acutely and chronically post-injury, but also because they show a protracted period of maturation and development (Diamond, 2013; Friedman et al., 2016; Perone et al., 2018) and are therefore likely to have an extended window of vulnerability to the effects of injury (Krasny-Pacini et al., 2017). EF was assessed for pTBI patients at approximately 24-months post injury (M(SD) = 754(80) days post-injury). EF was assessed in controls relative to their MRI scan (M(SD) = 367(135) days post-MRI). EF was assessed both using performance-based neuropsychological testing and a parent-reported measure.

Several standard 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.

Table 2

Neuropsychological tests and subtests used to calculate EF scores.

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

Performance scores for the neuropsychological test batteries were converted to age-scaled scores (M = 10, SD = 3). To provide a summary score for common EF performance, we summed these age-scaled scores across subtests, with higher scores representing better performance. The EF summary score was calculated for 80 subjects (TBI n = 52, controls n = 28) who had data for all subtests available, as well as sufficient data to calculate the DDI. This summary score was used for two main reasons; firstly, due to a limited sample size and the use of correlational analyses, we have limited power to look at each domain separately. Secondly, due to the fact we are using a whole-brain measure of developmental–divergence, it is likely that the measure is too coarse to capture the nuances across multiple sub-domains of executive functioning.

The Behaviour Rating Inventory of Executive Function (BRIEF, Gioia et al. (2000) measures EF in daily life, using purposeful, goal-directed behaviours to solve and adapt to problems (Donders and DeWit, 2017). The current study specifically uses the ‘Global Executive Composite’ T-score (GEC; M = 50, SD = 10), with higher scores representing greater difficulties in behavioural EF (measured in TBI n = 52, controls n = 32). By using two differing measures of EF (performance-based vs behavioural/parent report) we are able to assess concordance of our results across multiple measures.

Central Executive Network (CEN)

The DDI represents the divergence of a subject’s morphology from the typical SCN across all cortical nodes/ROIs. However, this may reduce the signal to noise ratio when looking at associations between DDI and EF, as not all regions may be relevant to subsuming EF. Thus, we also investigated DDIinv. and DDIage across a subgraph of the SCN, specifically regions within the CEN. The CEN is a neural network that shows heightened activity during typical tasks of EF (Seeley et al., 2007; Sherman et al., 2014; Thomason et al., 2011). We defined the CEN anatomically as per the Desikan-Killany (Desikan et al., 2006) atlas regions identified in Ryan et al. (2017), which comprised regions of dorsolateral pre-frontal cortex and posterior parietal cortex. Specifically, regions were bilateral caudal and rostral middle frontal gyrus, inferior and superior parietal lobule, precuneus and superior frontal gyrus. These regions have been identified (amongst others) as supporting common EF activation in adolescence and childhood (Horowitz-Kraus and Holland, 2016; McKenna et al., 2017) and adulthood (Niendam et al., 2012).

Statistical analysis

All analysis were performed in R (R Core Team, 2016) using the ‘stats’ packages. Analyses were planned a-priori as follows. Due to the non-normal distribution of the DDI metrics, both DDIinv. and DDIage were compared between patients and controls from the TBI cohort using a one-tail Mann-Whitney test (with the alternative hypothesis of the location shift of mean DDI from controls to the patient group is greater than 0). Pearson’s correlations were calculated between DDI measures (both DDIinv and DDIage) and each EF measure EF (EF score and BRIEF). This was calculated for both DDI calculated on the whole network and calculated on the CEN sub-graph. Correlations were calculated for the whole sample, and independently for patient and control groups. The sample sizes for the current study were larger than many current pTBI MRI studies (King et al., 2019) however, we acknowledge that this could still influence statistical analyses. Thus, both resampling approaches and false discovery rate (FDR) correction were used to mitigate these risks. The bootstrapped (100 iterations) 95 % confidence intervals (CIs) were calculated for all point estimates of correlation coefficients. Raw p-values calculated using a permutation resampling approach (5000 permutations, calculated using the jmuOutlier package in R version 2.2) are reported. Significance was assessed using an FDR correction (Benjamini and Hochberg, 1995). Results are presented using the ‘ggplot2′ (Wickham, 2009) and ‘ggpubr’ (Kassambara, 2018) packages.

Post-hoc analyses

A number of analyses were conducted post-hoc to assess the robustness of the approach. Firstly, split-half analyses were conducted to assess the internal reliability off the DDIinv across different subsets of the normative reference group. Briefly, the ABIDE dataset was randomly split into two groups (n = 164 & 163) and DDIinv was calculated for all pTBI patients using both the 1 st and 2nd halves of the ABIDE sample and the Pearson’s correlation between these is reported. This was repeated for 500 random split halves. Additional comparisons investigated whether DDIinv and DDIage differed as a function of injury severity. To maintain statistical power, mild-complex, moderate and severe injury classifications were grouped into a ‘Moderate/Severe’ group for comparisons. Clinical presentation between injury severities is very different and thus treating the patient group as a single cohort in patient vs control analyses of the divergence index may miss clinically meaningful differences. Finally, partial correlations (Pearson’s) were conducted between whole-brain DDIinv/DDIage and EF/BRIEF whilst controlling for age at scanning (yrs), to control for potential age-related biases in these measures and also simultaneously controlling for both age at scanning (yrs) and the interval between MRI and EF assessment (days).

Results

Age-invariant network and DDIinv

Median DDIinv. for the TBI and control group were 5.16e-05 (min = 1.08e-05, max = 8.47e-04) and 3.97e-5 (min = 1.37e-05, max = 1.90e-04), respectively. Violin plots of DDIinv. for each group can be seen in Fig. 1a. The difference of DDIinv. from the TBI group compared to the control group was not significantly greater than zero (W = 1046, p = .890). Given that divergence from the whole-group reference SCN may be due to the difference between age of subjects and the median age of the reference network, we plotted this absolute difference against DDIinv. No apparent relationship was found (in the TBI or control group), as can be seen in Fig. 1b. No significant association was found between DDIinv. and age at injury in the patient group. In terms of association with EF at two years post-injury, DDIinv. was significantly negatively correlated with EF performance across the whole sample (r = −.300, p = .009), but specifically in the TBI population (r = −.319, p = .024), and not controls. DDIinv. was significantly, positively correlated with BRIEF GEC in the whole sample (r = .277, p = .021). No significant relationships were found with the BRIEF GEC in the TBI group (see Fig. 1c and d).

Fig. 1

a) Violin plots of DDIinv. for both TBI and control groups, b) correlation between age at injury and DDIinv., and correlations between EF and DDIinv., specifically c) executive function score and d) BRIEF GEC, for both the TBI and control groups.

Age-matched network and DDI

Median DDIage for the TBI and control group were 4.583e-03 (min = 9.75e-04, max = 7.33e-02)_and 4.14e-03 (min = 1.38e-03, max = 1.56e-02), respectively. Mean absolute difference between age at scanning of the subject and the median age of the window that they were matched to was .30 yrs for both TBI (SD = .50) and control (SD = .43) groups. Violin plots of DDIage for each group can be seen in Fig. 2a. The difference between DDIage from the TBI group and the control group was not significantly greater than zero (W = 1181, p = .648). No significant association was found between DDIage and age at injury in the patient group. In terms of association with EF at two years post-injury, DDIage was significantly negatively correlated with EF performance across the sample (r = −.308, p = .007), and in the TBI population (r = -.330, p = .021), but not controls. Significant positive relationships were found with the BRIEF GEC (r = .277, p = .021) but not the TBI or control populations (see Fig. 2c and d).

Fig. 2

a) Violin plots of DDIage for both TBI and control groups, b) correlation between age at injury and DDIage, and correlations between EF and DDIage, specifically c) executive function score and d) BRIEF GEC, for both the TBI and control groups.

Due to the non-normal distribution of our DDI measures, for visualisation purposes, Fig. 3 displays the relationship between our DDI measures and the indexes of EF where the log of the DDI variables are plotted, rather than the observed values.

Fig. 3

Scatter plots of the log of the DDI measures (top row DDIinv, bottom row DDIage) plotted against the EF measures (first column EF score, second column BRIEF).

DDI of the CEN

When investigating DDI calculated from a subgraph consisting of regions of the CEN, the difference of DDIinv. from the TBI group compared to the control group was not significantly greater than zero (W = 1146, p = .730). The difference of DDIage from the TBI group compared to the control group was not significantly greater than zero (W = 1302, p = .335). For the CEN, across DDIinv and DDIage, no tested association with EF was found to be significant, as seen in Fig. 4.

Fig. 4

Correlation between DDIinv. of the CEN, and a) EF and b) BRIEF, and correlation between DDIage of the CEN, and c) EF and d) BRIEF for both TBI and control groups.

Post-hoc exploratory analyses

We conducted post-hoc analyses to assess robustness of these findings. When the DDIinv was calculated using 500 random split halves from the ABIDE data, there was considerable agreement between DDI calculated from the first and second halves of the sample (DDIinv mean pearson’s r = .988, mean spearman’s rho = .981).

For both DDIinv and DDIage, we compared our DDI measure between controls and injury severity groups (mild and moderate/severe). Across the DDIinv and DDIage calculated for both the whole brain and CEN, no significant differences are reported. These results are seen in the Supplementary materials. Partial correlations between whole-brain DDIinv/DDIage and EF/BRIEF controlled for multiple factors that may have biased analyses. When controlling for age at MRI, correlation coefficients remained qualitativelty similar to those found previously. This was also true when simultaneously controlling for age at injury and interval between MRI and EF assessment. These results are seen in the supplementary materials.

Discussion

Previous research suggests a TBI during childhood can result in the deviation of the brain from the typical developmental trajectory (Dennis et al., 2017a,b; King et al., 2019; Mayer et al., 2015; Wilde et al., 2012a,b; Wu et al., 2010) and that these deviations may act as a marker of brain health, neurological disorders and cognitive functioning (Cole and Franke, 2017; Erus et al., 2015). We aimed to quantify. The current study utilises a modified version of Saggar et al. (2015) add-one-patient approach, termed the developmental divergence index, to calculate individual-level divergence from the typical SCN (estimated from a large paediatric dataset) for a cohort of patients who have experienced a TBI during childhood. This was a proxy measure of the level of divergence, with greater divergence hypothesised to be associated with poorer functional outcome. For the first time, the current study combined both this measure of divergence with a ‘sliding-window’ approach to generate developmentally-appropriate, age-matched, reference networks.

The current study found significant correlations between an index of divergence, calculated both against a general paediatric reference group but also an age-matched reference group, and executive functioning, measured with both performance and behavioural measures. These were in the expected direction; greater distance from a typical reference network was related to worse executive function skills and increased behavioural-problems related to poor EF. We found these relationships in both the whole sample and the subgroup of TBI patients only and not the control group. These results are in spite of the considerable heterogeneity in the neuropathology which occurs as a result of TBI (Dennis et al., 2017a,b), and the global, whole brain nature of the DDI metric. This may highlight the benefit of considering the broader impact of the injury and subsequent development beyond the individual regions.

The greater strength of association seen in the TBI group is somewhat unsurprising. Whilst in the patient group developmental divergence due to injury is likely to explain much of the variation in EF outcomes, in the control sample, it is likely that other individual differences explain a greater proportion of variance. The magnitude of these relationships between the DDI and EF are small. However, due to the limited sample size, estimating accurate point estimates of the correlation coefficients is difficult, as seen in the confidence intervals listed in Table 3. These wide confidence intervals also prohibited investigating whether the correlational relationships were significantly different between controls and TBI patients. Given the brain-behaviour relationships being seen in the whole sample, it is important to tease apart whether the DDI measure represents ‘normative’ development in the TBI group, rather than informing us how neuropathological effects (such as developmental divergence) are potentially disrupting the development of cognitive skills.

Table 3

Pearson’s correlation coefficients (r), 95 % bootstrapped confidence intervals and associated permutation-based p-values for each group and the sample as a whole.

DDI Measure	DV	Regions	TBI Patients				Controls				Whole Sample
			r	Lower CI	Upper CI	p	r	Lower CI	Upper CI	p	r	Lower CI	Upper CI	p
DDIinv	Age at Injurya	Whole Brain	.007	−.203	.168	.956	–	–	–	–	–	–	–	–
		CEN	.054	−.206	.261	.659	–	–	–	–	–	–	–	–
.	EFb	Whole Brain	−.319	−.589	.085	.024*	−.040	−.388	.322	.839	−.300	−.565	.081	.009*
		CEN	−.011	−.371	.241	.939	−.003	−.428	.337	.987	−.037	−.246	.187	.746
	BRIEFc	Whole Brain	.272	−.221	.623	.053	.421	.155	.631	.020*	.299	−.130	.635	.013*
		CEN	.111	−.089	.384	.418	−.004	−.261	.280	.981	.108	−.093	.336	.316
DDIage	Age at Injurya	Whole Brain	.041	−.198	.202	.730	–	–	–	–	–	–	–	–
		CEN	.093	−.137	.259	.426	–	–	–	–	–	–	–	–
	EFb	Whole Brain	−.330	−.604	.203	.021*	−.054	−.426	.304	.786	−.308	−.563	.082	.007*
		CEN	.077	−.178	.363	.586	−.027	−.564	.339	.893	.035	−.147	.217	.764
	BRIEFc	Whole Brain	.260	−.309	.648	.058	.309	.035	.622	.083	.277	−.144	.661	.021*
		CEN	−.023	−.276	.219	.878	.138	−.045	.426	.418	.018	−.150	.206	.865

Cases for correlation is n = 75.

Complete cases for correlation are n = 52 for TBI group and n = 28 for controls.

Complete cases for correlation are n = 52 for TBI group and n = 32 for controls. Upper and lower 95 % confidence intervals for the correlation coefficients are calculated using a bootstrap approach with 100 iterations, CIs which do not cross zero are highlighted in bold. p values are raw, uncorrected values calculated using a permutation approach with 5000 resamplings.

However, no significant differences were found between controls and patients within the TBI cohort in estimated DDI, across DDIinv. and DDIage for both whole brain and in the CEN. This was despite optimising our window-size to maximise between group differences (Supplementary materials – Appendix C). The sample of pTBI patients used for the current study was recruited across all injury severities, from mild to severe, with the majority of cases falling within mild injury. Whilst there is evidence for morphometric change due to injury across moderate to severe injury classifications (King et al., 2019) there is less evidence for this difference in mild injury cases (i.e. Ryan et al. (2017)). We therefore compared DDI metrics between injury severity groups in a post-hoc analyses, and yet no differences were found. Whilst there are no significant group differences in DDI, even at this very early stage post-injury, the DDI measure showed predictive validity with regard to executive functions. It is important to note the timings of both the MRI (<90 days post-injury) and neuropsychological assessment (24 months post-injury). The existing literature shows that neuroanatomical changes that occur post-injury persist over time (King et al., 2019). Given this, and the fact that we are still able to find these significant relationships (despite their weak magnitude) between relatively acute neuroanatomy and chronic functional outcome, one explanation is that these acute changes to the brain in response to injury, seemingly have a persistent effect which may guide the subsequent neurodevelopment required to subsume these executive functions. However, it is important to remember the evidence presented here is not causal in nature, but it does provide strong grounds upon which to further explore these relationships in independent cohorts. Overall, the current study shows that early imaging can assist in prognosis for cognition and therefore guide early intervention planning.

Cognitive-skills are particularly vulnerable to dysfunction due to damage during the period of skill-maturation (Ewing-Cobbs et al., 2004; Krasny-Pacini et al., 2017). Thus, the protracted development of EF (Diamond, 2013; Friedman et al., 2016; Perone et al., 2018) is likely to result in an extended window of vulnerability of EF to brain insult (Krasny-Pacini et al., 2017). Mechanistically, this vulnerability is likely due to damage within still-developing brain networks that subsume EF development (Khundrakpam et al., 2013). Essentially, a key principal is that, developmental processes happening at the time of insult are those which are the most vulnerable (Spencer-Smith and Anderson, 2009).

Structural covariance has an ongoing developmental trajectory throughout the neonatal period, childhood and adolescence (Alexander-Bloch, Raznahan, et al., 2013; Fan et al., 2011; Khundrakpam et al., 2016, 2013; Raznahan et al., 2011a,b; Váša et al., 2017), and structural covariance across association-cortex networks such as those supporting EF has a yet more protracted development (Khundrakpam et al., 2013). Thus, we investigated whether deviation of SC across regions of the CEN (which are commonly reported as supporting common EF activation in adolescence and childhood (Horowitz-Kraus and Holland, 2016; McKenna et al., 2017)), was related to executive dysfunction. The fronto-parietal regions included in the CEN are commonly affected by pTBI (King et al., 2019; Wilde et al., 2005), likely due to unique biomechanics of injury in the context of the paediatric brain (Pinto et al., 2012). Also, cross-sectional differences in cortical thickness of dorsolateral pre-frontal cortex have been found acutely post-injury (McCauley et al., 2012, 2010; Urban et al., 2017) with significant correlations between CT of frontal brain regions and BRIEF (Wilde et al., 2012a). Despite these findings, divergence from the age-appropriate structural covariance in the CEN was not associated with later EF. Overall, our findings support previous conclusions that the integrity of development in the entire brain is necessary for achieving age-appropriate, intact EF (Anderson et al., 2010), rather than early vulnerability due to specific damage to the networks that subsume EF development (Anderson, 2002). Taken together, these findings underscore the importance of considering metrics of connectivity when attempting to understand how brain insults impact on functional outcomes in a developmental context.

We used a composite measure of EF scores to explore structure-function relationships and this may contribute to the patterns of results reported. We adopted this approach to mitigate the relatively small sample size and the need to preserve statistical power. Thus, we were unable to investigate these skills with more granularity by examining discrete sub-components of EF. Such an approach would enable us to uncover whether regional / network deviations explain variance in specific EF impairments and future research should consider these more complex relationships. This is especially important given the variability in the age at which these different sub-domains of functioning (i.e. inhibitory control) come on-line during childhood (Miyake et al., 2000) and thus may differentially ‘react’ at different ages at which the injury occurs.

If we make the assumption that there are critical periods of vulnerability to the mechanical and pathological effects of injury then we might assume that greater divergence may be seen at one age versus another (Anderson et al., 2011; King et al., 2019). This may be due to the effects of injury differentially interacting with the myriad of developmental process that occur at different points throughout childhood brain development. Interestingly however, we found no linear relationship between age at injury and our proxy measure of brain perturbation. This is inconsistent with the idea of critical-periods of vulnerability, with no age at injury showing greater propensity to greater developmental divergence. Previous research investigating potential ‘age at injury’ effects post-TBI, do not primarily consider the magnitude of the perturbation the brain post injury (i.e. (Resch et al., 2019)). Thus, the current research opens up new opportunities in this area, offering a quantitative measure of brain perturbation (the DDI) by which we can investigate the individual and potentially interactive effects of both age at injury and magnitude of injury. This will better inform our understanding of critical-periods of vulnerability to TBI.

Generating an SCN allows the investigation of population-level covariance in neuroanatomy (Alexander-Bloch, Raznahan, et al., 2013). The individual contribution metric (proposed by Saggar et al. (2015)), enables an estimate of the distance of a patient from a group-level, reference SCN, to allow subject-level analyses. Previous studies show that greater divergence from the ‘typical’ SCN is related to worse neuropsychological performance (Saggar et al., 2015; Watson, 2016a).

In the context of the current study of pTBI, we ‘rebrand’ this metric as a measure of ‘developmental divergence’. This focus is primarily based upon our approach adopting a novel analytic framework whereby we use developmentally-appropriate control groups to calculate a reference network for the typically-developing SCN, using a sliding window approach. With the advent of large-scale, publicly-accessible neurodevelopmental studies such as ABIDE (Di Martino et al., 2014), ABCD (Casey et al., 2018) and HCP-development (Somerville et al., 2018) we are able to better understand the normative variation in brain development across age. The current study capitalised on this by calculating age-appropriate reference networks using MRI of typically-developing children from the ABIDE dataset.

The benefits of this are two-fold. Firstly, the variance of age within the window still allows us to better capture typical developmental variance within age bands, which here means that our reference groups from the ABIDE data captures variation due to individual differences in morphometry. By using discrete windows, which act similarly to age-bins, we also account for non-linearity in the changes to the SCN over time, as opposed fitting a continuous/linear reference trajectory. Previous studies have used a single control group to calculate Saggar et al.’s (2015) individual contribution metric, potentially conflating ‘normal’ differences due to discrepancies in age between the participant and the reference network with what is proposed to be ‘pathologically’-related divergence.

However, there is a limited number of cases at much younger ages in ABIDE. Thus, estimation of the DDI at these younger ages may be less reliable. A further limitation of the window-based approach is the small number of subjects with which each window was constructed (n = 26), given the size of the correlation matrix being estimated (68 × 68). The size of this window was selected empirically, based on maximising the between group-difference and the recommendations of Saggar et al. (2015). Future research could use a larger reference group to allow ‘denser’ age-windows to be generated with more subjects. However, this could result in the ‘mean’ network generated from the age-matched window being highly robust to the addition of new participants, and thus, based on the addition of the patient (AOP), the distance between and would be minimal, and the DDI measure is therefore likely to scale with the size of the reference group. This makes between study comparisons difficult.

We posited that deviation from a developmentally-appropriate reference group represents developmental-divergence and that in the context of a preceding brain injury, this reflects a negative perturbation or abnormality of expected brain development at a macroscopic-level. However, compensatory responses to brain injuries may also contribute to observed measurements of developmental divergence. The potential capability of the brain to experience adaptive or compensatory morphometric change, due to mechanisms such as neural plasticity, could potentially lead to restitution of function (Anderson et al., 2011; Bigler and Wilde, 2010). Therefore, one potential limitation of the DDI methodology is that it fails to disentangle change due to pathology and that which is compensatory and assists in recovery. Because MRI scans were conducted acutely (<90 days post-injury), divergence from typical morphometry, at this stage post-injury, is likely to be related to injury mechanisms, rather than recovery mechanisms. However, previous research observes both a persistent morphometric difference from controls, even at 10 years post-injury (Beauchamp et al., 2011), but also an ongoing neurodegenerative effect of injury (Keightley et al., 2014), typically related to worse cognitive performance (King et al., 2019). Therefore, we believe that the majority of variance in DDI is due to injury-related change. Future research may also investigate DDI pre- and post- neurorehabilitation, in order to investigate the role of divergence from typically-developing reference groups as a potential indicator of positive divergence supporting recovery of function. Differences in pre-processing steps used in our own experimental sample and that of the ABIDE reference group may influence the pattern of findings we observed. The ABIDE data was pre-processed using Freesurfer version 5.1 whilst our data was processed using the newer 6.0 release. Previous studies (and the Freesurfer developer community) recommend not comparing morphometric results between versions, with significant differences in measures being found for the same MRI scans (Chepkoech et al., 2016; Gronenschild et al., 2012). However, these differences will be systematic across all participants, in which case the DDI measure will comprise of a combination of systematic version-error and the ‘true’ divergence. Also, no direct comparisons have occurred between the morphometric measures calculated on different versions. The SCN networks were produced from the inter-correlations of these measures and then these SCNs are then compared, rather than the raw data. Future research may wish to consider this as a potential area of concern needing greater study.

Conclusions

We calculated individual-level divergence from the SCN (estimated from a large paediatric dataset) for a cohort of pTBI patients and found an association whereby greater divergence from the normative SCN was related to poorer executive functioning two years later. By investigating the CEN we took a neural-systems perspective to cognitive dysfunction, on the assumption that ‘damage’ to the network of regions supporting EF will relate to executive dysfunction (Anderson, 2002). However, the lack of correlation between CEN DDI and executive dysfunction in the TBI group highlights the nuanced role of immature networks subsuming neuropsychological functioning in childhood and that whole-brain integrity is required for age-appropriate EF abilities.

We propose that the DDI of the whole cortex may provide unique insights into the effects of brain injury on typical neurodevelopmental outcomes following early life brain injuries, and could be used in predictive models that seek to identify more accurately those children at greatest risk of long-term difficulties.

Declaration of Competing Interest

None.