Resting Network Plasticity Following Brain Injury

1Traumatic brain injury (TBI) is a debilitating neurological disorder defined as an injury from an external source resulting in a period of altered consciousness and deficits in physical, cognitive, and/or psychosocial functioning. While examination of behavioral deficits associated with head trauma has a very long history, spanning over 60 years, only recently have the consequences of TBI received attention via functional imaging methods.

2“Activation” studies using functional MRI and positron emission tomography have been used to examine TBI-related deficits in episodic memory,, working memory –, and executive control. While the results of these studies have generated working hypotheses regarding how plasticity is expressed in disrupted neural systems, discrete regions of interest and localized “activation” results have not been interpreted in the context of an integrated neural network.

3There has been recent emphasis in studies using BOLD fMRI to approximate brain activity, to incorporate baseline or “resting” measurements of the BOLD signal. Systematic examination of baseline BOLD signal was first introduced by examining motor cortex and has recently received significant attention resulting in demonstration of a discrete system of networks that are non-task or “default mode”,. Emanating from these early findings has been a wellspring of studies examining resting brain states in the context of cognitive, sensory, and motor functioning. Most recently, these methods have been applied in cross-sectional work examining resting BOLD states in the clinical neurosciences.

4Resting state fMRI has thus provided unique information about the behavior of voxels (or “networks”) in the absence of direct stimulation. What has not been examined in this relatively new literature is if resting states are plastic, and, in particular, if they are changing after neural disruption. The current study aims to document change in resting neural networks during recovery from brain injury by examining macro-level functional connectivity in the BOLD fMRI signal via graph theory (described below). To date, there has been no work using serial MRI to examine changes in neural connectivity during recovery from neurological insult and such methods may provide additional insight into how neural plasticity is expressed in the injured brain. Such analyses may offer insight into how networks adapt to neurological disruption. For example, it remains unclear if the neural recruitment observed almost universally in cross-sectional “activation” studies of working memory deficit is due to formal brain reorganization or is indicative of neural inefficiency during periods of cognitive challenge,. What appears to be a critical element in making this determination is the nature of this neural recruitment over time and if the number of neural connections is altered during recovery. Activation studies in clinical samples can be methodologically challenging – and one potential method for examining how plasticity is expressed in the injured brain during recovery is to first document how resting networks are altered. Thus, it is an important aim to determine if, during recovery, networks become more elaborate, including the creation of additional connections, or if early recruitment of resources reflects use of available auxiliary networks, that later give way to diminished connectivity and greater neural efficiency. The whole-brain network analyses conducted here afford this opportunity.

5Graph theory provides a powerful framework for the mathematical treatment of complex systems such as neural networks in the brain –. The use of graph theory to date has led to important developments in how neural systems are understood, including small-world networks,,, scale-free network properties, “robustness for malfunction”,, and critical state in healthy subjects. Furthermore, early application of this type of network analysis to the study of neurological disorders, such as Alzheimer's disease, and schizophrenia,, has helped to characterize network abnormalities such as “neural disconnection”.

6In the current study, graph theory is used to examine changes in whole brain resting connectivity during recovery from moderate and severe TBI. To do so, we analyze “ resting” data in order to characterize changing network properties during a time period now well established as a critical window for recovery of cognitive functioning in TBI,.

7There are two over-arching goals in the current study. First, we aim to examine potential changes in resting connectivity during a period when behavioral recovery is known to occur after severe traumatic brain injury. To do so, resting data derived from working memory task data are analyzed at two time points early after injury. This approach affords the opportunity to observe the expression of neural plasticity during this critical window of recovery. Second, we use graph theory to examine changes in the resting BOLD response in order to examine whole-brain alterations in network connectivity over time. This approach represents a paradigm shift away from traditional cross-sectional work examining region-of-interest (ROI) changes associated with behavioral deficits. The inherent advantage to this approach is that whole-brain data analysis of “resting state data” effectively circumvents many of the methodological pitfalls associated with examining “activation” changes in clinical samples (e.g., task difficulty, performance inequalities) (see, ). Thus, the current study examines network change after injury via graph theory analyses in order to identify the how plasticity is expressed during recovery. Ultimately, these data may offer additional context for understanding the findings in the current cross-sectional literature using fMRI to examine deficits in working memory and rapid decision making.

8The averaged functional connectivity matrices for Time 1 and Time 2 were calculated by averaging absolute correlation matrices in the TBI sample (Fig. 1). The strength of functional connectivity diminished from Time 1 to Time 2, more closely approximating the results observed in the HC sample.

