Functional outcome is tied to dynamic brain states after mild to moderate traumatic brain injury.

The current study set out to investigate the dynamic functional connectome in relation to long-term recovery after mild to moderate traumatic brain injury (TBI). Longitudinal resting-state functional MRI data were collected (at 1 and 3 months postinjury) from a prospectively enrolled cohort consisting of 68 patients with TBI (92% mild TBI) and 20 healthy subjects. Patients underwent a neuropsychological assessment at 3 months postinjury. Outcome was measured using the Glasgow Outcome Scale Extended (GOS-E) at 6 months postinjury. The 57 patients who completed the GOS-E were classified as recovered completely (GOS-E = 8; n = 37) or incompletely (GOS-E < 8; n = 20). Neuropsychological test scores were similar for all groups. Patients with incomplete recovery spent less time in a segregated brain state compared to recovered patients during the second visit. Also, these patients moved less frequently from one meta-state to another as compared to healthy controls and recovered patients. Furthermore, incomplete recovery was associated with disruptions in cyclic state transition patterns, called attractors, during both visits. This study demonstrates that poor long-term functional recovery is associated with alterations in dynamics between brain networks, which becomes more marked as a function of time. These results could be related to psychological processes rather than injury-effects, which is an interesting area for further work. Another natural progression of the current study is to examine whether these dynamic measures can be used to monitor treatment effects.

INTRODUCTION

Traumatic brain injury (TBI), and especially mild TBI (mTBI), affects millions of people around the world each year (Gardner & Yaffe, 2015). Approximately 30–40% of patients with mTBI do not fully recover at 6 months postinjury, although rates vary, depending among other things on the selected cutoff value for good recovery (McMahon et al., 2014; van der Naalt et al., 2017). Various studies have been published on prediction of functional outcome after mild and moderate TBI (Einarsen et al., 2018; Jacobs et al., 2010; Lingsma et al., 2015; van der Naalt et al., 2017; van der Naalt, van Zomeren, Sluiter, & Minderhoud, 1999). Many factors (demographic, injury‐related, radiological, and psychological) have been identified as important predictors of outcome after TBI. Interestingly, measurement of psychological factors early after injury has been shown to add significantly to the prediction of outcome (van der Naalt et al., 2017). Despite these studies, there is still a paucity of research on the actual neural mechanisms underlying recovery after TBI. Many functional MRI (fMRI) studies have reported relationships between persistent symptoms and neural networks (Mayer, Mannell, Ling, Gasparovic, & Yeo, 2011; Palacios et al., 2017; van der Horn et al., 2017; Wilson, Pettigrew, & Teasdale, 1998). However, the neural basis of functional outcome, such as that measured using the Glasgow Outcome Scale Extended (GOS‐E), may be different. To the best of our knowledge, this topic has not yet been thoroughly investigated.

Studying separate brain areas and networks is still at the center of gravity in resting‐state fMRI research, while results vary extensively across studies (Mayer, Bellgowan, & Hanlon, 2015). In fact, various research has shown that traumatic brain injury, and mTBI in particular, may be a disorder of global rather than local brain connectivity, with affected local brain areas being quite heterogeneous (Eierud et al., 2014; Kaushal et al., 2019; Pandit et al., 2013; van der Horn et al., 2016). Therefore, by focusing on separate brain areas, one may not fully apprehend the neural mechanisms related to recovery. This approach also obligates researchers to apply multiple comparison corrections, which may give rise to type II errors, especially in such a heterogeneous condition as mTBI (Mayer et al., 2015). Additionally, up until now, most studies have focused on static functional connectivity (connectivity between regions) or static functional network connectivity (FNC, connectivity between temporally coherent networks) measures, while there is preliminary evidence showing that time‐related patterns of FNC (i.e., dynamic FNC) may provide more detailed information about injury effects and recovery after TBI (Hou et al., 2019; van der Horn et al., 2016; Vergara, Mayer, Kiehl, & Calhoun, 2018). The brain's neural landscape is a dynamic system, which changes its configuration in response to ongoing cognitive and emotional demands (Du, Fu, & Calhoun, 2018). Therefore, it may be right to assume that neurological and psychological recovery after TBI is captured more accurately using dynamic instead of static connectivity measures. Furthermore, the neural underpinnings of recovery may change as a function of time, and may not be fully revealed using only a single brain scan. To the best of our knowledge, no previous study has applied dynamic functional connectivity analyses to longitudinal fMRI data of patients with mild to moderate TBI.

In the current longitudinal resting‐state fMRI study, patients with mild to moderate TBI were scanned on average 1 and 3 months postinjury. The main objective of the study was to investigate whether specific alterations in static and dynamic functional network connectivity explain poor long‐term functional outcome at 6 months post‐TBI, as measured with the GOS‐E. In addition, it was questioned whether neural patterns associated with recovery would become more apparent with the passing of time, and whether longitudinal patterns of connectivity are informative of outcome.

METHODS

Participants

