Connectomic assessment of injury burden and longitudinal structural network alterations in moderate-to-severe traumatic brain injury.

Traumatic brain injury (TBI) is a major public health problem. Caused by external mechanical forces, a major characteristic of TBI is the shearing of axons across the white matter, which causes structural connectivity disruptions between brain regions. This diffuse injury leads to cognitive deficits, frequently requiring rehabilitation. Heterogeneity is another characteristic of TBI as severity and cognitive sequelae of the disease have a wide variation across patients, posing a big challenge for treatment. Thus, measures assessing network-wide structural connectivity disruptions in TBI are necessary to quantify injury burden of individuals, which would help in achieving personalized treatment, patient monitoring, and rehabilitation planning. Despite TBI being a disconnectivity syndrome, connectomic assessment of structural disconnectivity has been relatively limited. In this study, we propose a novel connectomic measure that we call network normality score (NNS) to capture the integrity of structural connectivity in TBI patients by leveraging two major characteristics of the disease: diffuseness of axonal injury and heterogeneity of the disease. Over a longitudinal cohort of moderate-to-severe TBI patients, we demonstrate that structural network topology of patients is more heterogeneous and significantly different than that of healthy controls at 3 months postinjury, where dissimilarity further increases up to 12 months. We also show that NNS captures injury burden as quantified by posttraumatic amnesia and that alterations in the structural brain network is not related to cognitive recovery. Finally, we compare NNS to major graph theory measures used in TBI literature and demonstrate the superiority of NNS in characterizing the disease.

INTRODUCTION

Traumatic brain injury (TBI) is a global public health problem with 69 million new cases estimated to occur worldwide each year (Dewan et al., 2018). Primarily caused by motor vehicle accidents, falls, and sports concussions, TBI has claimed more than 50,000 lives in the US alone in 2014 (Centers for Disease Control and Prevention, 2019), and frequently leads to long‐term disabilities (Humphreys et al., 2013).