9Figure 2(a) shows group averaged probability distributions of absolute partial correlation coefficients for both groups. When comparing Time 1 and Time 2, lower correlation values ranging from r = 0.05 to 0.5 increased over time, but correlation values ranging from r = 0.6 to 0.95 decreased during recovery. The mean correlation coefficients were significantly different between time-points ( = 0.40±0.05 for Time 1 and 0.31±0.01 for Time 2), and values for Time 2 approximated those observed in the HC sample ( = 0.30±0.03).

10Though the sparsity of the network showed a non-significant difference between time points (Fig. 2(b)), the strength S significantly decreased at Time 2 (Fig. 2(c)), closely approximating the results observed in healthy adults. These results indicate that the frequency of strong functional connectivity decreased during recovery, but the overall number of connections in the network remained relatively stable from Time 1 to Time 2. The specific values for all indices are summarized in Table 1.

11Table 1 summarizes network properties of un-weighted functional brain networks. All indices characterizing network statistical properties did not show any significant differences between time points. Local efficiency and clustering coefficient C were higher at Time 1 relative to HC.

12The correlation matrix revealed varying strengths in functional connectivity dependent upon the brain region pairings (see correlation matrices, Fig. 1) and the distributions of correlation coefficients was influenced by time point of measurement (see Fig 2(a)). These findings suggest that the strength of functional connectivity may be altered during recovery from TBI. Because un-weighted network analysis does not consider region-specific connectivity strengths and relies solely upon a binary description of connectivity (e.g., 1 or 0), emphasis in the current study was placed upon results from weighted networks. Figure 3 and Figure 4 show network indices of weighted functional connectivity network. The global efficiency (Fig. 3(a)) and local efficiency (Fig. 3(b)) showed significantly higher values at Time 1 compared with Time 2 and Control. These findings are consistent with an overall reduction in “strength” in the connected network. Because of its inverse relationship to, path length L (Fig. 3(c)) significantly increased during recovery. Thus, as network strength (and reciprocally global and local efficiency) decreased during recovery and these index changes from Time 1 to Time 2 became more consistent with healthy adults.

13The results of indices which were normalized by producing randomized networks are shown in Fig 4. The components making up “small worldness”, changed from Time 1 to Time 2, approximating what was observed in healthy adults over time (see Fig. 4(a)). The small-worldness index was significantly lower at Time 1 compare with Control, indicating the loss of small-worldness following injury. Of note, the two estimates of “length” in connectivity, path length L and, were significantly different in opposite directions. Because is a measure of path length standardized via the randomized networks, further discussion focuses on this finding. The specific values of all indices for weighted network analysis are summarized in Table 2.

14One concern regarding network “weighting” is that findings may be at least partially attributable to statistical thresholding effects. When considering alternative statistical thresholds, the current results were observed to be quite robust; the results are not altered when comparing significant levels of p<0.05 and p<0.01 (see Tables S2 and S3).

15For both groups the cumulative degree distribution followed a non-Gaussian distribution with stretch-exponential-like heavy tails (Fig. 5). The distribution of Time 2 showed the similar decay as that of the HC group (Fig. 5(b)). By contrast, the distribution of Time 1 showed slower decay than that of Time 2 and the HC group, indicating an increase in the number of regions highly interconnected regions at Time 1 (Fig. 5(a)). Furthermore, the distribution of Time 1 was similar to the form of power-law with exponential decay, compared with Time 2 and the data derived from HCs.

16In order to quantify change in the distributions from Time 1 to Time 2, the stretched exponential function was fitted to each, using a range from degree k = 30 to 80. The resultant averaged fitting value of was 3.65±0.41 for Time 1 and 4.58±0.18 for Time 2 (Table 1). The value of was 4.78±0.16 for Control group. The change in these distributions was significant between time points and Time 2 more closely approximating the shape of that observed in healthy adults.

17We examined changes in functional connectivity at separate time points during a known period of recovery after TBI,. To do so, we employed graph theory analyses to document changes in resting BOLD signal between time points. To the best of our knowledge, this is the first report examining functional connectivity during recovery from neurological insult using graph theory.

18Un-weighted networks here revealed no significant changes between time points in neural networks during recovery. There is compelling evidence, however, that networks “weighted” according to a partial correlation matrix maintain important advantages over “un-weighted” data for understanding network characteristics. That is, weighted networks maintain information not only about suprathreshold connections, but the magnitude of those connections which potentially offers greater sensitivity to network change during a period of functional recovery. For these reasons, the following discussion focuses on the results of weighted network analysis.