Sixty‐eight patients with mild (n = 63; 92.6%), or moderate (n = 5; 7.4%) TBI were prospectively enrolled between March 2013 and February 2015 at a level 1 trauma‐center in the Netherlands (University Medical Center Groningen). Mild TBI was defined by loss of consciousness of maximally 30 min, a Glasgow Coma Scale (GCS) score of 13–15 afterward, and a period of posttraumatic amnesia (PTA) no longer than 24 hr according to the American Congress of Rehabilitation Medicine criteria (Kayd et al., 1993). Moderate TBI was defined by loss of consciousness of 30 min or more, a GCS score of 9–12, or posttraumatic amnesia lasting more than 24 hr (Einarsen et al., 2018; Godoy, Rubiano, Rabinstein, Bullock, & Sahuquillo, 2016; Malec et al., 2007). The majority of patients sustained their injuries due to traffic accidents and/or falls. Fourteen patients (nine mTBI and all moderate TBI) had lesions on computed tomography (CT) scans made at time of admission to the emergency department, with one mTBI patient requiring neurosurgery to evacuate an epidural hematoma. Supporting Information S1 contains a table with lesion characteristics for this subgroup of patients.

Exclusion criteria for this fMRI‐study were neurological or psychiatric comorbidity, hospital admission for prior TBI, drug or alcohol abuse, insufficient comprehension of the Dutch language, intellectual disability, and contraindications for MRI (any implanted ferromagnetic devices and objects, pregnancy, and claustrophobia). In addition to the patient‐group, 20 healthy controls (HC) without a history of TBI were enrolled. The HC‐group was matched to the TBI‐group with respect to age, biological sex, education, and handedness.

The study was approved by the local Medical Ethics committee of the University Medical Center Groningen. All participants provided written informed consent. Study procedures were conducted according to the declaration of Helsinki. Part of the data was used in previous analyses (van der Horn, Liemburg, et al., 2016).

Clinical measures

The main purpose of this study was to examine FNC related to outcome after TBI. At 6 months postinjury, functional outcome was determined using the GOS‐E, which is an eight‐point outcome measure (8 = upper good recovery; 7 = lower good recovery; 6 = upper moderate disability; 5 = lower moderate disability; 4 = upper severe disability; 3 = lower severe disability; 2 = vegetative state; and 1 = death) (Wilson et al., 1998). The GOS‐E interview contains questions about consciousness, independence at and outside home, work, social and leisure activities, and family and friendships. For analyses, a dichotomy was made between complete (CR; GOS‐E = 8), and incomplete recovery (ICR; GOS‐E ≤ 7). In addition to the GOS‐E, persistent posttraumatic symptoms were measured using the Head Injury Symptom Checklist at 6 months postinjury (de Koning et al., 2016; van der Naalt et al., 1999).

At 3 months postinjury (at the time of the follow‐up scan), a neuropsychological assessment was performed comprising: the Trailmaking test A (processing speed) and B (executive functioning; Reitan & Wolfson, 1985), Stroop I test (verbal speed; Jammes & Hammes, 1971), Digit Span test backwards (working memory; Wechsler, 2001), Dutch version of the Rey Auditory Verbal Learning test (RAVLT; immediate and delayed verbal memory; Rey, 1964), and the Controlled Oral Word Association test (COWAT; verbal fluency; Benton, de Hamsher, & Sivan, 1983). These tests encompass cognitive domains that match symptoms that are often reported after mTBI. Raw scores were used for analyses, corrected for age and education. To account for underachievement, the Amsterdam Short Term Memory test (ASTM; cutoff <85; Schmand, de Sterke, & Lindeboom, 1999), and a brief version (neurological impairment and amnestic disorders subscales) of the Structured Inventory of Malingered Symptomatology (SIMS; cutoff >5; Smith & Burger, 1997) were also administered.

MRI‐acquisition protocol

Patients were scanned at approximately 4 weeks (M = 37.2 days; range 22–69), and at 3 months (M = 96.5 days; range 61–207) postinjury. Of these 68, 63 patients (92.6%) returned for follow‐up scanning. The average time between first and second scan was 59 days. Healthy controls were scanned once. All MRI‐scans were collected using a 3 Tesla Philips Intera Achieva MRI‐scanner (Philips Medical Systems, Best, The Netherlands) equipped with a 32 channel SENSE head‐coil. An anatomical (T1‐weighted) image was made for referencing purposes (repetition time [TR] 9 ms, echo time [TE] 3.5 ms, flip angle 8°, field of view [FOV] 256 × 232 × 170 mm, reconstructed voxel size 1 × 1 × 1 mm). For resting‐state fMRI, 300 T2*‐weighted volumes were acquired with slices aligned to the anterior commissure (AC)–posterior commissure (PC) plane in descending order (TR 2000 ms, TE 20 ms, FOV 224 × 224 × 136.5 mm, reconstructed voxel size 3.5 × 3.5 × 3.5 mm). Subjects were instructed to keep their eyes closed and to stay awake during scanning.

Functional MRI preprocessing

Preprocessing was performed using statistical Parametric Mapping (SPM12; Wellcome Department, University College London, England) implemented in MATLAB (version 2017b; MathWorks, Natick, MA). Volumes were reoriented to the AC, slice time corrected, realigned, coregistered to the T1‐anatomical image, normalized to the Montreal Neurological Institute (MNI) space using an echo planar imaging (EPI)‐template (3 × 3 × 3 mm isotropic voxels) (Calhoun et al., 2017), and smoothed using a 8 mm full width at half maximum (FWHM) kernel. The first five volumes were discarded to account for lack of T1‐equilibrium. Before smoothing, voxel time courses were orthogonalized with respect to the six realignment parameters, their derivatives, and to linear, quadratic, and cubic trends (Vergara, Mayer, Damaraju, & Calhoun, 2017; Vergara, Mayer, Damaraju, Hutchison, & Calhoun, 2017).