A major characteristic of TBI is the shearing of axons across the white matter, induced by external mechanical forces. Diffuse axonal injury (DAI), as it is called, causes disruptions in the connectivity between brain regions throughout the network (Adams et al., 1982; Hayes et al., 2016), leading to cognitive deficits (Fagerholm et al., 2015) that often require rehabilitation for recovery (Chua et al., 2007). Another important characteristic of TBI is its heterogeneity in many aspects including cause, mechanism, and severity of injury, as well as recovery rate and burden of chronic symptoms (Maas, 2016; O'Brien et al., 2020). In treatment and rehabilitation planning, this heterogeneity poses a big challenge that makes subject specific approaches necessary (Huie et al., 2020; Stocchetti et al., 2017). Despite the recent increase in connectomics research focusing on subject specific analysis (Chamberland et al., 2021; Jolly et al., 2020), neuroimaging biomarkers for individualized diagnosis have been very limited for TBI. Considering these two characteristics of TBI, network level analysis of connectivity disruptions in TBI is essential to provide subject specific measures quantifying injury burden of individuals (Attye et al., 2021), which would help in achieving personalized treatment, patient monitoring, and informing the patient and caregivers regarding the potential long‐term progression of the disease (Huie et al., 2020; Wilson et al., 2017).

Advancements in neuroimaging within the last decades have enabled analysis of connectivity disruptions in TBI with modalities such as functional (Hillary et al., 2014; Olsen et al., 2014) and structural MRI (Hutchinson et al., 2018; Levine et al., 2008; Xiong et al., 2014). Diffusion MRI (dMRI), a structural MRI method measuring the diffusion of water molecules in the tissue, has especially been promising in the analysis of TBI as it is shown to be sensitive to axonal injury at a microstructural level, that is not captured well in conventional MRI (Irimia et al., 2014; Niogi et al., 2008). Most of the dMRI based studies investigate axonal injury either locally in isolated brain regions (Singh et al., 2010) or across certain white matter tracts (Wang et al., 2011), by using dMRI measures such as fractional anisotropy, mean diffusivity, or radial diffusivity (Hulkower et al., 2013; Irimia et al., 2014, Vijayakumari et al., 2021). Analyses involving such microstructural measures, however, fall short in capturing the impact of TBI on overall network topology.

Analysis of structural connectomes, that is, connectivity maps derived from dMRI data quantifying connections between brain regions, enables evaluation of the brain as a network (Sporns et al., 2005). Despite TBI being considered as a “disconnection syndrome” due to damaged structural pathways connecting brain regions (Hayes et al., 2016), analysis of structural connectivity disruptions evaluated on connectomes and longitudinal change in network organization is surprisingly scarce (Imms et al., 2019; Irimia et al., 2014; Kim et al., 2014). The majority of studies investigating structural connectivity in TBI utilize graph theoretical measures, reporting increase in shortest path length (Kim et al., 2014) and small‐worldness (Yuan et al., 2015), and decrease in global efficiency, clustering coefficient (Raizman et al., 2020), betweenness centrality, and eigenvector centrality (Fagerholm et al., 2015). Changes in such graph theoretical measures are also explored for subnetworks of the brain, such as working memory and reasoning networks, demonstrating significant differences in patients along with a correlation with cognitive scores (Jolly et al., 2020). While such measures provide insights into the mechanisms of change of the brain's network structure in TBI, each measure captures a specific aspect of connectivity alteration in the network, which are limited in capturing the overall topological change representing injury burden (Caeyenberghs et al., 2014; Dennis et al., 2017; Raizman et al., 2020). As they are mathematical constructs that are defined for networks at large without any special consideration for brains, interpretation of ensuing results poses further challenges. Additionally, in the absence of a hypothesis that defines the nature of TBI induced change in network topology, it is common to explore a large set of graph theoretical measures that are available in the literature to find those that would demonstrate statistical significance with the data. This exploratory approach, however, suffers from multiple comparison issues (Poldrack et al., 2017), affecting TBI studies more than other neuroscientific research due to small sample sizes in the domain. Hypothesis driven studies that suggest markers for TBI by taking the characteristics of the disease into account, on the other hand, are scarce (Kuceyeski et al., 2019; Solmaz et al., 2017).

Of special relevance to the current study, Irimia et al. (2014) have explored longitudinal changes of connectivity in three TBI patients within 6 months post injury at interhemispheric level as well as individual connections. Being an early example of connectomic analysis of TBI, this study is limited in statistical analysis due to the limited sample size. More recently, Kuceyeski et al. (2019) have reported increased network segregation in structural and reduced integration in functional connectivity of TBI patients. Being a hypothesis driven longitudinal study, this study explored only mild TBI patients. Solmaz et al. (2017) proposed a specialized network level score quantifying connectivity disruptions as a weighted sum of frequently damaged edges across patients for capturing injury burden on moderate‐to‐severe TBI, which they showed to correlate with injury severity. This study, however, was limited in being cross‐sectional and edge centric. Overall, longitudinal analysis of network level change in moderate‐to‐severe TBI assessed with a subject specific measure leveraging the characteristics of the disease is still lacking.

In this article, we present a longitudinal analysis of brain connectivity changes in moderate‐to‐severe TBI patients using a novel measure that we call network normality score (NNS). NNS is designed to capture the integrity of structural connectivity in patients by leveraging two major characteristics of the disease, that are, diffuseness of the injury and the heterogeneity of the disease. Diffuseness of the injury can be best captured by a connectome‐level measure that is sensitive to the global effects of local connectivity disruptions. Heterogeneity of the disease, on the other hand, can be best captured by a normative measure that compares each patient with a reference healthy control sample. Taking a graph matching based approach by extending our prior connectomic measures (Osmanlıoğlu et al., 2018; Osmanlıoğlu et al., 2019), we define NNS as the overall network similarity of moderate‐to‐severe TBI patients relative to a healthy control sample. We hypothesize that NNS captures the injury burden of individuals with TBI, which we test by calculating correlation between NNS and posttraumatic amnesia scores of patients. We evaluate our measure on a cohort of 34 patients with moderate‐to‐severe TBI, who underwent dMRI and cognitive assessment at 3, 6, and 12 months postinjury, as well as 35 age‐ and sex‐matched healthy controls. In our analysis, we investigate cross‐sectional and longitudinal relationships between the NNS and injury severity, as well as cognitive outcome. We also investigate longitudinal changes in network topology of patients relative to controls as quantified by NNS, and evaluate its relationship with the change in cognitive scores over time. Finally, we compare NNS with standard graph theoretical measures that are commonly reported in TBI literature, in their relationship with injury severity and cognitive outcome.

MATERIALS AND METHODS

Participants

The data used in this study was acquired as part of a larger project investigating the neuroimaging correlates of functional recovery after diffuse TBI (PI: JJK). All participants provided informed consent directly or via a legally authorized representative. Study procedures were approved and overseen by the Institutional Review Board at the Moss Rehabilitation Research Institute, Elkins Park, Pennsylvania, and the University of Pennsylvania. The cohort investigated in this study consists of 40 participants with moderate‐to‐severe TBI and 35 healthy controls (HC) (Solmaz et al., 2017). Inclusion criteria for TBI participants were being in the age range 18 to 64 and diagnosis of non‐penetrating moderate‐to‐severe TBI, indicated by at least one of the following: (i) Glasgow Coma Scale score <13 in the emergency department (ED not due to sedation, paralysis, or intoxication), (ii) documented loss of consciousness for more than 12 h, (iii) Prospectively documented PTA >24 h. Exclusion criteria for TBI participants were (i) history of prior TBI, CNS disease, seizure disorder, schizophrenia, or bipolar disorder, (ii) history of long‐term abuse of alcohol or psychostimulants that could have resulted in neurologic sequelae, (iii) pregnancy, (iv) inability to complete MRI scanning due to ferromagnetic implants, claustrophobia, or restlessness, (v) nonfluency in English; or (vi) a level of disability preventing completion of testing and scanning by 3 months postinjury. TBI participants with total estimated volume of focal intraparenchymal lesions larger than 5 cm3 for subcortical lesions and larger than 50 cm3 for cortical lesions were also excluded to ensure that the TBI was predominantly diffuse. Healthy controls recruited were comparable in age, sex, and education to TBI subjects. Exclusion criteria for HCs were the same with TBI participants with the addition of exclusion for any history of TBI resulting in alteration or loss of consciousness.

Cognitive assessment and dMRI scans were obtained for HCs once and for patients three times at approximately 3, 6, and 12 months postinjury. Imaging data was not available for some of the patients at certain time points due to either the patient not attending a follow up session or the data being removed from the dataset because of MRI quality issues such as segmentation problems arising from lesions in the brain. Imaging quality assessment (QA) was undertaken by manual inspection of every volume of the dMRI scans for diffusion artifacts including motion‐induced signal drop‐out and slice‐wise intensity artifacts (Soares et al., 2013). If more than 33% of the unweighted volumes (i.e., 2 of b = 0 volumes) or more than 10% of the weighted volumes (i.e., 3 of b = 1000 volumes) were affected by artifacts, the entire scan was excluded; otherwise, any artifact‐affected volumes were removed from the data for further analysis. In our analysis, we removed six patients from the dataset whose imaging data failed the QA at 3 months postinjury, leaving 34 patients (12 f) to be analyzed for the study. Among these patients, 27 (10 f) had dMRI data available at 6 and 12 months. We note that dMRI data of only 22 (8 f) of the patients had passed the imaging QA at all three time points. In order to increase the power of the analysis, we used all patient data available at follow up sessions rather than doing the analysis with the patients that have data at all time points. Demographics of the participants are detailed in Table 1.

TABLE 1

Demographics of the moderate‐to‐severe TBI dataset with healthy controls.

		Healthy controls		Patients
		Male	Female	Male	Female
Count		26	9	22	12
Age	Avg.	36.7	30.0	35.5	34.0
	SD	9.4	10.8	14.7	15.6
PTA	Avg.	NA	NA	31	21.1
	SD	NA	NA	23.6	18.6
GCS a	Avg.	NA	NA	8.73	12.2
	SD	NA	NA	4.68	2.7

GCS score was not assessed for 10 patients that were sedated or intubated at admission to ED. Average GCS score is reported for the remaining 24 patients.

Data acquisition, preprocessing, and connectome construction

Structural MRI scans were acquired on a Siemens 3 T TrioTim whole‐body scanner with an 8‐channel array head coil (single‐shot, spin‐echo sequence, TR/TE = 6500/84 ms, b = 1000 s/mm2, 30 directions, six unweighted [b = 0] volumes, flip angle = 90°, resolution = 2.2 × 2.2 × 2.2 mm). High‐resolution T1‐weighted anatomic images were also obtained using a 3D MPRAGE imaging sequence with TR = 1620 ms, TI = 950 ms, TE = 3 ms, flip angle = 15°, 160 contiguous slices of 1 mm thickness, FOV = 192 × 256 mm2, 1NEX, resolution = 1 × 1 × 1 mm. T1 images were preprocessed using the FreeSurfer 5.3.0 recon‐all pipeline (http://surfer.nmr.mgh.harvard.edu) (Fischl, 2012) and registered to the FA using rigid followed by deformable SyN registration in ANTs (Avants et al., 2008) with the deformation constrained to the anterior–posterior direction to correct for the EPI distortions in the dMRI. 86 regions of interests from Desikan atlas (Desikan et al., 2006) were extracted to represent the nodes of the structural network. Five‐tissue‐type images for anatomically constrained tractography (ACT) (Smith et al., 2012) were created from Freesurfer outputs. 500 seeds for tractography were placed at random inside each voxel of the mask of the gray‐matter white‐matter interface (GMWMI). Single‐shell, multi‐tissue constrained spherical deconvolution (Dhollander et al., 2016; Jeurissen et al., 2014; Tournier et al., 2007) was performed in mrtrix3 (Tournier et al., 2019) to fit a fiber orientation distribution (FOD) at every voxel in the brain. Probabilistic tractography was performed using the iFOD2 algorithm (Tournier et al., 2010) with angle curvature threshold of 60°, step size of 1 mm, and minimum and maximum length thresholds of 25 and 250 mm, respectively. Connectomes were then generated as an 86 × 86 adjacency matrix of weighted connectivity values, where each element represents the number of streamlines between regions. Each connectome was subsequently normalized by the GMWMI volume of the individual.

Behavioral and cognitive measures

TBI patients underwent behavioral assessment at each time point to yield multiple outcome measures. Duration of posttraumatic amnesia (PTA), calculated as the number of days between the TBI and the first of two occasions within 72 h that the patient was fully oriented, was used as a sensitive behavioral index of the injury severity (Benson et al., 2007; Povlishock & Katz, 2005). Full orientation was defined as a score above 25 on the Orientation Log (Jackson et al., 1998), or documentation of consistent orientation for 72 h in the medical record.

Cognitive outcome was measured using a test battery assessing three neuropsychological domains: information processing speed (PS), verbal learning (VL), and executive functioning (EF). To create these domain scores, all cognitive test scores were transformed to T‐score units based on available normative data. We used T‐score transformed Processing Speed Index from the Wechsler Adult Intelligence Scale‐IV (Wechsler, 2008a) to represent the PS domain. To operationalize the VL domain, the T‐score for immediate recall trials I‐V was used from The Rey Auditory‐Verbal Learning Test (Rey, 1958). For the EF domain, we built a composite score to reduce type I error and increase signal‐to‐noise ratio, which was defined as an average of the T‐score transformed scores obtained from the following five tests: Controlled Oral Word Association Test (Benton et al., 1994), Trail Making Test‐Part B (Reitan & Wolfson, 1985), the Color Word section of the Color‐Word Interference Test from the Delis‐Kaplan Executive Function System (Delis, 2001), and Digits Backward and Letter‐Number Sequencing subtests from the Wechsler Memory Scale IV (Wechsler, 2008b). Further details on building domain scores can be found in (Rabinowitz et al., 2018).

Network normality score

In order to evaluate change of brain's network organization in TBI patients over time, we consider graph matching (Foggia et al., 2014; Osmanlıoğlu, 2016) to quantify connectomic similarity as it accounts for changes in the overall topology of the network rather than focusing on local changes in individual connections. Previously, we have successfully applied graph matching in deriving similarity between connectomes for quantifying injury severity in TBI patients (Osmanlıoğlu et al., 2018; Shen et al., 2020), evaluating subject‐wise structure–function correspondence (Osmanlıoğlu et al., 2019), and investigating connectomic stability within and across subjects (Osmanlıoğlu et al., 2020). In this study, we extend our previous approach by adopting a different use of graph matching to provide a normative connectomic similarity measure.

A graph matching based measure to quantify connectomic similarity

Here, we first provide a brief overview of graph matching and its use in connectomics, and refer the reader to (Osmanlıoğlu et al., 2019; 2020) for further information where we first devised the approach. Given two graphs A and B that are deemed to have a similar topology, the aim of graph matching is to find the optimal mapping between the two graphs by assigning each node of A to a node of B that structurally resembles it the most. Given a cost function c: A ⨉ B → ℝ determining the cost of assigning each node in A to a corresponding node in B, graph matching can be formulated as a combinatorial optimization problem where the aim is to calculate a one‐to‐one mapping f: A → B between the nodes of A and B by minimizing the objective function . On connectomes, we regarded the cost function c as the Euclidean distance between the k‐dimensional feature vectors of nodes encoding their connectivity signature relative to other nodes in a parcellation with k ROIs. We obtained the desired mapping by solving the optimization problem using the Hungarian algorithm (Kuhn, 1955).

Since brain structure has commonalities across healthy people with a relatively small variation compared with patients, and the parcellation that yielded graph representations of brains are the same across subjects, the resulting mapping would match nodes of A with their corresponding nodes in B among healthy subjects (i.e., the matching nodes should correspond to the same ROI). We call such a mapping a correct match. On the other hand, the brain structure of patients would have a larger variation and dissimilarities relative to healthy controls as well as among each other due to disease induced alterations, resulting in connectivity patterns of nodes that vary too much between the two graphs. This would lead the mapping to have incorrect matching of some of the nodes where nodes in A will be assigned to nodes in B that are not their counterparts. Consequently, we regarded network similarity (NS) as the percentage of correct matches relative to total number of nodes (Figure 1a), with larger values indicating higher similarity. Using this graph matching based measure in quantifying network similarity allows capturing the similarity of overall network organization since matching between the nodes are obtained through the solution to an optimization problem.

FIGURE 1

Network normality score (NNS) quantifying similarity of structural network organization of a subject's brain relative to a sample. (a) Taking two connectomes representing the structural connectivity of two subjects as input, the similarity between their graph representation is calculated using graph matching, yielding a binary matching matrix. Similarity between the connectomes is determined as the proportion of nodes which were correctly matched. Using graph matching (GM) as the measure of network similarity, we calculate network normality score of (b) each healthy control relative to the rest of the healthy controls (NNSH), (c) each patient at a certain time point relative to healthy control sample (NNSPH), and (d) each patient relative to the rest of the patients at the same time point (NNSP).

Normative connectomic similarity: similarity of a subject relative to a sample

Having defined NS as the similarity measure between two connectomes, we next define the network normality score (NNS) as a normative measure consisting of the mean NS of the subject relative to the reference sample (Figure 1b–d). Taking healthy controls as the reference, we first calculated network normality score among them to provide a basis for evaluation (Figure 1b). We then calculated similarity of patients at a certain time point (such as 3 months) relative to the healthy (Figure 1c), to quantify trauma induced network alterations in TBI patients. In order to evaluate heterogeneity and the course of relative changes in network topology among patients, we calculated a third anomaly score quantifying similarity of patients relative to the rest of the patients within the same time point (Figure 1d). For the sake of clarity, we refer to these three scores as NNSH, NNSPH, and NNSP in the rest of the paper, where subscripts denote NNS among the healthy, NNS of patients relative to healthy, and NNS among the patients, respectively.

Statistical analysis

Group level analysis: In order to evaluate group differences in the structural network organization cross‐sectionally, we ran the Mann–Whitney U test between NNSPH (or NNSP) and NNSH, and the Wilcoxon signed‐rank test between NNSPH and NNSP. Non‐parametric tests were preferred in calculating group difference since NNS scores demonstrated a non‐Gaussian distribution per D'Agostino's K‐squared normality test (NNSH: p < 10−6, NNSPH: p < .01 for patients at 3 and 6 months, p = .07 for patients at 12 months, NNSP: p = .08, .006, .17 for patients at 3, 6, and 12 months, respectively). We quantified the amount of change in scores using the following effect size formula: where are the sample size of groups, and , , and , and U is the test statistic for the Mann–Whitney U test, whereas for the Wilcoxon signed‐rank test. Effect size is regarded as small if |ES| ≥ 0.1, medium if |ES| ≥ 0.3, and large if |ES| ≥ 0.5.

Cross‐sectional linear model analysis: In devising a subject specific measure to quantify injury severity with a potential to serve as a biomarker, it is necessary to evaluate its relationship with cognitive scores of individuals cross‐sectionally at certain time points. To achieve this, we utilized a linear model (LM) format that controls for age and sex as follows: where diseaseRelatedScore is replaced by one of the cognitive scores or PTA, while tp indicates one of 3, 6, or 12 months time points. Analyses were done in R using the nlme package (Pinheiro et al., 2021).

Longitudinal linear mixed effect model analysis of NNS: In order to investigate whether the network organization of patients demonstrates a linear change over time when considered altogether, we evaluated the longitudinal change in their network normality scores (NNSPH and NNSP are evaluated separately). Since imaging data was not available at all time points for some subjects, we used linear mixed effects (LMEM) analysis with the following model: where we estimated NNS as a linear function of the fixed variables days since injury (DSI), PTA, age, and sex, along with the random intercept. Analyses were done in R using the lme4 (Bates et al., 2015) and lmerTest (Kuznetsova et al., 2017) packages. In our LM and LMEM analysis, we scaled the values of variables. Thus the estimated values of independent variables (e.g., DSI, PTA, age, etc. in Equation 2) can be interpreted as their correlations with the dependent variable (e.g., NNS in Equation 2).

Analysis of the trajectory of change: In order to evaluate whether there exists a relationship between the cognition and NNS in how they change over time, we calculated the rate of change as the slope of the line connecting measurements between two time points for each score type. We then used the following linear model that also accounts for age, sex, and PTA as baselines to determine the relationship:

Longitudinal linear mixed effect model analysis of cognition relative to NNS: We used a linear mixed effects model to longitudinally evaluate the relationship between cognitive scores and network similarity score of patients relative to healthy controls, days since injury, age, and sex across three time points using the following formulation:

Standard graph theoretical measures

In order to further highlight the efficacy of the proposed score in characterizing TBI and evaluating whether its global optimization based approach better captures TBI, we evaluated our moderate‐to‐severe TBI cohort using standard graph theoretical measures that are cited in the TBI literature. We considered node betweenness centrality, eigenvector centrality, clustering coefficient, small worldness, characteristic path length, global efficiency, and modularity, as these are the measures commonly used in the prior work. We used the Python implementation (bctpy, version 0.5.2, https://pypi.org/project/bctpy/) of Brain Connectivity Toolbox (Rubinov & Sporns, 2010) to calculate the measures over the connectomes. The statistical analysis for NNS was repeated for each of these graph theory measures individually (see Supporting Information S5 for further details on graph theory measures and their analysis).

RESULTS

Group level analysis of network similarity between patients and controls

In order to evaluate whether the proposed measure captures structural connectivity alterations, we performed a group‐level analysis between patients and healthy controls (Figure 2). We observed significantly lower network similarity scores for patients (NNSPH) compared with healthy controls (NNSH) at 3 months (ES = 0.42, p < 10−3), 6 months (ES = 0.41, p = 10−3), and 12 months (ES = 0.59, p < 10−4). This result shows that the NNS captures TBI induced alterations of the network topology to distinguish structural connectivity of patients from that of healthy controls up to 12 months postinjury.

FIGURE 2

Group level analysis of network normality score across patients and controls. We evaluated network normality scores (NNS) of patients relative to healthy controls (NNSPH), NNS among patients within the same time point (NNSP), and NNS among healthy controls (NNSH). We observed that network topology of patients is significantly dissimilar to that of the healthy (NNSPH < NNSH, purple lines at the top), showing that trauma induced injury introduced alterations across the network. We also observed that NNSH > NNSP with statistical significance (green lines at the top), highlighting a larger variance of network topologies among patients than controls. We then observed that NNSPH > NNSP (black lines at the top), indicating that patients resemble the healthy more than they resemble other patients. These results further show that the heterogeneity of the disease is captured at the structural brain network topology of patients (Note that lines at the top between pairs of sample groups show effect size for significant group differences with p < .05, results are FDR corrected).

We then investigated whether NNS captures the heterogeneity of the disease at the network level. We observed that patients had significantly lower within‐group network similarity scores (NNSP) compared with healthy controls (NNSH) at 3 months (ES = 0.72, p < 10−6), 6 months (ES = 0.76, p < 10−6), and 12 months (ES = 0.82, p < 10−6). This result underlines a higher heterogeneity in structural network topology among patients than that among healthy controls, indicating that the injury affecting each patient differently leads to a unique network organization. We also observed a significant group difference between NNSPH and NNSP at 3 months (ES = 0.49, p < 10−4), 6 months (ES = 0.54, p < 10−4), and 12 months (ES = 0.60, p < 10−4), which indicate that network structures of patients resemble that of healthy controls more than they resemble that of other patients.

Relationship between network similarity and injury severity

A significant negative association between PTA and NNSPH was observed (see Equation (1) for the LM) at 3 (p = .016, estNNS = −0.51), 6 (p = .016, estNNS = −0.48) and 12 months (p = .016, estNNS = −0.52) (Figure 3), while no significant association was observed for age and sex (see Table 2). This result indicates that more severely injured patients have lower network similarity in reference to healthy controls.

FIGURE 3

Relationship between network normality score and injury severity. Evaluating whether injury severity (PTA) can be described cross sectionally as a function of NNSPH, age, and sex using an LM, and we observed a significant relationship between NNSPH and PTA at 3, 6, and 12 months postinjury (p values are FDR corrected). This result indicates that trauma induced alterations at network topology captures injury severity.

TABLE 2

Results of fitting a linear model to evaluate the relationship between injury severity (PTA) and NNSPH, age, and sex.

Time point	Adj.R 2	p NNS	estNNS	p age	estage	p sex	estsex
3 Months	0.202	.016	−0.509	.810	−0.083	.270	0.568
6 Months	0.223	.016	−0.480	.810	0.045	.454	0.336
12 Months	0.176	.016	−0.521	.810	−0.178	.454	0.276

Note: We note that since the scores were scaled for the LM analysis, the estimated values provided in the table (columns labeled with “est”) indicate correlation of corresponding variables with PTA (see Equation 1 for LM, p values are FDR corrected for each variable across three models).

Change in network normality score over time

Group level analysis of network similarity scores shown in Figure 2 demonstrated an increase in effect size between patients and controls from 3 to 12 months, suggesting that the structural connectivity in patients becomes more unlike healthy controls over time. A further longitudinal analysis of NNSPH using LMEM (Equation 2) showed that the similarity score is a function of days since injury (DSI), PTA, and age (p DSI < 10−3, p PTA = .004, p age = .019, p sex = .581, Adj. R 2 = 0.316), with a negative association (estDSI = −0.144, estPTA = −0.429, estage = −0.33) (see Table S1a for further details). This result indicates that patients become significantly unlike healthy controls in their structural network connectivity as time progresses postinjury up to 12 months (Figure 4).

FIGURE 4

Analysis of change in network normality score of patients. Using an LMEM, we evaluated the change in NNSPH score as a function of days since injury, PTA, age, and sex, observing a significant decline in NNSPH with time. This result indicates that the structural network topology of the patients becomes unlike that of healthy controls over time.

Repeating the same analysis for network similarity among patients, we observed a significant decrease of NNSP over time (p DSI < 10−4, p PTA = .007, p age = .066, p sex = 0.662, Adj. R 2 = 0.298) with a slope of −0.306 for DSI, indicating a steeper decline when compared with change in NNSPH with a slope of −0.144 (see Table S1b for further details). This result indicates that although patients deviate from “normalcy” as defined by the network topology of the healthy, they do not converge to an alternate normal that would be common among patients either.

Relationship between the network similarity score and cognitive scores

We next investigated whether NNSPH captures information regarding cognitive function. Before evaluating the relationship between NNSPH and cognition, we first did a group level and LMEM analysis of cognitive scores to evaluate their change over time and their relationship with PTA. We observed that the patients perform significantly lower than controls at 3 months for each cognitive score type (Figure S2, top) and that their performance in each category improves over time significantly to reach the level of healthy controls at 12 months (Figure S2, bottom). We also observed a significant negative correlation between each cognitive score and PTA, with verbal learning (VL) having a marginal p value (see Supplementary Information  S2 for further details). These results show the presence of cognitive recovery in patients up to 12 months postinjury and demonstrate that cognitive performance is related to injury severity.

Observing a disparity between cognitive recovery and increasingly abnormal network topology in patients, we evaluated whether there exists a meaningful relationship between the two virtually diametrical trends (Figure 5). Calculating the relationship between rates of changes between consecutive time points for the NNSPH and cognitive scores separately using the linear model (Equation 3), we observed no significant relationship at any of the time intervals (i.e., 3–6 months, 3–12 months, or 6–12 months, p > .05 for all variables. See Supplementary Information S4 for details), indicating the lack of an association between cognitive recovery and change in structural connectivity organization of patients.

FIGURE 5

Change in NNS and cognitive scores across time. Plotting individual trajectories of change in NNSPH and cognitive scores for each patient, we observed a steady increase (green lines) for cognitive scores in most cases indicating recovery. In contrast, we observed several cases of decrease (red lines) for NNSPH indicating deviation from normalcy in terms of network topology. Calculating the correlation between the slopes of lines in NNSPH with the slopes of lines in each cognitive score separately, we did not observe any significant relationship. This result indicates that the rate of change in NNS is not associated with cognitive recovery.

Despite the lack of a significant relationship between the rate of change in cognitive scores and NNSPH, we evaluated whether there exists a relationship between the actual scores. Using an LMEM (Equation 4), we observed that executive function (EF) and processing speed (PS) are significantly and positively related with NNSPH and DSI (EF: p NNS < 10−3, p DSI < 10−4, R 2 m = 0.206, PS: p NNS = .006, p DSI < 10−4, R 2 m = 0.226) while verbal learning did not reveal any significant relationship with NNSPH (p NNS = .086, pDSI = 10−4, R 2 m = 0.141) (see Table S3 for further details). The positive correlation between NNSPH and EF and PS indicates that patients with structural connectivity more similar to healthy controls demonstrated better cognitive function.

Evaluation of the cohort with standard graph theory measures

In our analysis of graph theory measures, we first evaluated the association between NNSPH and graph theory measures longitudinally and cross‐sectionally, and observed no significant relationship (see Tables S5a and S5b, p values are FDR corrected for multiple comparison correction, see Supplementary Information S5 for a detailed explanation of the analysis). We then evaluated the association between graph theory measures and PTA using LM (see Equation S5b) showed statistical significance only for node betweenness centrality at 6 months (p NBC = .002, p age = .962, p sex = .984, Adj. R 2 = 0.504) (see Table S5c). Finally, we evaluated the association between cognitive scores and graph theory measures using a LMEM analysis (Equation S5c). After FDR correction, no significant association was observed (see Table S5d).

DISCUSSION

Traumatic brain injury is considered a disconnectivity syndrome (Hayes et al., 2016) due to the diffuse injury of axons across the brain tissue, leading to structural connectivity disruptions among brain regions. While local microstructural changes in the brain (Hutchinson et al., 2018; Levine et al., 2008; Xiong et al., 2014) as well as functional connectivity alterations (Hillary et al., 2014; Olsen et al., 2014) are well studied in TBI, literature focusing on the structural connectivity changes in the brain has been very limited (Imms et al., 2019). This small body of work has two main limitations: First, most of these studies utilize standard graph theoretical measures in their analysis, which are limited in capturing the diffuse characteristics of the injury. Second, although cross‐sectional studies are abundant, longitudinal analysis of structural changes in the brain's network topology and its relationship with cognitive function of TBI patients are scarce. In this study, we proposed a novel measure called Network Normality Score (NNS) that is tailored to capture the two established characteristics of TBI, the diffuseness of the injury (Adams et al., 1982) and the heterogeneity of the disease (O'Brien et al., 2020). In a moderate‐to‐severe TBI cohort, we demonstrated that the NNS captures injury‐induced structural connectivity alterations by quantifying the connectivity differences at the network level. This highlighted a significantly different network topology among patients relative to healthy controls. Our results also show that the heterogeneity of the disease is observable in the network topology of the patients as quantified by the NNS. We further observed that the network structure of the patients becomes more unlike that of healthy controls over time, despite cognitive recovery over the same interval. As we did not observe any significant relationship between the change in cognitive scores and the change in network similarity of patients over time, these results highlight a mismatch between structural change and cognitive recovery. Finally, we demonstrated that the NNS captures characteristics of TBI that are not captured by standard graph theory measures as there was no significant association between the NNS and any of the graph theory measures. We also observed that only node betweenness centrality demonstrated a significant association with injury burden at 6 months, and none of the measures showed a significant association with cognitive scores, as the results did not survive multiple comparisons correction. Overall, our results point to a new direction of research in the analysis of structural network alterations in TBI, involving similarity measures that are designed to capture the characteristics of the disease such as heterogeneity and diffuse injury.

Overall network similarity of TBI patients relative to the healthy, captures injury induced alterations in the structural connectivity

The negative correlations between PTA and NNS indicate (Figure 3) that patients with more severe brain injuries (high PTA score) have network topologies that are less like healthy controls (low network similarity score). When considered with the group level differences of network topologies between patients and controls (Figure 2), these results highlight the efficacy of the NNS in capturing trauma induced alterations.

The direct relationship between diffuse axonal injury and the disruptions in structural connectivity among brain regions underlines the potential of a network topological analysis in quantifying injury burden of TBI patients. Interestingly, the number of studies examining this relationship longitudinally on moderate‐to‐severe TBI patients is not too many (Caeyenberghs et al., 2014; Raizman et al., 2020). We were able to identify four major studies that considered graph theoretical measures to evaluate network abnormalities of patients and evaluated their relationship with injury severity, two of which reported a lack of a significant relationship in moderate‐to‐severe adult (Caeyenberghs et al., 2014) and pediatric (Dennis et al., 2017) TBI patients. Although the third study reported a positive correlation for node strength and global efficiency scores, it is worth noting that their cohort included mild TBI patients as well as moderate‐to‐severe TBI patients (Raizman et al., 2020). Hypothesis driven longitudinal study of Kuceyeski et al. (2019) have reported increased network segregation in structural and reduced integration in functional connectivity of TBI patients, which they evaluated as a function of characteristic shortest path length. They also demonstrated a positive relationship between the change in structural and functional network topology and cognitive recovery. This study also was conducted on mild TBI patients. A recent study by our group proposed the Disruption Index of the Structural Connectome (DISC) as a specialized network level score for capturing injury burden on TBI, which demonstrated a significant correlation with injury severity of patients (Solmaz et al., 2017). However, this study was limited in being cross‐sectional and the connectivity disruptions being quantified on the basis of edges, rather than at network level.

The lack of significant associations of graph theory measures with PTA and cognitive scores (except for node betweenness centrality at 6 months with PTA) along with lack of a significant relationship between the NNS and any of those graph theory measures, indicate the novelty and superiority of our measure over standard graph theory measures in characterizing TBI. We note that standard graph theory measures are mathematical constructs that are designed to evaluate any graph structure such as social networks or airline route maps, without any specific consideration for brain networks. In the absence of a hypothesis on which measure to use as a biomarker, exploratory analysis that investigates several graph theory measures becomes inevitable. This, however, reduces statistical power of the study due to multiple comparisons correction, which is already limited in TBI studies due to small sample sizes. Interpretation of ensuing results is a further challenge due to measures not being disease specific. Designed specifically to capture well known characteristics of TBI, on the other hand, our proposed measure has two major strengths over standard graph theory measures. First, it focuses on leveraging the diffuse characteristic of the injury by taking a graph matching approach. Since graph matching quantifies similarity through solving an optimization problem, it considers connectivity differences across the network altogether, rather than summarizing connectivity differences on the basis of individual edges. Second, it is a normative score that is calculated relative to healthy controls that leverage the heterogeneity of the disease.

Heterogeneity of the disease is observable in the structural connectivity among brain regions

A major characteristic of TBI is its heterogeneity in various aspects including the cause of initial injury (e.g., fall or motor accident), mechanism (e.g., direct impact or acceleration/deceleration), pathology (e.g., focal and/or diffuse axonal injury), severity (e.g., mild, moderate, or severe), ensuing cognitive deficits, and treatment of the disease (Iraji et al., 2016; Maas, 2016) as well as outcomes in cognitive recovery (Felmingham et al., 2004; O'Brien et al., 2020). Lower NNS of patients relative to controls show that network topology of TBI patients differs from the healthy control population at varying degrees (Figure 2). Network similarity among patients being even lower than their similarity relative to healthy controls further supports the previous result, highlighting that injury affects each patient in different ways, potentially due to heterogeneity of the disease in its etiology, mechanism, and severity. In combination, these results demonstrate that the heterogeneity of TBI is also observable at structural brain connectivity of patients.

Revisiting structural plasticity in TBI

Diffuse axonal injury is one of the major characteristics of TBI, which causes disruptions in the connectivity between brain regions (Adams et al., 1982), leading to cognitive deficits especially in moderate‐to‐severe cases (Scheid et al., 2006). Rehabilitation is known to improve cognitive functions of patients (Oberholzer & Müri, 2019; Cernich et al., 2010). Neuroplasticity, that is, the adaptive changes of structural (Schmidt et al., 2020) and functional (Olafson et al., 2021) neural circuitry in terms of molecular, synaptic, and cellular changes, is commonly cited as a potential explanation for the cognitive and functional recovery (Sophie Su et al., 2016). Although axonal sprouting and functional rewiring post TBI is reported (Castellanos et al., 2010; Nakagawa et al., 2013), the underlying mechanism of change in white matter structural connectivity over time at the network level is still unclear (Zatorre et al., 2012).

The significant decline in network similarity of patients relative to healthy controls over time (Figure 4), may be indicative that the connectivity alterations happening in the network are mainly degeneration in connectivity rather than a recovery. This is in line with consistent neurodegeneration and neuronal loss that is widely reported in the TBI literature, which starts with injury and continues decades postinjury (Farbota et al., 2012; Graham & Sharp, 2019; Johnson et al., 2013). An alternative explanation for connectivity alterations in favor of structural recovery could be that the network topology of patients reorganizes to converge to a new normal unlike that of healthy controls to regain the network integrity. The decline of longitudinal change in the similarity of patients among themselves being steeper (Sections 3.3 and Supplementary Information S1) than that of their similarity relative to healthy controls (Figure 4), however, contradicts this alternative, further supporting the point that the alterations in the white matter network are not a recovery but a degeneration.

In contrast to the decline in their NNS, the cognitive recovery of patients over time (Figure S2) highlights an interesting disparity. When considered together with the lack of a significant association between the rate of change in NNS and cognition (Figure 5), it can be inferred that the structural changes in the network topology do not directly translate into cognitive recovery. Considering that TBI is a complex disease with multiple, potentially opposing, mechanisms at work simultaneously (Veenith et al., 2009), there might be several reasons for this apparently paradoxical disparity between structural connectivity degeneration and cognitive recovery (Farbota et al., 2012). One possible explanation is that neuroplasticity happens at the gray matter in terms of axonal sprouting more than white matter plasticity such as myelination. Supporting this perspective, axonal rewiring and sprouting in cortical gray matter are reported to happen in mice post TBI (Niogi et al., 2008). Several studies on functional MRI, which investigate connectivity of gray matter regions, reported network reorganization after TBI which correlates with cognitive recovery, providing further evidence to that option (Castellanos et al., 2010, Castellanos et al., 2011, Han et al., 2020). Complementing this perspective of synaptic plasticity, another mechanism at play could be that structural connectivity is disrupted at the time of injury, leading to cognitive deficit, due to axonal damage. Although those injured axons do not get repaired and are practically nonfunctional, some are captured by MRI as healthy fiber tracts connecting brain regions due to the coarse resolution of imaging data. This makes the network topology of a patient look similar to that of a healthy control. As the debris of the damaged axons gets removed from the network, on the other hand, network similarity of patients declines. Since injured axons do not function following the injury, their removal from the network does not have any effect on the cognitive scores of patients as it does not introduce any further disconnectivity into the network. Another potential explanation could be that the decline in network similarity of patients is due to a mixture of degeneration of neurons and strengthening of compensatory pathways, where the latter leads to cognitive recovery (Bengtsson et al., 2005).

We note that the positive correlation between the NNS and EF and PS do not contradict the earlier observations of network similarity declining over time while cognitive scores improve. Since higher NNS values indicate a lesser injury, better cognitive performance would be expected from such individuals as disconnectivity between regions will be lesser. Thus, the negative correlation between injury severity as quantified by PTA and both the NNS and cognitive scores support a positive correlation between the NNS and cognition.

LIMITATIONS, FUTURE DIRECTIONS, AND CONCLUSIONS

Although this study investigates a unique longitudinal TBI dataset with dMRI data and cognitive assessment acquired at three timepoints and uses an advanced graph theoretical technique, certain limitations should be acknowledged. First, diffusion MRI is known to have inaccuracies in determining connectivity between regions, such as its limitations in characterizing white matter in complex regions where fibers intersect (Jones, 2010). In the case of TBI, axonal injury causing the degeneration of one of the crossing fibers, for example, can result in increased FA over the other fiber, which in turn results in increased connectivity between two regions (Hayes et al., 2016). Diffuse injury can potentially introduce further noise to the structural connectivity, that in turn affects the resulting network, leading to a structure that is dissimilar to healthy controls as assessed by NNS score. Since such shortcomings are inherent to dMRI based analysis, and the proposed approach does not make an assessment of the integrity of the brain graph, the results presented here should be considered accordingly. Second, as typical of TBI studies, statistical significance of our results is limited by the sample size of TBI cohort (Wang et al., 2018). Also, our study lacks mild TBI patients, and it should be noted that the results may not translate to a lower injury severity. Similarly, as patients with large lesions were excluded from our analysis and the proposed measure was devised for diffuse injury, efficacy of the measure on TBI with large lesions requires further evaluation. In order to evaluate the trajectory of structural change in the acute as well as chronic phase of the disease across the injury spectrum, re‐evaluation of results presented here on a larger dataset (such as TRACK‐TBI, https://tracktbi.ucsf.edu/; Yue et al., 2013) is left as a future work. We also see a potential in evaluating the changes in functional network topology and its association with cognitive recovery on TBI, as such a recovery was recently reported on stroke patients using graph matching (Olafson et al., 2021). Finally, evaluation of the measure on data obtained by alternative preprocessing techniques as well as investigation of the change at functional systems and nodal levels is to be carried out on a future study with a larger cohort.

In conclusion, our results demonstrate that the structural brain networks of patients with moderate‐to‐severe TBI differ from those of healthy controls by 3 months and become increasingly different up to 1 year postinjury. It also demonstrates the efficacy of our network normality score (NNS) as a principled measure for evaluating severity of diffuse injury, which can have potential uses in creating diagnostic and prognostic biomarkers of the disease when evaluated on larger datasets. Moving forward, we will expand our method to investigate changes in network topology of functional connectivity in TBI patients, in order to explore mechanisms of cognitive recovery with an overall network analysis perspective.

CONFLICTS OF INTEREST

The authors certify that they have NO affiliations with or involvement in any organization or entity with any financial interest (such as honoraria; educational grants; participation in speakers' bureaus; membership, employment, consultancies, stock ownership, or other equity interest; and expert testimony or patent‐licensing arrangements), or nonfinancial interest (such as personal or professional relationships, affiliations, knowledge or beliefs) in the subject matter or materials discussed in this manuscript.

AUTHORS CONTRIBUTIONS

The authors confirm contribution to the paper as follows: study conception and design: Yusuf Osmanlıoğlu, Drew Parker, Ragini Verma; methodology and software: Yusuf Osmanlıoğlu; formal analysis: Yusuf Osmanlıoğlu; data collection: Junghoon J. Kim, John Whyte; data curation: Drew Parker, Jacob A. Alappatt; analysis and interpretation of results: Yusuf Osmanlıoğlu, Drew Parker, James J. Gugger, Ramon R. Diaz‐Arrastia, John Whyte, Junghoon J. Kim, Ragini Verma; writing original draft: Yusuf Osmanlıoğlu, Drew Parker, Ragini Verma; writing review and editing: Yusuf Osmanlıoğlu, Drew Parker, Jacob A. Alappatt, James J. Gugger, Ramon R. Diaz‐Arrastia, John Whyte, Junghoon J. Kim, and Ragini Verma. All authors reviewed the results and approved the final version of the manuscript.

APPENDIX S1 Supporting information

Click here for additional data file.