19The current study analyzed “resting” data (as opposed to BOLD data associated with a known task) that are the result of a WM task. To provide context for understanding the resting connectivity data analyzed here, the “activation” findings achieved using a block design task revealed decreased involvement of relevant neural substrates (in particular prefrontal cortex) from Time 1 to time 2,. The findings here are consistent with what we have also observed during tasks of rapid decision making in TBI where task acclimation and faster reaction times result in diminished task-induced activation.

20In a growing functional imaging literature examining deficits in clinical samples, there remains little work examining how neurological insult influences connectivity and information transfer between nodes. That is, important questions remain regarding how neural networks respond to disruption and how plasticity is expressed in these networks during recovery. Overall, the results here indicate that recovery from TBI is associated with three important changes in resting networks: 1) a reduction in the number of highly significant connections in the overall network (strength S), or reduced “cost”, but little to no change in the number of connections, or sparsity; 2) a transition of the degree distribution from a power law with exponential decay at three months post injury to a stretched exponential decay at six months post injury; 3) increasing “small-worldness” during recovery, but residually diminished small-worldness compared to the HC sample. These findings are highlighted below.

21The degree distribution is a fundamental component of the indices that characterize a small-world network. Again a “degree” is the connectedness of any node to the network and network distributions take any of several forms including Gaussian, which is commonly observed in random networks or non-Gaussian forms including distributions with a fat tail (e.g., scale free networks). In the current study, the distribution more closely approximated power-law with exponential decay at Time 1 in the TBI sample, indicative of greater flexibility and reduced restriction on network change. This distribution shifted to becoming a stretched exponential decay distribution by Time 2 which more closely approximated what was observed in HCs (see Fig. 5).

22Thus, changes from Time 1 to Time 2 reflect a change from power-law with exponential decay to a distribution indicative of a more restrictive network, such as stretched exponential decay, which maintains greater constraint on overall network plasticity. From the viewpoint of statistical physics, such a transition in the tail of the degree distribution during recovery may be related to critical phenomena, and, indeed, recent work has reported observing characteristics consistent with criticality in the brain,,. Thus, this change in the degree distribution during recovery from neurological insult may operate to: 1) reduce “cost” in the network, 2) formalize the network, reducing its malleability and maximizing its efficiency. Earlier after injury, a distribution that more closely approximating power law might work to permit greater plasticity in order to accommodate challenge/disruption in the network. This early brain state may be neurally “expensive”, however, later giving way to a reduction in the number of highly significant connections and the overall strength of the network and increased small-worldness after recovery (i.e., Time 2). Overall, the current findings demonstrate that the degree distribution may be informative regarding the changing flexibility in global brain states.

23In “activation” studies in the clinical working memory literature a nearly universal finding is that neurological insult results in recruitment of additional neural resources during periods of cognitive challenge (see ). In the current resting data, the change occurs in the context of consistent “sparsity”, indicating that the magnitude of the relationships between nodes is reduced but the number of connections remains stable. If it can be assumed that not only the number but also the nature of the connections is essentially maintained from Time 1 to Time 2, it appears that existing support networks are recruited to tolerate network disruption and that the demand for these resources diminishes after recovery. These findings complement what is observed in cross-sectional studies of TBI,, and in healthy adults during periods of high task load – whereby increased involvement of neural resources may be highly dependent upon performance and does not represent fundamental changes in the network. That is, the alterations in neural connections observed during recovery may not signify formal brain reorganization (e.g., creation of novel connections). Instead these changes represent utilization of existing support resources early after neural disruption and this demand on auxiliary resources later remits resulting in a less costly network with greater neural efficiency.

24The current data demonstrate that early after injury, resting networks may be composed of more highly significant connections, creating a network higher in local clusters with shorter path length and a degree distribution that approximates the power law with cut-off. One characteristic of this degree distribution is that it permits greater neural dynamics and flexibility for adjusting the strength of connections within the network. One interpretation is that, over time, recovery results in resting connectivity that is diminished in overall strength and increased small-worldness approximating what is observed in healthy adults. Thus, recovery following TBI does not appear to result in elaboration and permanent inclusion of additional connections to the network, which may have important implications for understanding the neural recruitment observed in activation studies to date.