Two patients and one HC showed extreme head movement (defined as >3 SD on at least 2/6 of the mean framewise displacement [FD] measures; Mayer, Franco, Ling, & Cañive, 2007; Mayer, Ling, et al., 2015). To reduce the possibility of spurious findings, the healthy subject with extreme head movement was excluded from all analyses. As the two patients showed extreme movement only during follow‐up scanning, they were excluded from the follow‐up analyses only. Total mean FD did not significantly differ between recovery and HC subgroups (one‐way ANOVA for patients with ICR, CR and HC, p = .36). Supporting Information S2 contains the mean FD values for all subjects.

Independent component analysis

The Group ICA of fMRI Toolbox (GIFT) version 4.0b (implemented in MATLAB) was used to perform group‐information guided independent component analysis (GIG‐ICA; Calhoun, Adali, Pearlson, & Pekar, 2001; Du et al., 2016; Du & Fan, 2013). Seventy components from an independent study on mTBI were used as templates for gig‐ICA (Vergara et al., 2018). Back‐reconstruction using gig‐ICA has been shown to result in increased sensitivity to group‐differences as compared to spatiotemporal regression (Salman et al., 2019). The end result of GIG‐ICA was a set of 48 resting‐state networks (RSNs) that were organized into nine network domains: subcortical (SBC), cerebellum (CER), auditory (AUD), sensorimotor (SEN), visual (VIS), salience network (SAL), default mode network (DMN), executive control network (ECN), and language (LAN). Figure 1 displays the functional network analyses pipeline.

Figure 1

Functional network analysis pipeline

Static FNC analyses

Static FNC was computed using GIFT. First, RSN's time courses were despiked using the Analysis of Functional NeuroImages (AFNI) 3D Despike software (AFNI;, 1995), and filtered using a fifth‐order Butterworth bandpass filter (0.01–0.15 Hz). Subsequently, pairwise correlations were computed, resulting in 1128 correlation values (48*[48-1]/2), which were transformed using a Fisher's Z‐transform prior to statistical analyses.

Dynamic FNC analyses

The dynamic FNC toolbox v1.0a in GIFT was used to perform dynamic FNC (dFNC) analyses. Time courses were despiked and filtered using a fifth‐order Butterworth bandpass filter (0.01–0.15 Hz), and pairwise correlations were calculated using a sliding window approach for rectangular windows (22 TR = 44 s) convolved with a Gaussian (σ = 3TR) in steps of 1 TR, resulting in 273 windowed correlation matrices (Damaraju et al., 2014; Miller et al., 2016; Vergara, Abrol, & Calhoun, 2019). This procedure was followed by computation of average sliding window correlations (ASWC) using a window of 25 TR (50 s) in MATLAB, resulting in 248 final correlation matrices with each containing 1,128 component pairs (Vergara, Abrol, & Calhoun, 2019). Subsequently, k‐means clustering was performed (maximum iterations of 4,000, correlation distance method, 33 replicates) using both the ASWC as well as their first derivatives of all subjects (derivatives were normalized for every subject to match the variance of the correlations). Including the first derivatives has been shown to result in higher sensitivity for capturing group‐differences in dynamic FNC (Espinoza et al., 2019). To estimate the optimal number of clusters (k), the following cluster validity indexes were computed for a range of k (1–10) using GIFT: Elbow method, Bayesian information criterion (BIC), Dunn's index, Gap statistic, Akaike information criterion (AIC), and Silhouette method (Supporting Information S1). The optimal number of clusters was determined to be five, which was the mean of the optimal numbers for aforementioned cluster validity indices rounded to the nearest integer greater than or equal to that value. From now on, clusters will be referred to as dFNC states.

State clustering analyses

The k‐means clustering algorithm provides a state‐membership index for every correlation window for each subject (i.e., the cluster that most closely matches the dFNC matrix of that window; this method is also referred to as hard clustering). These indices were used to calculate mean dwell time (MDT; i.e., average time spent in a particular state before changing to another state), fraction of time spent per state (FT; i.e., time spent in one state divided by the total number of windows), and number of transitions (NT; i.e., how often a subject changes states) for every subject (using GIFT; Calhoun, Miller, Pearlson, & Adalı, 2014). Note that it can occur that a subject does not visit at certain state, in this case mean dwell time, and fraction of time spent were classified as zero for that particular state.

Meta‐state analyses

These additional analyses provide a more detailed view of dFNC states (Miller et al., 2016). Meta‐state analysis was run in GIFT, again by using a k‐means clustering algorithm on ASWC together with the first derivatives (maximum iterations of 4,000, correlation distance method, 33 replicates). Subsequently, weighted k‐means time courses were computed using the k‐means correlation patterns. Instead of assigning one particular state‐membership index to a particular dFNC window, the distance of a particular dFNC matrix to every of the five states is computed (i.e., soft clustering). The resulting five‐element vector is called a meta‐state (note that because for every window a vector reflects the distances of that window to every state, no zero values are present in these meta‐state vectors). Measures similar to those that were mentioned earlier were calculated using these meta‐states (i.e., state span, number of state‐changes). In addition, the total distance traveled through the state‐space (sum of distances between successive meta‐states), the number of states (i.e., unique meta‐states), and the number of hub‐states (a unique meta‐state that is visited ≥4 times throughout the windows) were computed.