25While the current findings are the first to use graph theory and weighted network analysis to examine neural recovery in humans, this report is not without limitations. The most important shortcoming to the current study is the limited sample size, which does influence the generalizeability of the findings. Future work should endeavor to increase the sample size so that sample subtypes may be examined (e.g., injury severity, magnitude of behavioral improvement). Even so, there are two reasons why we believe the sample size here can support the conclusions made in this study. First, the effect of time on network analysis was large and, therefore observable (and statistically significant) within this small sample size. Second, for all indices, meaningful change from Time 1 to Time 2 resulted in networks more closely approximating those observed in the healthy adult sample. In sum, the current findings represent an important first look at how networks adapt to disruption during a critical window of recovery.

26Participants recruited for the study included 8 healthy adult control participants (HCs) (3 females, ages 19–51) and 8 individuals (3 females) sustaining severe TBI between the ages of 19 and 55 underwent MRI scanning at separate time-points. Two subjects with TBI were later excluded following data collection due to significant head movement (n = 1) and poor compliance with task demands (n = 1). Clinical descriptive information regarding patient injury is available in Table 3. The severity of TBI was defined using the Glasgow Coma Scale (GCS) in the first 24 hours after injury and GCS scores from 3–8 were considered “severe”. The BOLD fMRI data were collected at approximately 3 months and 6 months after resolution of posttraumatic amnesia; these 3-month and 6-month follow-up time-points are henceforth referred to as Time 1 and Time 2. Of note, subjects were selected based upon observable injury to prefrontal cortex which is known to subserve WM dysfunction, (see Fig. S1), thus providing the opportunity to directly examine the influence of disruption in neural networks related to the task. Prior to official enrollment in the study, all subjects signed IRB-approved informed consent forms approved by the Human Subjects Projection Office (HSPO) at Hershey Medical Center (HMC) in Hershey, PA. All study procedures complied with the HSPO office at HMC and with The Health Insurance Portability and Accountability Act (HIPAA) standards.

27This study employed a simple non-verbal working memory task requiring delayed-response scenarios. This task primarily requires sustained visuo-spatial attention, working memory rehearsal, and speeded processing. The detail information regarding the cognitive task and data analysis methods are identical to those presented previously (see ). In brief, the task requires memory for included exposure to one, two, or four black and white images of male and female faces. This paradigm required the subject to examine a target slide for 3000 msecs. After a delay of 3000 msecs requiring focus on a fixation point, the subject is provided a final target stimulus, or a single face presented in one of four quadrants. At the time of presentation of the Target, the subject is required to make a yes/no decision about the identity (match/no match) and the location of a single face presented in one of the four quadrants. In this delayed response task, the subject is required to examine a Stimulus Slide for 3000 msecs and after a delay of 3000 msecs (fixation point), the subject made a determination regarding a final target stimulus (3000 msecs).

28For all scans, data were acquired using a Philips 3T system (Philips Medical Systems, Best, the Netherlands) in the Department of Radiology, Hershey Medical Center, Hershey, PA. High resolution anatomical images (MPRAGE) with isotropic spatial resolution of 1.2 mm×1.2 mm×1.2 mm were acquired. Other imaging parameters of the MPRAGE sequence consisted of: 468.45 ms/16.1 ms/18, TR/TE/flip angle, 250×200 mm2 FOV, and a 256×180 acquisition matrix. Imaging parameters for echo planar imaging (EPI) consisted of: 2000 ms/30 ms/89, TR/TE/flip angle, mm2 FOV, 256×180 acquisition matrix. In this block design paradigm, 166 images were obtained for each run for a total scanning time of 6 minutes and 32 seconds.

29AFNI was used to perform the following preprocessing steps, including motion correction, spatial smoothing, and mean-based intensity normalization. As a final preprocessing step, each individual's time series was spatially normalized by registration to the MNI152 template (Montreal Neurological Institute), with 2 mm3 resolution, using a 12 degrees-of-freedom affine transformation. Parcellation divides each cerebral hemisphere into 56 anatomical regions of interest (ROI). Regional mean time series were estimated for each individual by averaging the fMRI time series over all voxels in each of 112 regions. The resultant signals ware filtered with a bandpass filter ([0.01, 0.1] Hz).

30Recently, Fox et. al reported that task-evoked activity are linearly superimposed on underling spontaneous activity (resting data),. This funding suggests that if the effect of task-induced activity is adequately removed from task-related design data, the remaining residual component should represent resting state data. The deterministic task-related brain activity can be properly modeled using a linear model, and can be removed by regression approaches. Here, we employed an orthogonalization method to remove task-related effects from intrinsic BOLD fluctuations.

31Functional connectivity is defined as the temporal correlations between spatially remote neurophysiological events, –. The correlation between a given pair of regions is typically used to index functional connectivity,,, and the resultant pair-wise correlation matrices R are then thresholded to generate a statistically significant functional connectivity network.

32In this study, the partial correlation coefficients were used as a measure of functional connectivity between given pairs of BOLD signal from regions i and j (i, j = 1,2,…,M; here M is a number of ROIs). The partial correlation matrix R is a symmetric matrix in which each off-diagonal element is the correlation coefficient between a pair of regions after filtering out effects of all other brain regions. The procedure for obtaining partial correlation values here is consistent with other work,.

33For the current data, construction of a functional connectivity network was achieved by testing the null hypothesis that the partial correlation coefficient was significantly different from zero between any region pairs. Two regions were considered functionally connected if their partial correlation coefficient was significant at p<0.05 level, resulting in an adjacency matrix A containing information about the connectivity structure of the network. From resultant matrix, the sparsity, which is the connection density in the network, and the strength S, which is an average of all significant correlation coefficients and a measure of expensiveness in the network, were evaluated.

34For the purposes of this study, both weighted and un-weighted networks were examined. In constructing un-weighted networks, each component of the adjacency matrix was labeled 1 if a corresponding was significant, and 0 in other cases. For weighted networks, a number of potential “weight” definitions can be used. Consistent with previous work using graph theory to examine neural systems, networks were defined with the definition if the corresponding was significant and if it was non-significant. Such analysis results in adjacency matrices: graph consisting of a set of nodes (brain regions), a set of edges (significant connectivities), and a set of associating weights between node i and j. Of note, alternative weight definitions such as have been applied to examine other networks (e.g., world wide web),. This weight definition was used in the current data and the results were largely similar to those presented in this manuscript. We have chosen to focus here on the results using as a weight definition because of its prior use in neural systems.

35In Graph theory, the clustering coefficient C, and characteristic path length L of any graph G are key indices characterizing the statistical properties of a network. For weighted networks the shortest path length between two nodes i and j is the smallest sum of the weights for all the possible paths from node i to node j. The path length L of a graph G is the mean of the shortest path lengths over all possible pairs of nodes: L thus represents the extent and average connectivity or overall routing efficiency of a network.

36The clustering coefficient C i of a node i with degree is calculated as the geometric average of subgraph weights:where is the normalized weight by the largest value of the weight in the network. The clustering coefficient C of the network is an average of C i over all nodes. The clustering coefficient is an index of local structure, and has been interpreted as a measure of resistance to failure and the extent of the local density or cliquishness of the network.

37In random networks, characteristic path length is typically short, and clustering coefficient is typically small,. In small-world networks, by definition, we expect the ratio and the ratio. A scalar summary of small-worldness can be defined as the ratio, which is typically >1 for a small world network, or greater than random clustering and near random path length,. In order to obtain random networks, we used a surrogate method where generated networks preserve the same number of nodes, edges, and degree (the degree of each node is the number of nodes directly connected to the node i distribution as the real network). For all data sets here, including each session per subject, 100 random networks were generated for each network and the averaged indices and were calculated.

38Application of Graph theory also provides the opportunity to examine network efficiency on global and local scales,,,. The global efficiency of a graph G is defined as;The global efficiency is an indicator of global efficiency of parallel information transfer in the network,,.

39The local efficiency of each node can be defined as follows:where G i is a set of nodes directly connecting to the ith node,,,. The measure is a mean local efficiency of information transfer in the immediate neighborhood of each node. Note that both and can be calculated for both un-weighted and weighted networks as well as path length and clustering coefficient.

40The degree distribution likely represents the most fundamental component of the various indices that characterize a graph. The degree of each node is defined as the number of nodes directly connected to the node i in the network, so the degree distribution is the observed distribution of the number of direct connections of nodes in the network. Here we use a cumulative distribution in the place of. The degree distribution furthermore classifies small-world networks into different structural classes including: (a) scale-free networks, characterized by a degree distribution that decays as a power law, i.e, (b) broad-scale networks, characterized by a degree distribution that has a power law regime followed by a sharp cutoff, and (c) single-scale networks characterized by a degree distribution with a fast decaying tail.

41Wilcoxon matched-pairs test was used to test significant differences between different time points. For Control group there were no significant differences in any network indices between point Time 1 and Time 2, thus we used the averaged values of both points for all indices to increase reliability in the observed networks. Dannet's test was used to determine the differences in network indices between the TBI sample and the Healthy Control (HC) sample. In the data presented here, error bars indicate the standard error of mean and an asterisk in the graphs indicates a significant difference between time points.