State‐transition analyses

The ASWC and their first‐order derivatives were used for state‐transition analysis, which estimates the probability of transitioning from one state to another. Mathematically, state transitions can be written as P{, with representing the k‐means state membership index before transition, and representing the one after transition (Vergara et al., 2019). These probabilities were estimated in a group‐wise manner using the membership‐indices per window. The resulting probability matrices were analyzed for the presence and disruption of cyclic patterns. Under normal circumstances, a transition is followed by another transition that returns to the initial membership. For example, a transition from membership one to membership two (1 → 2) is followed by a two to one transition (2 → 1). We call this set of transitions an attractor. Attractors consist of connectivity patterns that oscillate, or in other words, dynamic FNC states that will transit to each other with the highest probability. States within an attractor have identical ASWC, but different derivatives, reflecting different patterns of increasing and decreasing (oscillating) functional connectivity over time. These time‐varying patterns of connectivity strength have been shown to orbit around a centroid, which is the center of an attractor. Supporting Information S1 contains a detailed summary of the state transition analyses.

In addition to the patient subgroups, the group of HC was examined to identify healthy transitions between states.

Statistical analyses

For clinical measures, analyses were performed using the Statistical Product and Service Solutions (SPSS, version 20, IBM Corp., Armonk, NY). Student's t tests and analyses of variance (ANOVA) and covariance (ANCOVA) were used for group comparisons of continuous data. For neuropsychological data, age and education were included as covariates. For nominal and ordinal data, Chi‐square tests were used. Alpha was set at 0.05 (two‐tailed).

Since there was a significant difference in age between patients with CR and ICR, variance related to age was regressed out of the data prior to performing statistical analyses of network measures. Regarding transition analyses, estimated probabilities rely on data from a whole group and there is no subject‐wise information. For these reasons, age could not be regressed out of these data. Therefore, the two groups of patients with ICR (n = 20 and 20 for Visits 1 and 2, respectively) and CR (n = 29 and 27) were compiled in such a way that they were age‐matched at both visits, with sample sizes as large as possible. We also ensured groups remained gender‐ and education‐matched.

Permutation tests were performed using MATLAB (10,000 permutations) to compare static connectivity values (correlation values of each component‐pair), state clustering and meta‐state measures. Permutation tests for one‐way ANOVA were used to test for the presence of any group‐differences between HC, patients with CR, and patients with ICR (α = 0.05, two‐tailed; Anderson, 2001). In case of a significant group‐effect, post hoc pairwise permutation tests were performed. To control type I errors, FDR‐corrections were applied. For static FNC, corrections were applied for the number of component pairs (i.e., 1,128 tests). For the dynamic measures MDT and FT, corrections were applied for the number of states (i.e., five tests). To assess longitudinal effects, paired permutation tests were performed for the CR and ICR groups, and delta‐values expressing longitudinal change were compared between groups using two‐sample permutation test.

In addition to testing for statistical significance, eta‐squared (η2; for ANOVA), and Hedges’ g (for two‐sample comparisons) effect sizes were computed using the Measures of Effect Size Toolbox in MATLAB (Hentschke & Stüttgen, 2011).

To examine the influence of measures of injury severity (GCS and/or PTA) on our findings, group analyses were conducted with and without inclusion of patients with moderate TBI. Supplementary analyses were conducted to assess the influence of CT‐lesions by comparing network measures between patients with and without CT‐lesions.

RESULTS

Participants’ characteristics

A total of 57 patients (84%) completed the GOS‐E (37 with CR and 20 with ICR), and 54 patients (79%) of these patients were included in the longitudinal analyses (34 with CR and 20 with ICR).

Table 1 lists demographic and neuropsychological test scores for patients with ICR, CR, and HC. Age was significantly different between patients with CR and ICR (M = 35.1 vs. 44.3; p = .024; t = 2.32). Gender and education were similar for both groups. Mean GOS‐E‐score for patients with ICR was 6.2 (range 5–7). Number of symptoms for patients with CR (M = 3.2 [0-14]) and ICR (M = 8.4 [3-18]) differed significantly at 6 months postinjury (p < .0001; t = 4.49). It has to be noticed that despite full recovery, patients with CR may still report symptoms (e.g., 15 patients still reported 3 or more symptoms).

Table 1

Participant characteristics

	ICR (n = 20)	CR (n = 37)	HC (n = 19)	p‐value; effect size
Age, years	44.3 (22–61)	35.1 (19–64)	35.9 (18–61)	.06; F = 2.89
Sex, % male	60	76	68	.47; χ2 = 1.53
Education level, median (range)a	5.5 (4–7)	6 (2–7)	6 (5–7)	.45; χ2 = 7.84
GCS‐score, median (range)	14.5 (9–15)	15 (13–15)	N/A	.24; χ2 = 4.23
Days between injury and first visit	36.3 (25–69)	36 (22–67)	N/A	.93; t = 0.09
Days between injury and second visit	99.7 (61–207)	93.5 (76–127)	N/A	.31; t = 1.03
Percentage (n) with moderate TBI	15 (3)	3 (1)	N/A	.08; χ2 = 3.01
Percentage (n) with CT‐lesions	25 (5)	16 (6)	N/A	.42; χ2 = 0.64
Neuropsychological tests b	n = 16	n = 31	n = 18
TMT‐A	28.1 (14–49)	28.7 (14–74)	26.9 (12–43)	.61; F = 0.50
TMT‐B	58.5 (34–99)	57.1 (27–115)	60.3 (28–114)	.70; F = 0.37
Stroop‐I	45.5 (30–60)	45.5 (31–72)	45.2 (26–71)	.71; F = 0.34
Digit‐span backwards	5.6 (3–9)	5.2 (2–8)	5.3 (3–8)	.56; F = 0.58
RAVLT
Immediate	44.8 (17–56)	47.6 (31–69)	47.4 (32–71)	.99; F = 0.01
Delayed	9.3 (4–15)	10.0 (2–15)	9.9 (3–15)	.98; F = 0.02
COWAT	36.4 (20–61)	38 (11–59)	41.8 (21–62)	.53; F = 0.65

Note: Values are expressed as mean (range), unless stated otherwise. All statistical analyses for (raw) neuropsychological test scores included age and education level as covariates.

Abbreviations: COWAT, Controlled Oral Word Association Test; CR, complete recovery; HC, healthy controls; ICR, incomplete recovery; MRI, Magnetic Resonance Imaging; N/A, not applicable; GCS, Glasgow Coma Score; RAVLT, Rey Auditory Verbal Learning test; TMT, Trailmaking test.

Education level was based on a Dutch classification system, according to Verhage (1964), ranging from 1 to 7 (highest).

Subjects who underachieved based on the SIMS and ASTM tests were excluded.

Our previous work has shown that posttraumatic symptoms are not related to microhemorrhagic lesions on T2‐gradient echo and susceptibility‐weighted imaging scans (van der Horn, de Haan, Spikman, de Groot, & van der Naalt, 2018). In the current study, outcome was not related to the number of microhemorrhagic lesions in patients who did not have lesions on CT (i.e., uncomplicated mTBI, ICR vs. CR: p = .75, t = 0.32).

The mean static FNC matrices for subgroups (HC, ICR, and CR) are depicted in Figure 2. In all matrices, it is noticeable that networks for sensorimotor, visual and auditory processing (CER, AUD, SEN, VIS) are strongly intra‐ and inter‐connected, as are networks for higher cognitive processing (SAL, DMN, ECN), while between these two sets of domains there are mostly anti‐correlations.

Figure 2

Mean static FNC matrices for patients with ICR, CR, and HC

A permutation test for one‐way ANOVA revealed several group‐effects, however, differences for none of the component‐pairs survived FDR‐corrections during post hoc group comparisons. There were no significant longitudinal effects within the ICR nor CR group. The matrices containing the log‐transformed uncorrected p values with sign of the associated group differences show a trend toward weaker correlations within and between networks for sensorimotor, visual and auditory processing (CER, AUD, SEN, and VIS) in patients with ICR (particularly during the second visit) compared to CR and HC, while these networks were more strongly interconnected with networks for higher cognitive processing (SAL, DMN, ECN; Figure 3).

Figure 3

Group differences in static FNC displayed as sign(mean difference in correlation) * −log10(uncorrected p‐values). Blue vertices in the matrices show trends toward decreased static FNC, red vertices reflect trends toward increased FNC

Supporting Information S1 contains additional analyses comparing patients with and without lesions on CT.

Dynamic FNC states

Cluster centroids for ASWC as well as the first derivatives for every of the five dFNC states are depicted in Figure 4. It is evident that state one and two show similar ASWC patterns, with mostly anti‐correlations between the CER/AUD/SEN/VIS domains on the one hand, and the SAL/DMN/ECN on the other, resembling the segregated static FNC patterns. Interestingly, one DMN component (posterior cingulate) and two ECN components (superior parietal) in states one and two show positive correlations with the CER/AUD/SEN/VIS domains. States one and two exhibit opposite derivative patterns, reflecting either increasing or decreasing correlations over time. The ASWC matrices for States 3, 4, and 5 look similar, with higher between‐domain correlations as compared to States 1 and 2; their derivatives show different patterns. Supporting Information S3 contains ASWC‐ and derivative‐matrices for every state for patients with ICR and CR (both visits), and HC.

Figure 4

Centroids of the five dFNC‐states for the total group of participants (TBI and HC). Per state the number of windows as well as the percentage of total windows assigned to that particular state across all participants is depicted. The top right graph shows the number of subjects visiting each state for the TBI (first and second visit) and HC groups

Permutation ANOVA revealed a group‐effect regarding MDT (p = .041, F = 3.31, η2 = .09) and FT (p = .045, F = 3.19, η2 = .08) for State 2 during the second visit. Post hoc tests showed that FT for State 2 was significantly lower in patients with ICR compared to patients with CR (p = .007, g = −0.78; Figure 5). This effect was also observed at trend‐level for MDT (p = .015, g = −0.7, p FDR = .075). No statistically significant differences in MDT or FT were found between either of the patient subgroups and HC. There was a trend toward a significant difference in delta (Visit 1 − Visit 2) regarding MDT for State 1 between patients with ICR and CR (p = .013; g = 0.71; p FDR = .063); however, within these groups there were no longitudinal changes (p = .08 and .11 for ICR and CR, respectively). There were no significant effects for NT. No significant group‐differences or trends toward significance were found for the first visit (with small effect sizes [η2 = 0.01] for MDT and FT in State 2).

Figure 5

Fraction of time spent per state for patients with ICR and CR during the second visit, and HC. Mean, 95% confidence interval, and one standard deviation are shown. Asterisk (*) indicates a significant group difference after FDR correction. Note that because age was regressed out of the data, values can be negative

Supporting Information S1 contains additional analyses comparing patients with and without lesions on CT.

A group‐effect was found for total distance (p = .047, F = 3.08, η2 = .08) and number of meta‐state changes (p = .045, F = 3.13, η2 = .08) during the second visit (Figure 6). Post hoc comparisons showed a significantly lower total distance (p = .035, g = −0.69) and number of meta‐state changes (p = .034, g = −0.69) in patients with ICR as compared to HC. This was also observed, although at trend‐level, in patients with ICR compared to CR (p = .06 and g = −0.54; p = .086 and g = −0.49, respectively). There were no significant longitudinal effects. Regarding number of meta‐states, state span and number of hub‐states, there were no significant group effects. No significant group‐differences or trends toward significance were found for the first visit.

Figure 6

Total distance traveled through meta‐state space for patients with ICR and CR during the second visit, and HC. Mean, 95% confidence interval, and one standard deviation are shown. Note that because age was regressed out of the data, values can be negative. Asterisk (*) indicates a significant group difference after FDR correction. Pound sign (#) indicates a trend (p < .1)

Supporting Information S1 contains additional analyses comparing patients with and without lesions on CT.

Analysis of the group of HC revealed the presence of two attractors: (a) consisting of dFNC States 1 and 2, and (b) consisting of States 3, 4, and 5 (Figure 7). In patients with ICR, these attractors were disrupted during the first and second visit, with a more chaotic pattern of transition, and crosslinks between the two attractors. In patients with CR, the attractors themselves resembled the ones in HC, although there was a crosslink between the two during both visits.

Figure 7

State transition probabilities and schematic representations of attractors for patients with ICR and CR and HC. Additional rendered brain images show aggregated component T‐maps weighted based on the values in the attractor's ASWC‐ and derivative‐centroids (i.e., component's T‐map * [sum of all values for that component in the centroid, with subsequent variance normalization]). Red and blue colors in the attractor's center image reflect components with positive and negative correlations, respectively. In the surrounding state images, red and blue reflect time‐related (i.e., over windows) increasing and decreasing correlations, which are based on the derivatives. Thus, there are oscillating patterns of connectivity, orbiting an attractor's center. These patterns are disrupted in patients with ICR

Supporting Information S1 contains additional analyses comparing state transition patterns between patients with and without lesions on CT.

Analyses regarding injury severity measures

For static FNC, results remained consistent after exclusion of patients with moderate TBI. Also, state clustering results remained significant after exclusion of patients with moderate TBI, and patients with ICR now showed a significant longitudinal decrease in MDT for State 1 (p = .009, g = −0.43), with a significant difference in delta between ICR and CR (p = .003, g = 0.89). There were still no significant effects for NT. Permutation ANOVA for total distance and number of meta‐state changes did no longer show statistically significant group‐effects for the second visit after exclusion of patients with moderate TBI (both p = .12), although there were still medium effect sizes (η2 = .06 for both total distance and meta‐state changes), and the post hoc group‐differences remained similar.

State transition results also remained the same after exclusion of patients with moderate TBI, except for attractor B during the first visit for patients with ICR, which may seem more disrupted (Supporting Information S1).

DISCUSSION

To the best of our knowledge, this is the first longitudinal resting‐state fMRI‐study showing that dynamic functional network connectivity analyses provide information on functional outcome after TBI. The dynamic neural blueprint of incomplete recovery was characterized by an increased global connectivity pattern, which stood out in the early chronic phase postinjury. More specifically, patients with incomplete recovery spent less time in a dynamic brain state that was characterized by the presence of two main functional subdomains, one for sensory and motor function, and one for higher cognitive processing. Furthermore, they exhibited fewer shifts between dynamic brain states, and had more chaotic attractor dynamics.

In addition to static functional connectivity analyses, the current study examined specific time‐related patterns of functional connectivity, defined as dynamic brain states, in patients with mainly mild TBI. Previous studies have already shown the potential of dynamic connectivity analyses in the investigation of posttraumatic symptoms, and classification of patients with mild TBI (Hou et al., 2019; van der Horn, Liemburg, et al., 2016; Vergara et al., 2018). In the present study, k‐means clustering was used to partition the time‐varying functional connections into patterns of connectivity, defined as dynamic brain states. The emerging five states were further classified as attractors using a relatively novel technique focusing on transition probabilities (Vergara, van der Horn, et al., 2019). Findings from our healthy control group showed the presence of two main attractors, which were further studied in relation to outcome after TBI. Our study also included meta‐state analyses, which encompasses a soft clustering approach providing more detailed time‐varying patterns of functional connectivity and dynamic state behavior (Miller et al., 2016). To the best of our knowledge, this method has not been applied so far to study patients with TBI.

Cluster analyses showed that patients with incomplete recovery spend less time in the segregated dynamic brain state two. The corresponding attractor A was also found to be disrupted, with cross‐transitions from attractor A to B (Figure 7). Dynamic brain state one and two were characterized by high connectivity within the sensorimotor, auditory and visual domains on the one hand, and within the salience, default mode, and executive domains on the other, with anti‐correlations between these two. Interestingly, our static functional connectivity results also pointed toward a reduction in segregation of these anti‐correlated network domains in patients with incomplete recovery at trend level. Altogether, these findings suggest a shift toward more interconnected brain dynamics, which could reflect adaptive processes, or at least processes that are aimed at adaptation. It could also be the case that if these “adaptive” mechanisms fail or are employed for a too long period of time, this might have detrimental psychological effects (leading to a more “rigid” mental state), although we acknowledge the hypothetical nature of these explanations. Previous studies have reported increased connectivity (Hillary et al., 2014; Kaushal et al., 2019; Mayer et al., 2011), and less efficient organization of brain networks in patients with TBI relative to healthy subjects (Pandit et al., 2013). Our findings are also supported by a recent study showing that patients with persistent symptoms at 6 months after mTBI spend less time in efficient brain states (Hou et al., 2019). Although patients with incomplete recovery in the present study had significantly more symptoms than those with complete recovery, it is important to note that having persistent symptoms is not equal to incomplete recovery as measured with the Glasgow Outcome Scale Extended. In fact, our data show that patients can still report symptoms while functionally recovered. Interestingly, Hou and colleagues also reported alterations within auditory and visual networks associated with persistent symptoms. Functional outcome, persistent symptoms, preinjury mental health, psychological mechanisms (e.g., coping), and injury effects are intricately linked in mTBI, and therefore, we realize that measuring functional outcome does not fully capture the patterns of recovery after mTBI (Mayer, Quinn, & Master, 2017; van der Horn et al., 2019). However, the Glasgow Outcome Scale Extended is widely used internationally for studying outcome after mTBI, and for the development of prediction models (Nelson et al., 2019; van der Naalt et al., 2017). Other outcome measures such as (persistent) posttraumatic symptoms have proven less reliable, because patients with symptoms may still resume their daily activities, symptoms are only partly specific for mTBI, and definitions of persistent symptoms vary (de Koning et al., 2016; Quinn, Mayer, Master, & Fann, 2018; Voormolen et al., 2018).

Other research has reported trend‐significant increased static functional connectivity between visual, auditory and sensorimotor networks in TBI patients relative to healthy controls, with decreased dynamic functional network connectivity, defined as the variability of component‐pair correlation strength (Mayer, Ling, et al., 2015). Moreover, Palacios and collaborators also reported a negative correlation between functional connectivity within visual networks and symptom scores at 6 months postinjury (Palacios et al., 2017). Although our results should be interpreted with caution, we propose that the loss of network segregation shown in our study might be a reflection of enhanced awareness and higher‐order cognitive processing of sensory stimuli, possibly explaining symptoms such as heightened sensitivity to light and noise, headaches, and neck pain, which are symptoms that are often reported and impede daily functioning. In addition, increased connectivity of the default mode and executive network with the motor system may reflect an increased tendency of patients to act in response to cognitive processes, for example, with the intention of changing their current situation, thereby possibly overstressing themselves and their adaptive brain systems. Furthermore, the lower number of meta‐state changes and lower meta‐state‐traveling‐distance in patients with incomplete recovery might be a manifestation of a less flexible, more rigid mental architecture, possibly underlying an exaggerated focus on symptoms or on decreased level of daily functioning. In other words, our results may represent the neural basis of “getting stuck” in disadvantageous mental processes. Again, these interpretations are theoretical, and no supporting evidence from other studies can be presented at this point. Further studies need to be conducted to identify whether symptoms, mental processes and behavior are truly linked to specific changes in dynamic neural systems.

Distorted brain dynamics within the first weeks after injury is probably due to the brain injury itself and may persist in those patients with disadvantageous psychological characteristics, leading to poor long‐term outcome (van der Horn et al., 2019). In the current study, patients were first scanned at approximately 4–5 weeks postinjury, which impedes drawing inferences about connectivity changes earlier after injury. However, there are various studies that have shown changes in static and dynamic functional connectivity within the first weeks after injury (Hou et al., 2019; Kaushal et al., 2019; Mayer et al., 2011). Our results showed that patients with complete recovery had few disruptions of connectivity as compared to patients incomplete recovery. Previous research has shown that changes in cerebral blood flow and functional connectivity after mild TBI reach a peak within the first 2 weeks, and therefore, it is possible that more disruptions of connectivity would have been present in our study sample within the first weeks after injury (Kaushal et al., 2019; Meier et al., 2015).

Interestingly, a recent study has shown a negative relationship between time spent in a segregate brain state (comparable to states one and two in our study), and neuroticism scores in patients with major depressive disorder (Wu et al., 2019). This state was also characterized by strong connectivity of the default mode network with sensorimotor, visual and auditory networks, which we also observed for a posterior default mode network component in states one and two. According to the authors, the lack of this brain state in patients with high neuroticism scores could be related to a failure to process sensory stimuli and associated initiation of mood change, and to a reduced emotion regulation capacity. Our results might be viewed as manifestations of psychological rather than injury‐related processes. Patients with full recovery may still experience symptoms, which means that a considerable proportion of patients are able to adapt and cope with their symptoms, and may return to their normal daily activities, while others fail to do so. These patients are likely to have more beneficial psychological characteristics compared to patients with symptoms that show incomplete recovery.

Remarkably, regarding state clustering and meta‐state measures differences between patients with incomplete recovery and complete recovery, and healthy controls, seem to be most pronounced at 3 months postinjury. This is further illustrated by the observation that mean dwell time for state one seems to decrease over time in patients with incomplete recovery. It goes without saying that recovery manifests with time, and it is also likely that psychological effects (e.g., rumination about symptoms and reduced daily functioning) become more apparent with passing of time. We acknowledge that these explanations of delayed connectivity changes may be speculative. In fact, state transition analyses showed clear differences during both visits, which indicates that functional outcome has an early neural substrate. We have reason to believe that this state transition method is able to capture more subtle changes in network dynamics, possibly due to a lower susceptibility to noise, and because it provides a more integrated analysis than metrics on individual states by moving closer to a blend of dynamical systems and data driven approaches. Regarding the longitudinal pattern of functional connectivity changes after mTBI, studies so far are inconsistent. It has been reported that functional connectivity changes (probably related to the injury itself) may be most prominent within the first 1–2 weeks after injury, but delayed injury effects beyond this period have also been reported (Kaushal et al., 2019; Meier, Bellgowan, & Mayer, 2017). Our study specifically focused on outcome, and showed persistent changes in network behavior up to 3 months postinjury, that may be related to psychological factors, although the presence of injury‐related effects cannot be ruled out. The fact that crosslinks between attractor A and B were found patients with complete recovery, may be subtle manifestations of the latter.

The aim of the present study was to investigate outcome in patients across the entire spectrum of mild TBI, including patients with lesions on CT (i.e., complicated mTBI). Previous studies indicated the marginal added value of CT‐scans to prediction models of outcome after mTBI (Jacobs et al., 2010; van der Naalt et al., 2017). Supplementary analyses in the current study revealed alterations of state transition patterns in patients with lesions across both visits, while the group without CT‐lesions showed a restored attractor pattern at follow‐up. The group with CT‐lesions in our study was relatively small, but our results show that macroscopic lesions may impact dynamic brain states, thereby possibly influencing outcome. Despite the fact that more than 90% of the sample consisted of patients with mild TBI, the current study also included five patients who had sustained a moderate TBI. As injury severity is a composition of indicators such as loss of consciousness, Glasgow Coma Scale score and posttraumatic amnesia, unambiguous classification of injury severity is not always feasible. Grading of TBI is a continuum, with inherent difficulties of grading systems (Malec et al., 2007; Teasdale et al., 2014). For example, single indicators of injury severity, such as the Glasgow Coma Scale score can be underestimated (i.e., lower scores) due to confounding factors (e.g., intoxication, sedation effects, hearing impairment). Even though five patients were classified as moderate TBI, the majority of these patients were clinically close to the severe end of the mild TBI spectrum (Glasgow Coma Scale score in mild TBI range with posttraumatic anterograde amnesia lasting longer than 24 hr, Glasgow Coma Scale score in the moderate TBI range, but with speedy recovery of this score, posttraumatic amnesia lasting shorter than 24 hr, and/or no or short admission to the Intensive Care Unit). Furthermore, outcome scores for this group were not biased toward the lower end, and most importantly, results were mostly consistent after exclusion of patients with moderate TBI, and (meta‐)state measures for these individuals were within the range of the group they were part of (complete or incomplete recovery).

This study has limitations and caveats that need to be addressed. First, the group of healthy controls was relatively small, and they were only scanned once. This one scan was used for comparisons with patients’ data from both the first and second visits. Ideally, healthy controls would have also returned for follow‐up scanning. Second, the group of patients with incomplete recovery is relatively small compared to those who fully recovered, which may have affected our results due to reduced power. Although effect sizes were small, a larger sample size might have resulted in finding more (subtle) differences between complete and incomplete outcome during the first visit. Sample size was also one of the reasons for dichotomizing the outcome measure. In case of larger sample size, with a larger variability of outcome scores, ordinal regression may uncover valuable information. Third, the proportion of patients with moderate TBI is quite small compared to the ones with mild TBI. In future studies, we plan to include a larger sample of moderate TBI, to make better comparisons between these groups. Fourth, due to controlling for underachievement not all of the neuropsychological test data could be analyzed, which resulted in suboptimal power to discern differences, even though no differences were expected.

To conclude, the current study demonstrated alterations in dynamic brain states that are associated with recovery processes in patients with mild to moderate TBI. In our opinion, future studies should aim to also investigate these brain states at an earlier stage postinjury, to get a better sense of longitudinal changes related to outcome. This will lead to a better understanding of causal mechanisms of poor outcome, especially when it is coupled with measurement of psychological characteristics, such as personality traits (e.g., neuroticism). Furthermore, measurement of dynamic brain states may provide a future tool to identify patients who are at‐risk for poor outcome, which may facilitate the start of early interventions, and eventually, to monitor treatment effects.

Appendix S1. Supporting Information 1.

Click here for additional data file.

Appendix S2. Supporting Information 2.

Click here for additional data file.

Appendix S3. Supporting Information 3.

Click here for additional data file.