It has been previously shown that traumatic brain injury (TBI) is associated with reductions in metastability in large-scale networks in resting-state fMRI (rsfMRI). However, little is known about how TBI affects the local level of synchronization and how this evolves during the recovery trajectory. Here, we applied a novel turbulent dynamics framework to investigate whole-brain dynamics using an rsfMRI dataset from a cohort of moderate to severe TBI patients and healthy controls (HCs). We first examined how several measures related to turbulent dynamics differ between HCs and TBI patients at 3, 6, and 12 months post-injury. We found a significant reduction in these empirical measures after TBI, with the largest change at 6 months post-injury. Next, we built a Hopf whole-brain model with coupled oscillators and conducted in silico perturbations to investigate the mechanistic principles underlying the reduced turbulent dynamics found in the empirical data. A simulated attack was used to account for the effect of focal lesions. This revealed a shift to lower coupling parameters in the TBI dataset and, critically, decreased susceptibility and information-encoding capability. These findings confirm the potential of the turbulent framework to characterize longitudinal changes in whole-brain dynamics and in the reactivity to external perturbations after TBI.

Significant advances in nonlinear dynamics and computational modeling have opened up the possibility of studying how whole-brain dynamics may be impacted by traumatic brain injury (TBI). One hypothesis is that whole-brain dynamics may show differential spatiotemporal patterns during the recovery trajectory. We used a novel turbulence framework to examine how several measures related to turbulent dynamics differ between healthy controls and TBI patients at 3, 6, and 12 months post-injury. This framework revealed a significant reduction in these empirical measures after TBI differentially affecting long distances, with the largest change at 6 months post-injury. The alterations observed at network level, however, showed a certain degree of recovery after 1 year. In addition, the Hopf whole-brain model demonstrated decreased susceptibility and information-encoding capability after TBI. The clinical implications of this work are discussed.

There has been a growing interest in the development and application of large-scale, nonlinear brain modeling methods to resting-state functional magnetic resonance imaging (fMRI) data (Breakspear, 2017; Cabral et al., 2014; Deco et al., 2015; Deco et al., 2018; Deco et al., 2019; Deco et al., 2021; Deco & Kringelbach, 2014; Kringelbach & Deco, 2020), as it offers whole-brain coverage and can be acquired in patient populations who are unable to perform cognitive tasks (Demertzi et al., 2015; Escrichs et al., 2022; López-González et al., 2021). Importantly, these methodological advancements have enabled reconceptualizing existing measures of nonlinear dynamics in order to explain the complex dynamics of human brain activity. One prominent example is turbulence, which has been recently demonstrated in the healthy brain (Cruzat et al., 2022; De Filippi et al., 2021; Deco & Kringelbach, 2020) and in patients with a disorder of consciousness following brain injury (Escrichs et al., 2022). In this study, we combine a turbulent dynamics framework (model-free) with a computational modeling approach (model-based) to characterize whole-brain dynamics in moderate to severe traumatic brain injury (TBI) patients during a 1-year recovery period.

TBI is associated with a wide range of cognitive deficits, such as attentional or memory impairments or executive dysfunction, that are often persistent, contribute to poor functional outcomes, and have a profound impact on the quality of life of the patients (McAllister et al., 2004; Ponsford et al., 2014). The emergence of cognitive function in the brain is contingent on the integrity of the structural connectome, especially long-range white matter connections that contribute to the small-world topology underlying the spatiotemporal dynamical patterns observed in a healthy regime (Bassett & Bullmore, 2006). In addition to focal lesions, diffuse axonal injury has been considered one of the main mechanisms of TBI (Jolly et al., 2021), leading to disturbances in long-range connections and ultimately disrupting the spatiotemporal properties of functional brain networks (Caeyenberghs et al., 2014; Hellyer et al., 2013; Jilka et al., 2014; Kinnunen et al., 2011). Such structural and functional disruption results in long-term cognitive impairment (Bonnelle et al., 2011; Jilka et al., 2014; Kinnunen et al., 2011) and may play a role in the progressive neurodegeneration associated with TBI (Sharp et al., 2014). Nevertheless, reliable biomarkers accurately predicting cognitive outcomes after TBI have not been recognized (National Academies of Sciences, Engineering, and Medicine, 2022). Here, resting-state fMRI (rsfMRI)-based measures like turbulence could be clinically relevant as they can be assessed repeatedly without a specific cognitive task (De Filippi et al., 2021; Deco & Kringelbach, 2020; Escrichs et al., 2022). Despite its clinical relevance, a better understanding of how TBI affects whole-brain dynamics longitudinally is still needed.

The shift towards describing the impact of brain injury in terms of changes in whole-brain dynamics holds great promise to bestow a common framework not only to characterize cognitive deficits after TBI but, more generally, to provide insights into the mechanistic principles that support brain function. One widely used approach to studying whole-brain dynamics is based on metastability, a dynamical regime where brain regions engage and disengage flexibly over time, showing transient dynamical patterns (Friston, 1997; Shanahan, 2010; Tognoli & Kelso, 2014). In the context of TBI, recent evidence has shown that patients had reduced global and network-level metastability compared with controls (Hellyer et al., 2015). Critically, this reduction was associated with cognitive impairment and damage to structural connectivity, suggesting that metastability may act as a conceptual bridge between brain structure and behavior (Tognoli & Kelso, 2014). However, the concept of metastability measured as the variability over time of the global level of synchronization of the brain, commonly known as the global Kuramoto order parameter of a dynamical system (Kuramoto, 1984), is limited as it cannot dissociate the dynamical patterns occurring at local spatial scales.

Recently, we have proposed an extension of metastability to account for such limitation based on turbulence theory developed first for fluid dynamics and since extended to other domains. As shown by Kuramoto (Kuramoto, 1984), coupled oscillators can be used to describe synchronization globally through metastability (Cabral et al., 2014; Deco, Kringelbach, et al., 2017) and to capture turbulence using a measure of local metastability (Kawamura et al., 2007). By combining Kuramoto’s framework with Kolmogorov’s concept of structure functions (Kolmogorov, 1941), turbulence can be calculated locally at different spatial scales in an analogy to the rotational vortices found in fluid dynamics. By taking this approach, we have found long-range correlations in higher order task-specific regions that could be controlling a functional resting-state core involving primary sensory regions (Deco & Kringelbach, 2020). Further, the results from the whole-brain computational model indicated that turbulence is a good index of the information-processing capability of the brain, as both reached the maximum level at the optimal fitting of the model. When applied to a population of patients with disorders of consciousness—which is often altered in severe cases of TBI—turbulence was significantly reduced at large and intermediate distances, as was the sensitivity of the brain to respond to external perturbations (Escrichs et al., 2022). Therefore, the question arises as to whether similar alterations in whole-brain dynamics could be found after TBI, across different stages of recovery.

Here, by combining empirical and computational approaches, we investigated how whole-brain dynamics are affected by TBI measured at three time points (3, 6, and 12 months post-injury) across 1 year. To do so, we used a publicly available rsfMRI dataset of patients with moderate to severe TBI (N = 14 in each session) and age-matched healthy individuals (N = 12) and extracted the time series from each of the 1,000 parcels in the fine-grained Schaefer parcellation (Schaefer et al., 2018). We tested the hypotheses that whole-brain amplitude turbulence would (a) accurately differentiate patients with TBI from healthy controls (HCs) and (b) progressively approach the values of HCs across the 1-year recovery trajectory. Furthermore, we expected to see other measures derived from the turbulent dynamics framework to discriminate between TBI patients and HCs. Finally, we anticipated that the whole-brain computational model would demonstrate that the brain’s capability to respond to external in silico perturbations is compromised after TBI.

Model-Free Approach

Participants and behavioral data.

The source of data for this study (Figure 1A) was obtained from the open-access repository OpenNeuro (https://openneuro.org/datasets/ds000220/versions/00002). Please refer to the original publication (Roy et al., 2017) for a full description of the participants included. Briefly, the dataset included 14 patients with moderate to severe TBI between the ages of 18 and 36 with education ranging from 12 to 18 years, and 12 HCs of comparable age and education. A Glasgow Coma Scale at time of injury was indicative of severe (3–8) and moderate (9–12) injury. The NIfTI images from 2 TBI patients were discarded because of incorrect labels (superior-inferior swap) as shown in FSLeyes. All subjects with TBI completed three separate sessions (3, 6, and 12 months following injury). Resting-state fMRI data from HCs at baseline was also included for group comparisons. The authors of the original publication also provided the z scores for the neuropsychological battery used to evaluate TBI patients at 3 months post-injury. These included the Visual Search and Attention Test, WAIS-III Digit Span (Forward and Backward), Trail Making Test B, and Stroop Color and Color-Word.

Figure 1.

Schematic of the model-free and model-based analysis pipelines. We used an open-access longitudinal dataset from a previous publication (Roy et al., 2017). T1-weighted (T1-W) and resting-state fMRI (rsfMRI) data from traumatic brain injury (TBI) patients were collected at three time points (3, 6, and 12 months post-injury). For group comparisons, T1-W and rsfMRI were collected from age-matched healthy controls in two sessions. (A) Model-free approach. Top left: rsfMRI data were preprocessed in the CONN toolbox and the time series were extracted for each of the 1,000 nodes in the Schaefer parcellation and the instantaneous phase calculated (Schaefer et al., 2018). Bottom left: Visualization of the change over time and space of the Kuramoto local order parameter in a healthy subject (left hemisphere). The Kuramoto local order parameter can define the turbulent signatures of brain activity at different spatial scales in the vortex space. Here, we computed four turbulent measures based on the Kuramoto local order parameter to characterize the brain’s information processing: amplitude turbulence, information cascade flow, information cascade, and information transfer (see the Methods section for more details). (B) Model-based approach and simulated attack. We built Hopf whole-brain dynamical models that described the intrinsic dynamics of each brain area by a Stuart-Landau nonlinear oscillator. Such local oscillators are linked through the underlying structural connectivity to simulate the global dynamics. We used the diffusion tensor imaging matrix with the exponential distance rule and long-range exceptions as structural connectivity to fit the empirical data as a function of the Euclidean distance. To account for the effect of focal lesions, we adapted a simulated attack approach (Medaglia et al., 2022). The overlap between the patient’s specific lesion mask and the Schaefer parcellation was calculated and used to create a group lesion mask thresholded at 1.5 and 2 SD from the group’s overlap. Next, a binary or weighted disconnection was applied to the group lesion mask. For each threshold and disconnection approaches, a model was computed by multiplying the structural connectivity by the group lesion mask. Then we applied in silico perturbations to assess the model’s reaction for each time point in the TBI group compared with healthy controls using the susceptibility and information-encoding capability model-based measures (bottom panels adapted with permission from Deco & Kringelbach, 2020).

Figure 1.

Schematic of the model-free and model-based analysis pipelines. We used an open-access longitudinal dataset from a previous publication (Roy et al., 2017). T1-weighted (T1-W) and resting-state fMRI (rsfMRI) data from traumatic brain injury (TBI) patients were collected at three time points (3, 6, and 12 months post-injury). For group comparisons, T1-W and rsfMRI were collected from age-matched healthy controls in two sessions. (A) Model-free approach. Top left: rsfMRI data were preprocessed in the CONN toolbox and the time series were extracted for each of the 1,000 nodes in the Schaefer parcellation and the instantaneous phase calculated (Schaefer et al., 2018). Bottom left: Visualization of the change over time and space of the Kuramoto local order parameter in a healthy subject (left hemisphere). The Kuramoto local order parameter can define the turbulent signatures of brain activity at different spatial scales in the vortex space. Here, we computed four turbulent measures based on the Kuramoto local order parameter to characterize the brain’s information processing: amplitude turbulence, information cascade flow, information cascade, and information transfer (see the Methods section for more details). (B) Model-based approach and simulated attack. We built Hopf whole-brain dynamical models that described the intrinsic dynamics of each brain area by a Stuart-Landau nonlinear oscillator. Such local oscillators are linked through the underlying structural connectivity to simulate the global dynamics. We used the diffusion tensor imaging matrix with the exponential distance rule and long-range exceptions as structural connectivity to fit the empirical data as a function of the Euclidean distance. To account for the effect of focal lesions, we adapted a simulated attack approach (Medaglia et al., 2022). The overlap between the patient’s specific lesion mask and the Schaefer parcellation was calculated and used to create a group lesion mask thresholded at 1.5 and 2 SD from the group’s overlap. Next, a binary or weighted disconnection was applied to the group lesion mask. For each threshold and disconnection approaches, a model was computed by multiplying the structural connectivity by the group lesion mask. Then we applied in silico perturbations to assess the model’s reaction for each time point in the TBI group compared with healthy controls using the susceptibility and information-encoding capability model-based measures (bottom panels adapted with permission from Deco & Kringelbach, 2020).

Close modal

Resting-state fMRI acquisition.

Five-minute rsfMRI data were acquired using a whole-brain EPI sequence (150 volumes, TR = 2,000 ms, TE = 30 ms, flip angle = 90°, 240 × 240 mm field of view, 80 × 80 acquisition matrix, voxel size 1.9 × 1.9 × 4 mm). In addition, a high-resolution whole-brain T1-weighted image was acquired (voxel size 1 mm3).

Resting-state fMRI preprocessing.

Preprocessing was conducted using Statistical Parametric Mapping software (SPM12, Wellcome Department of Cognitive Neurology, University College London) running under MATLAB Release 2021a (MathWorks, Inc., Natick, MA). The first five volumes were removed to control for signal instability. Preprocessing steps included slice timing correction, realignment, segmentation, functional to structural coregistration, spatial normalization (voxel size 2 mm3), and smoothing (8 mm Gaussian kernel).

We then applied a denoising pipeline in order to minimize the variability due to physiological, outlier, and residual subject-motion effects (Nieto-Castanon, 2020). The preprocessed output files from SPM were imported into the CONN toolbox version 20.b (https://www.nitrc.org/projects/conn), including the following subject- and session-specific files: T1 and rsfMRI scans; segmented gray matter (GM), white matter (WM), and cerebrospinal fluid (CSF) images; and six realignment parameters after rigid body motion correction. We applied the same pipeline as in our previous publication with TBI patients (Martínez-Molina et al., 2021; see Supporting Information for more details). The 1,000-node and 400-node Schaefer parcellations in MNI152 volumetric space were added as an atlas in the CONN toolbox to extract the time series from each node (Schaefer et al., 2018). The Euclidean distance between parcels’ centroids from the Schaefer parcellation was also calculated as in our previous publications (Deco & Kringelbach, 2020; Escrichs et al., 2022).

Probabilistic tractography analysis.

We used the diffusion-weighted and T2-weighted neuroimaging data from 32 participants of the HCP database, as previously reported (Deco & Kringelbach, 2020). In the HCP website, a description of the acquisition parameters is detailed (Setsompop et al., 2013). We used the Lead-DBS software package (https://www.lead-dbs.org/) for the preprocessing of tractography data (Horn et al., 2017), which, briefly, were processed by using a q-sampling imaging algorithm implemented in DSI Studio (https://dsi-studio.labsolver.org). Segmentation of the T2-weighted anatomical images and coregistering the images to the b0 image of the diffusion data using SPM12 produced a WM mask. For each individual, 200,000 fibers were sampled within the WM mask. Fibers were transformed into MNI space using the Lead-DBS software (Horn & Blankenburg, 2016). Finally, the structural connectivity (SC) using the Schaefer 1,000 parcellation (Schaefer et al., 2018) was extracted with the standardized methods in Lead-DBS. The SC was averaged across participants before the simulated attack.

Simulated attack.

In this approach, we used the binary lesion maps from the TBI patients provided by the authors of the original publication (Roy et al., 2017) in order to simulate connectome attacks relative to the SC derived from healthy participants (see previous section). Our goal was to create a lesion mask matrix based on the lesions in the TBI patients that could be applied to the healthy SC. A similar methodology has been recently implemented in patients with aphasia and provided a good approximation of the effects of the real lesions on anatomic networks as well as explaining a comparable amount of variance in behavioral data (Medaglia et al., 2022). Further, virtual lesions have allowed us to disentangle the local and recurrent components from transcranial magnetic stimulation (TMS)-EEG evoked potentials using a whole-brain connectome-based computational modeling (Momi et al., 2023).

The lesion masks for each TBI patient (N = 11, one patient did not have visible lesions) in native space were normalized to MNI space (after normalization one patient was discarded because of misalignment) and the overlap with each node from the Schaefer 1,000-node parcellation computed using the fslstats command in FSL v.6.0 (https://fsl.fmrib.ox.ac.uk/fsl/). This value was divided by the volume of the Schaefer node to account for variability in the size of the nodes. We created the lesion mask matrix to simulate the lesion on the healthy SC using four thresholds: 1.5, 2, 3, and 4 SD from the mean overlap volume of lesions in all TBI patients. Thus, for each TBI patient, we calculated the nodes to attack based on these four thresholds. This approach allowed us to calculate the frequency of lesion for a given node as well as find the unique nodes to attack, that is, nodes that were lesioned in any given patient. The lesion mask matrices thresholded at 3 and 4 SD resulted in only 10 and 7 nodes to attack and were discarded. The lesion mask matrices thresholded at 1.5 and 2 SD produced 41 and 22 nodes to attack. The simulated attack with the 1.5 and 2 SD lesion matrices was performed in two ways, leading to the disconnection of the lesioned node: (a) a complete deletion of the SC between the node and the rest of the brain (binary approach, value set to 0 in the lesion mask matrix); (b) a weighted downregulation of the SC between the node and the rest of the brain. The weight was computed by calculating the number of patients with the lesioned node divided by the total number of TBI patients. Then the value in the lesion mask matrix between that node and the rest of the brain was set to 1 - weight. For each masking threshold, the SC was multiplied by the lesion mask matrix before running the computational modeling. We took a dual approach, one binary (1) and one weighted (2), in order to ascertain the sensitivity of the whole-brain model perturbational measures to the effect of virtual disconnections. A list of the nodes to attack for the 1.5 SD threshold, the corresponding anatomical region, and the frequency of lesion can be found in Table S19 in the Supporting Information.

Amplitude turbulence.

We measured amplitude turbulence by first defining the Kuramoto local order parameter and then taking the standard deviation of the modulus across time and space. First, we defined the amplitude turbulence, Rλ (x¯, t), as the modulus of the Kuramoto local order parameter for a given brain region as a function of the time at each spatial scale λ:
Rλx¯teiϑλx¯t=kdx¯Gλx¯x¯eiφx¯t,
(1)
where φ(x¯, t) denotes the phases of the BOLD signal data, Gλ refers to the local weighting kernel Gλx¯=eλx¯, and k denotes the normalization factor dx¯Gλx¯x¯1. The BOLD time series were filtered with a second-order Butterworth filter in the range between 0.008 and 0.08 Hz, where we chose the typical high-pass cutoff to filter low-frequency signal drifts (Fox et al., 2005; Glerean et al., 2012) and the low-pass cutoff to filter the physiological noise, which tends to dominate the higher frequencies (Fox et al., 2005; Glerean et al., 2012). We then applied the Hilbert transform in order to obtain the complex phase of the BOLD signal for each brain node as a function of time.
Local levels of synchronization at a certain λ, as a function of space (x¯) and time (t), are defined by Rλ. We then measure the amplitude turbulence, D, as the standard deviation across nodes and time of the modulus of the Kuramoto local order parameter (Rλ):
D=Rλ2x,tRλx,t2,
(2)
where the brackets 〈〉x,t denote averages across space and time.

This measure captures the brain vortex space over time, motivated from fluid dynamics and the rotational vortices (Deco & Kringelbach, 2020). Specifically, we explored the level of Kuramoto amplitude turbulence over different λ values, that is, from 0.01 (∼100 mm) to 0.30 (∼3 mm), in steps of 0.03.

Information cascade flow.

The information cascade flow indicates how the information travels from one scale (λ) to a lower scale (λ − Δλ, where Δλ corresponds to a scale step) in successive time steps (t and t + Δt). It is calculated as the temporal correlation between the Kuramoto local order parameter in two successive scales and times as the following:
Fλ=corrtRλx¯t+ΔtRλΔλx¯tx¯,
(3)
where the brackets x¯ indicate averages across time and nodes. Then the information cascade is obtained by averaging the information cascade flow across scales λ, which captures the whole behavior of the information processing across scales.

Information transfer.

The information transfer captures how information travels across space at each λ scale. It is calculated as the slope of a linear fitting, in the log-log scale, of the temporal correlation between the Kuramoto local order parameter of two nodes at each scale as a function of its Euclidean distance (r) in the inertial subrange.
logcorrtRnλRpλr=Aλ*logr+Bλ,
(4)
where Aλ and Bλ are the fitting parameters for the λ scale, and r denotes the distance in the brain.

Local node-level metastability.

The node-level metastability, NLM, defines the variability of the local synchronization of the nodes and is computed as the standard deviation across time of the Kuramoto local order parameter as follows:
NLMnλ=Rnλt2tRnλtt2,
(5)
where brackets 〈〉t denote average across time.
Here, we estimated the discrete version of the node-level Kuramoto order parameter, with modulus R and phase ν, which represents a spatial average of the complex phase factor of the local oscillators weighted by the coupling:
Rnλteiνnt=pCnpλqCnqλeiφpt,
(6)
where φp(t) are the phases of the spatiotemporal data and Cnqλ is the local weighting kernel between node n and q, and λ defines the spatial scaling:
Cnp=eλrnp.
(7)

Model-Based Approach

Whole-brain Hopf model.

We built whole-brain dynamical models with Stuart-Landau oscillators based on the normal form of a supercritical Hopf bifurcation (Deco, Kringelbach, et al., 2017). This type of bifurcation can change the qualitative nature of the solutions from a limit cycle that yields self-sustained oscillations toward a stable fixed point in phase space. This model is characterized by model parameters that rule the global dynamical behavior. One control parameter is the multiplicative factor (G) representing the global conductivity of the fibers scaling the structural connectivity between brain areas, which is assumed to be equal across the brain (Deco, Kringelbach, et al., 2017; Deco, Tagliazucchi, et al., 2017). The other relevant parameters are the local bifurcation parameter (an), which rules the dynamical behavior of each area between noise induced (a < 0), self-sustained oscillations (a > 0), or a critical behavior between both (a ∼ 0) (Figure 1B). The model parameters were optimized to better fit the empirical functional connectivity as a function of the distance, r, within the inertial subrange. The models consisted of 1,000 brain regions from the resting-state atlas (Schaefer et al., 2018). The underlying structural connectivity matrix Cnp was added to link the brain structure and functional dynamics by computing the exponential distance rule (Equation 7). The local dynamics of each brain region were characterized by the normal form of a supercritical Hopf bifurcation, which emulates the dynamics for each region from noisy to oscillatory dynamics as follows:
dxndt=anxnxn2+yn2xnωnyn+νηnt,
(8)
dyndt=anynxn2+yn2yn+ωnxn+νηnt,
(9)
where ηn(t) is the additive Gaussian noise with standard deviation ν = 0.01. This normal form has a supercritical bifurcation at an = 0, such that for an > 0, the system is in a stable limit cycle oscillation with frequency fn = ωn/2π, while for an < 0 the local dynamics are in a stable point (the noisy state). The frequency ωn of each region was obtained from the empirical functional MRI data as the power spectrum peak.
Finally, the whole-brain dynamics were determined by the set of coupled equations as follows:
dxndt=anxnxn2+yn2xnωnyn+Gp=1NCnpxptxn+νηnt,
(10)
Dyndt=anynxn2+yn2yn+ωnxn+Gp=1NCnpyptyp+νηnt,
(11)
where Cnp is the structural connectivity weight between node n and p, and G is the global coupling weight with equal contribution between all the nodal pairs. The Cnp was derived using the exponential distance rule and long-range exceptions (Deco et al., 2021). We used λ = 0.18, as this was shown to produce the best fit. This parameter is different from the λ values used to calculate global amplitude turbulence at different spatial scales. For each simulated attack (four in total: 1.5 and 2 SD, binary and weighted approach), this Cnp was multiplied by the corresponding lesion mask matrix. The output of the model is a set of complex-valued BOLD-like signals with time-varying intrinsic frequencies (ω). The real part is considered to represent BOLD signal recorded in the fMRI session, whereas the imaginary part characterizes the hidden state of the oscillator unseen to the scanner.

Functional connectivity fitting.

The functional connectivity fitting of the BOLD signal data was assessed by applying Kolmogorov’s structure-function. The structure-function for a variable u characterizes the evolution of the functional connectivity (FC) as a function of the Euclidean distance between equally distant brain regions and is defined as follows:
Sr=ux¯+rux¯2x,t=2FC0FCr,
(12)
where FC(r) is the spatial correlations of two points separated by a Euclidean distance r, and is given by
FCr=ux¯+rux¯x,t,
(13)
where 〈〉x,t denotes average across the spatial x coordinates of the nodes and time. Then, we calculated the fitting as the Euclidean distance between simulated and empirical FC(r) within the inertial subrange (Deco & Kringelbach, 2020).

Susceptibility and information-encoding capability.

The susceptibility measure obtains the brain’s sensitivity to react to external stimulation. The Hopf model was perturbed for each G by randomly changing the local bifurcation parameter, an, in the range [−0.02:0] and was computed by measuring the modulus of the Kuramoto local order parameter as follows:
χ=R~λsx¯ttRλsx¯tttrialsx¯,
(14)
where R~λsx¯t corresponds to the perturbed case, the Rλsx¯t to the unperturbed case, and 〈〉t, 〈〉trials, and x¯ to the average across time, trials, and space, respectively.
The information-encoding capability (I) captures how the external stimuli are encoded in whole-brain dynamics. This measure was defined as the standard deviation across trials of the difference between the perturbed R~λs(x¯, t) and unperturbed Rλs(x¯, t) mean of the modulus of the local Kuramoto order parameter across time t, averaged across all brain regions n as follows:
I=R~λsx¯ttRλsx¯tt2trialsx¯R~λsx¯ttRλsx¯tt2trialsx¯,
(15)
where the brackets 〈〉t, 〈〉trials, and x¯ denote the averages defined as above.

Statistical Analyses

A 4 (Group) × 10 (Lambda) ANOVA was used to compare TBI patients and controls for global amplitude turbulence and information transfer. Information cascade flow differences were assessed with a 4 (Group) × 9 (Lambda) ANOVA. One-way ANOVA (four groups) was used to test significant differences in information cascade. In all cases, post hoc tests were performed to evaluate pairwise differences. The longitudinal effects within the TBI group were assessed with repeated measures ANCOVA with lesion volume as a nuisance covariate followed by post hoc comparisons. Differences in resting-state network (RSN) turbulence were assessed with one-way ANOVA (four groups) followed by post hoc tests. Differences in node-level metastability were examined with the Kolmogorov-Smirnov test. Pearson’s correlations were performed to assess brain-behavior correlations. A Wilcoxon signed rank test was used to determine differences within HCs or TBI groups in susceptibility and information-encoding capability. Between-group differences were assessed with a Wilcoxon ranked sum test. When needed, we corrected for multiple comparisons using false discovery rate (FDR) with the method of Benjamini and Hochberg (Benjamini & Hochberg, 1995).

Our findings are split into two different approaches: model-free and model-based. In the first approach, we computed empirical measures of whole-brain turbulent dynamics using a fine-grained parcellation with 1,000 regional phase time courses derived from rsfMRI in both HCs (N = 12) and patients with TBI (N = 12) at 3, 6, and 12 months post-injury (see Figure 1A and the Methods section for more details). Since there were not significant differences between Session 1 and 2 in HCs (see Table S1 in the Supporting Information), we averaged the values for both sessions. In the model-based approach, we built a whole-brain dynamical model for each condition and group to explore the sensitivity of TBI patients and HCs to react to external in silico perturbations. For each whole-brain model, we applied in silico perturbations by changing the local bifurcation parameter and computed the susceptibility and information-encoding capability measures. A whole-brain model was computed separately for each stage of recovery in TBI patients and for each simulated attack approaches as well as with the intact structural connectivity (see Figure 1B and the Methods section for more details). Given the significant differences between Session 1 and 2 in HCs for the perturbational measures (Supporting Information Table S1), these were not averaged across sessions.

Model-Free Measures of Synchronization

Empirical measures of turbulence are reduced after traumatic brain injury.

TBI patients showed differences compared with HCs for each turbulent measure. To examine global amplitude turbulence differences between HCs and TBI patients at the 10 spatial scales under study, we used a 4 (Group) × 10 (Lambda) ANOVA. The results from this analysis revealed a significant Group × Lambda interaction (F(27,440) = 2.534, p < 0.0001; Supporting Information Table S2, Figure 2A). Post hoc tests indicated that HCs’ amplitude turbulence was significantly higher than that of TBI patients across all stages of recovery for λ = 0.01 and λ = 0.03 (Table S2), which are the scales including larger spatial distances across nodes. For λ = 0.06, we also found significant differences between HCs and TBI patients at 3 and 6 months post-injury (Table S2). It is worth noting that the greatest mean difference was observed between HCs and TBI patients at 6 months at all significant λ values. The rest of the λ values did not show any significant difference between HCs and TBI patients post hoc.

Figure 2.

Turbulent dynamic measures differ between traumatic brain injury (TBI) patients and healthy controls (HCs) over time, peaking at 6 months post-injury. Significant differences in turbulent-like dynamics between HCs, averaged across sessions, and TBI patients at 3, 6, and 12 months post-injury were calculated using a 4 (Group) × 10 (Lambda) ANOVA. The Group × Lambda interaction was significant at long distances for both (A) global amplitude turbulence (λ = 0.01–0.06) and (B) information cascade flow (λ = 0.01–0.03). (C) We found a group effect for information cascade. Post hoc tests revealed that it was due to the difference between HCs and TBI patients at 6 months post-injury. (D) For information transfer, the Group × Lambda interaction was significant for intermediate and short distances (λ = 0.27–0.15). Bar plots depict mean values (+SEM, standard error of mean). The lower and upper edges of the boxplots indicate the first (25th percentile) and third quartile (75th percentile), and middle lines are medians. Asterisks represent the effect of the interaction for all measures except for information transfer, where the results from the post hoc test are shown (*** p < 0.0001).

Figure 2.

Turbulent dynamic measures differ between traumatic brain injury (TBI) patients and healthy controls (HCs) over time, peaking at 6 months post-injury. Significant differences in turbulent-like dynamics between HCs, averaged across sessions, and TBI patients at 3, 6, and 12 months post-injury were calculated using a 4 (Group) × 10 (Lambda) ANOVA. The Group × Lambda interaction was significant at long distances for both (A) global amplitude turbulence (λ = 0.01–0.06) and (B) information cascade flow (λ = 0.01–0.03). (C) We found a group effect for information cascade. Post hoc tests revealed that it was due to the difference between HCs and TBI patients at 6 months post-injury. (D) For information transfer, the Group × Lambda interaction was significant for intermediate and short distances (λ = 0.27–0.15). Bar plots depict mean values (+SEM, standard error of mean). The lower and upper edges of the boxplots indicate the first (25th percentile) and third quartile (75th percentile), and middle lines are medians. Asterisks represent the effect of the interaction for all measures except for information transfer, where the results from the post hoc test are shown (*** p < 0.0001).

Close modal

The information cascade flow also captured differences in information propagation between consecutive spatial scales. We found a significant Group × Lambda interaction (F(24,396) = 3.007, p < 0.0001; Figure 2B). Like the findings with amplitude turbulence, post hoc tests revealed group differences at larger spatial scales. Specifically, all group comparisons were significant for λ = 0.01 and between HCs and TBI at 6 months for λ = 0.03 (Table S3 in the Supporting Information). In the one-way ANOVA, we found a significant group effect for information cascade (F(3,44) = 7.061, p < 0.001; Figure 2C). As indicated by post hoc tests, a significant group difference was found between HCs and TBI at 3 months (T(22) = 2.833, p = 0.034; Table S4 in the Supporting Information) and 6 months (T(22) = 4.552, p < 0.001; Table S4).

The Group × Lambda interaction was also significant for information transfer (F(27,440) = 2.325, p < 0.001; Figure 2D). In this case, however, differences were found at larger and intermediate λ values (λ = 0.27–0.15), with all group comparisons differing for λ = 0.27 and λ = 0.24 (Table S5 in the Supporting Information).

We repeated these analyses with the 400-node Schaefer parcellation. Although the pattern of results regarding the spatial scales affected was not identical, the results also indicated that turbulent measures were systematically decreased in TBI patients as compared with HCs (see Tables S6S9 in the Supporting Information).

In an additional analysis, we assessed longitudinal differences within TBI patients and the potential confounding effect of lesion volume by using a repeated measures ANCOVA with time points as within-group factor and lesion volume as a nuisance covariate. We found a significant Group × Lambda interaction in global amplitude turbulence (F(18,180) = 2.073, p = 0.008; Figure S1a in the Supporting Information). Post hoc comparisons revealed that global amplitude turbulence was significantly decreased after 6 months compared with 3 months post-injury for λ = 0.01 (T(22) = 3.844, p = 0.018; Supporting Information Table S10). Information cascade flow also showed a significant Group × Lambda interaction (F(16,160) = 1.809, p = 0.034; Supporting Information Figure S1b). Post hoc, we found a significant reduction at the 6-month post-injury stage compared with 3 months for λ = 0.01 (T(22) = 4.733, p < 0.001; Supporting Information Table S11). Information cascade was significantly reduced at the 6-month post-injury stage compared with the 12-month stage (group effect: F(2,20) = 4.205, p = 0.003; T(22) = −2.906, p = 0.026; Figure S1c and Table S12 in the Supporting Information). No significant differences were found for information transfer in the within-group analysis (Supporting Information Table S13).

Reductions in RSN empirical turbulence at the 6-month stage show recovery at the 12-month stage.

Next, we examined differences in amplitude turbulence between HCs and TBI patients in the seven RSNs included in the Schaefer parcellation for the spatial scales where we found significant reductions in global amplitude turbulence between HCs and TBI patients at all stages of recovery (λ = 0.01 and λ = 0.03; see previous section). The results of the one-way ANOVA including amplitude turbulence in HCs (averaged across both sessions) and TBI patients at 3, 6, and 12 months post-injury revealed a significant group effect for λ = 0.03 in all seven RSNs except for the limbic network (Table S14, Figure 3A). Post hoc we found that this effect was in all cases attributable to the difference between HCs and TBI patients at 6 months post-injury. In other words, this implies that after a decrease in turbulent dynamics at 6 months most RSNs recovered healthy levels by 12 months at this spatial scale. No significant group effect was found for λ = 0.01.

Figure 3.

Amplitude turbulence in most resting-state networks (RSNs) shows a U-shape recovery trajectory with no significant differences between healthy controls (HCs) and traumatic brain injury (TBI) patients by 12 months post-injury (λ = 0.03). Nodes in the somatosensory, dorsal attention, and control networks show the highest difference in TBI patients at 6 months post-injury. (A) Amplitude turbulence is shown for each dataset and the seven RSNs included in the Schaefer 1,000-node parcellation. Asterisks indicate significance level for the group effect in the comparison with HCs, averaged across sessions. (B) Kolmogorov-Smirnov distance (KSD) distributions for the difference in node-level turbulence between HCs (averaged across sessions) and TBI patients at 3 (left), 6 (middle), and 12 (right) months post-injury. Of note, the largest KSD values were reached for lower λ values (large distances) at the 6-month post-injury stage for all spatial scales (λ). (C) Rendered brains represent the absolute difference of the node-level metastability between each dataset for scale λ = 0.03. (D) Radar plots representing the number of nodes on the top 15% quantile of the absolute difference for the HCs vs. TBI patients at 3-, 6-, and 12-month post-injury stage comparison for λ = 0.03. At 6 months post-injury, most of the nodes were ascribed to the somatosensory, dorsal attention, and control networks. Abbreviations: VIS, visual; SM, somatosensory; DAT, dorsal attention; VAT, ventral attention; LIM, limbic; CNT, control; DMN, default mode. Asterisks: * p < 0.05, ** p < 0.01 (RSN turbulence); * denotes significant differences (KSD; for specific p values, see Supporting Information Table S14).

Figure 3.

Amplitude turbulence in most resting-state networks (RSNs) shows a U-shape recovery trajectory with no significant differences between healthy controls (HCs) and traumatic brain injury (TBI) patients by 12 months post-injury (λ = 0.03). Nodes in the somatosensory, dorsal attention, and control networks show the highest difference in TBI patients at 6 months post-injury. (A) Amplitude turbulence is shown for each dataset and the seven RSNs included in the Schaefer 1,000-node parcellation. Asterisks indicate significance level for the group effect in the comparison with HCs, averaged across sessions. (B) Kolmogorov-Smirnov distance (KSD) distributions for the difference in node-level turbulence between HCs (averaged across sessions) and TBI patients at 3 (left), 6 (middle), and 12 (right) months post-injury. Of note, the largest KSD values were reached for lower λ values (large distances) at the 6-month post-injury stage for all spatial scales (λ). (C) Rendered brains represent the absolute difference of the node-level metastability between each dataset for scale λ = 0.03. (D) Radar plots representing the number of nodes on the top 15% quantile of the absolute difference for the HCs vs. TBI patients at 3-, 6-, and 12-month post-injury stage comparison for λ = 0.03. At 6 months post-injury, most of the nodes were ascribed to the somatosensory, dorsal attention, and control networks. Abbreviations: VIS, visual; SM, somatosensory; DAT, dorsal attention; VAT, ventral attention; LIM, limbic; CNT, control; DMN, default mode. Asterisks: * p < 0.05, ** p < 0.01 (RSN turbulence); * denotes significant differences (KSD; for specific p values, see Supporting Information Table S14).

Close modal

To complement this analysis, we computed the Kolmogorov-Smirnov distance (KSD) that quantifies the difference in the distributions of node-level metastability between HCs and TBI patients. We found that the KSD monotonically decreases (i.e., distributions are more similar) across scales in all group comparisons, whereas the value of λ increases. In other words, the KSD is maximal for lower values of λ, that is, long distances in the brain (Figure 3B). Remarkably, the highest KSD is found between HCs and TBI patients at 6 months post-injury (Supporting Information Table S15).

For illustrative purposes, Figure 3C displays the absolute difference in node-level metastability between HCs and TBI patients at 3, 6, and 12 months post-injury at λ = 0.03 rendered onto the brain cortex. Additionally, we computed the number of nodes on the top 15% quantile of the absolute difference in node-level metastability for the HCs versus TBI patients at the 6-month post-injury comparison (Figure 3D). The nodes with the largest differences in node-level metastability were located in the dorsal attention, somatomotor, and control networks.

Empirical measures of turbulent dynamics predict cognitive performance.

To determine whether empirical measures of turbulent dynamics relate to cognitive impairment in TBI patients, we performed Pearson’s correlations with the neuropsychological battery of tests and turbulent measures at 3 months post-injury as provided by the authors of the original publication. The battery included those tests more sensitive to detect deficits in processing speed and working memory that are common after TBI (see the Methods section). We found a positive association between Digit Span (Forward and Backward) and information cascade flow at λ = 0.01 (large spatial scale) (r = 0.754, p = 0.007; Figure 4A, Supporting Information Table S16). This result survived multiple-comparison correction with FDR across a number of battery tests (p-adj. = 0.042). In other words, the greater the propagation of information between consecutive spatial scales, the better the working memory performance. Furthermore, we found that performance in this working memory task was also positively associated with turbulence in the default mode network for λ = 0.03 (r = 0.744, p = 0.009; Figure 4B, Supporting Information Table S17). This result also survived multiple-comparison correction with FDR across a number of battery tests (p-adj. = 0.045). Lesion volume did not correlate with either cognitive performance or amplitude turbulence.

Figure 4.

Working memory performance is positively associated with turbulent measures in traumatic brain injury (TBI) patients at 3 months post-injury. Scatterplots representing the positive correlation between z scores in the WAIS-III Digit Span (Forward and Backward) and (A) information cascade flow for λ = 0.01 and (B) turbulence in the DMN for λ = 0.03. Abbreviations: F: forward; B: backward; DMN: default mode network.

Figure 4.

Working memory performance is positively associated with turbulent measures in traumatic brain injury (TBI) patients at 3 months post-injury. Scatterplots representing the positive correlation between z scores in the WAIS-III Digit Span (Forward and Backward) and (A) information cascade flow for λ = 0.01 and (B) turbulence in the DMN for λ = 0.03. Abbreviations: F: forward; B: backward; DMN: default mode network.

Close modal

Model-Based Framework

We defined a whole-brain model of coupled oscillators (Stuart-Landau oscillators) for each time point of the TBI patients’ and HCs’ rsfMRI data. To account for the effect of focal lesions on the SC, we adopted a simulated attack approach with lesion mask matrices thresholded at 1.5 and 2 SD and applied a binary or weighted attack, resulting in a total of four simulated attacks (see Figure 1B). Additionally, for each time point of the TBI patients’ dataset, a computational model was run with intact healthy SC to assess the sensitivity of the model to the simulated attack. In total, 11 different whole-brain computational models were computed for the 1.5 SD lesion mask. To optimize the fitting between the whole-brain model and the empirical rsfMRI data, we varied the global coupling parameter G from 0 to 3 in steps of 0.01, and for each G value, we ran 100 simulations with the same TR (2 s) and time duration (145 scans) as the empirical rsfMRI data. We then determined the optimal working point of each model as the minimum of the fitting level. Next, we calculated differences in susceptibility and information-encoding capability between groups (HCs Session 1 and 2; TBI patients at 3, 6, and 12 months post-injury) after in silico perturbations.

We first identified the working point G for each rsfMRI dataset, which represents the global coupling parameter. First, we found that the working point of the model was lower for the TBI patients at all time points (3, 6, and 12 months post-injury) as compared with the value for both sessions of HCs (Figures 5AC), regardless of whether an intact or lesioned SC was used as input to the model. Crucially, we detected the same pattern, that is, the global coupling parameter for TBI patients at 6 months was smaller than that at 3 and 12 months post-injury. Second, the global coupling parameter shifted to lower values after the simulated attack with both thresholds (1.5 and 2 SD). However, the difference with respect to the intact healthy SC was small. And the binary approach did not yield smaller values than the weighted approach (see Table S18 in the Supporting Information).

Figure 5.

Traumatic brain injury (TBI) patients show reduced reactivity to in silico perturbations over time as compared with healthy controls (HCs), with the lowest value at 6 months post-injury with intact structural connectivity (SC) and after the simulated attack. (A–C) Evolution of the error of fitting the whole-brain model to the rsfMRI data as a function of the global coupling strength, G. The optimal working point fitting the TBI data was lower than that for HCs in all conditions: with intact SC or after a simulated attack with binary and weighted lesion masks thresholded at 1.5 SD from the group’s mean volume of overlap with the Schaefer 1,000-node parcellation (see the Methods section). The shadow indicates the standard error of the mean. (D) Information-encoding capability and (E) susceptibility measures were significantly different between HCs and TBI patients across all three time points, with the lowest value shown by TBI patients at 6 months post-injury in all conditions. Further, the values with the binary approach for the simulated attack were also lower than with the weighted approach. The lower and upper edges of the box indicate the first (25th percentile) and third quartile (75th percentile), and middle lines are medians. Orange asterisks denote comparisons between HCs Session 1 and 2 and TBI patients (*** p < 0.0001, ** p < 0.01, results after Wilcoxon rank sum test comparing the distribution of values in 100 trials per condition). Black asterisks denote comparisons between intact SC and simulated attack approaches (*** p < 0.001, results after Wilcoxon signed rank test comparing the distribution of values in 100 trials per condition). Abbreviations: mo: months; les: lesion; 1.5: 1.5 SD threshold; B: binary; W: weighted.

Figure 5.

Traumatic brain injury (TBI) patients show reduced reactivity to in silico perturbations over time as compared with healthy controls (HCs), with the lowest value at 6 months post-injury with intact structural connectivity (SC) and after the simulated attack. (A–C) Evolution of the error of fitting the whole-brain model to the rsfMRI data as a function of the global coupling strength, G. The optimal working point fitting the TBI data was lower than that for HCs in all conditions: with intact SC or after a simulated attack with binary and weighted lesion masks thresholded at 1.5 SD from the group’s mean volume of overlap with the Schaefer 1,000-node parcellation (see the Methods section). The shadow indicates the standard error of the mean. (D) Information-encoding capability and (E) susceptibility measures were significantly different between HCs and TBI patients across all three time points, with the lowest value shown by TBI patients at 6 months post-injury in all conditions. Further, the values with the binary approach for the simulated attack were also lower than with the weighted approach. The lower and upper edges of the box indicate the first (25th percentile) and third quartile (75th percentile), and middle lines are medians. Orange asterisks denote comparisons between HCs Session 1 and 2 and TBI patients (*** p < 0.0001, ** p < 0.01, results after Wilcoxon rank sum test comparing the distribution of values in 100 trials per condition). Black asterisks denote comparisons between intact SC and simulated attack approaches (*** p < 0.001, results after Wilcoxon signed rank test comparing the distribution of values in 100 trials per condition). Abbreviations: mo: months; les: lesion; 1.5: 1.5 SD threshold; B: binary; W: weighted.

Close modal

Once we identified the working point G for each rsfMRI dataset, we introduced external perturbations by changing the local bifurcation parameter an of each brain node n. For each value of G, we perturbed the whole-brain model 100 times with random parameters for the local bifurcation parameter. Our model-based measures, susceptibility and information-encoding capability, allowed us to evaluate the reactivity to the external perturbations introduced (for more details, see the Methods section). Briefly, susceptibility was estimated by first measuring the perturbed and nonperturbed modulus of the Kuramoto local order parameter and then calculating the difference between the perturbed and nonperturbed cases averaged across time, trials, and brain nodes. As an extension of this measure, the information-encoding capability of the whole-brain model was calculated as the standard deviation across trials of the difference between the perturbed and unperturbed mean of the modulus of the local order parameter across time, averaged over all brain nodes.

Figures 5D5E show that both susceptibility and information-encoding capability were significantly reduced in TBI patients across the recovery trajectory compared with both sessions of HCs regardless of whether an intact or lesioned SC was used as input to the model. However, the values of these measures for the binary simulated attack were significantly reduced at all stages of recovery with respect to the values in the weighted simulated attack, suggesting that the model is most sensitive to focal lesions that greatly overlap. Since in this cohort of TBI patients the frequency of lesion per node was mostly 1 (Supporting Information Table S19), the results from the weighted attack may be more representative of the underlying average SC. With the 1.5 SD weighted lesion mask matrix, TBI patients at 6 months post-injury showed the lowest mean value across trials (Figures 5D5E, Supporting Information Tables S20S21). This pattern of results was replicated when using the weighted 2 SD lesion mask (Figure S2 and Tables S22S23 in the Supporting Information). Therefore, the results from our in silico stimulation suggest that, at 6 months post-injury, the brain in this cohort of TBI patients is working in a regime that is the least responsive to external perturbations.

Here we investigated how TBI impacts whole-brain resting-state fMRI dynamics using a turbulence-based framework (Deco & Kringelbach, 2020; Deco et al., 2021). Our findings revealed specific spatial patterns of reduction in global amplitude turbulence between HCs and TBI patients at all three time points (3, 6, and 12 months post-injury), differentially affecting long distances in the brain. Noteworthily, the greatest mean difference was observed between HCs and TBI patients at 6 months at all significant λ values. This was also the stage of recovery where we found differences in amplitude turbulence over time at long distances. Between-group differences were also found at long distances for the other turbulent measures except for information transfer, where short distances were affected and no differences were found in the longitudinal analysis. Importantly, when analyzing network-level turbulence between HCs and TBI patients for λ = 0.03 (long distances), we found that, after a significant decrease at 6 months post-injury, amplitude turbulence recovers healthy levels in all RSNs except for the limbic by 12 months. Finally, the results from the whole-brain model indicated that susceptibility and information-encoding capability are significantly reduced after TBI, with the largest difference again in the comparison between HCs and TBI patients at 6 months post-injury regardless of whether an intact or lesioned SC was used as input to the model. These results provide causal links needed for a better understanding of the recovery process 1 year after TBI. To the best of our knowledge, this is the first study to examine the longitudinal evolution of turbulent dynamics during recovery from TBI.

Our empirical finding of reduced turbulence at long distances after TBI as compared with HCs is in alignment with previous work showing a reduction in global metastability in TBI patients as compared with HCs (Hellyer et al., 2015). Furthermore, the turbulent dynamics framework adopted here lends indirect support to the graph theory results obtained in the same dataset (Roy et al., 2017) and extends those findings by offering a dynamical approach across spatial scales and a causal perspective on the information-processing deficits following TBI during 1 year. Previously, Roy et al. (2017) had shown that global functional hyperconnectivity, in terms of network strength, peaked at 6 months post-injury. Of note, the local analysis revealed hyperconnectivity in the left frontal DMN and temporoparietal attentional control networks across time points, but it was only associated with an increased cost at 6 months post-injury. Moreover, a cost-efficiency analysis showed that this elevated cost was mostly due to a higher number of medium-range connections rather than long-range connections. This fits well with our network-level finding, where turbulence is reduced at 6 months post-injury in most RSNs with no residual effects after 12 months, when there was no cost associated with hyperconnectivity. Remarkably, we also found that the 15% quantile nodes in the HCs versus TBI patients at 6 months were in the somatosensory, dorsal attention, and control networks, in accordance with increased cost found in the temporoparietal control network by Roy et al. (2017). Overall, our network-level findings suggest a shift in the functional plasticity of the brain in response to TBI, where the turbulent dynamics at large distances rebalance to a healthy regime during the period spanning from 6 months to 12 months post-injury.

Regarding brain-behavior relationships, we found that performance on the Digit Span task (Forward and Backward) was positively associated with information cascade flow and amplitude turbulence in the DMN at long distances. This task is used to evaluate working memory processes, which are commonly affected after TBI, and requires the rapid communication across large-scale brain networks involving sensory (patients are read a sequence of numbers) and higher order regions responsible for the manipulation of information for temporary use. Our findings suggest that the propagation of information between consecutive spatial scales and higher levels of amplitude in the DMN at rest might be an important dynamical mechanism to support working memory tasks. However, future work with larger sample sizes could better characterize additional associations between turbulent measures and performance in other cognitive domains, as this has been previously found for metastability and cognitive flexibility, information processing, and associative memory after TBI (Hellyer et al., 2015).

Our computational findings are consistent with our empirical results, demonstrating that TBI hinders the reactivity of the brain to external in silico perturbations. When fitting the Hopf model to each empirical condition, we found lower values in the optimal working point for TBI patients across time points regardless of whether an intact or lesioned SC was used as input to the model. However, significant differences between the binary and the weighted approach emerged when examining the susceptibility and information-encoding capability after our in silico perturbation. Indeed, the values of these measures for the binary simulated attack were significantly reduced at all stages of recovery with respect to the values in the weighted simulated attack. This is indicative that these measures are less affected by a distributed pattern of lesions, as the one in this cohort of TBI patients, and may be more sensitive to focal lesions that greatly overlap in the sample. Whether this result is contingent upon the anatomical identity (and connectivity profile) of the lesions remains an interesting question for future research. That said, we believe that this is an interesting finding on its own, highlighting the importance of the impact that focal lesions may have in other patient populations, such as stroke, where usually most patients present a lesion in the vascular territory of the middle cerebral artery and there is greater degree of overlap.

The results from the present study have some potential clinical implications. Turbulent measures might be used to improve clinical diagnosis and prognosis after TBI, as they demonstrated sensitivity to the longitudinal changes in whole-brain dynamics compared with HCs. However, further information regarding the relationship of these measures and TBI severity is still needed and calls for collaborative efforts to perform these analyses on large cohorts of patients. On the other hand, information-processing-based measures obtained from resting-state data could provide a meaningful and repeatable way to monitor the evolution of TBI, the recovery progress, and the patient’s possible response to treatments, as these remain prominent challenges in clinical settings to date (National Academies of Sciences, Engineering, and Medicine, 2022). Our in silico perturbation approach with the Hopf computational model also indicated reduced coupling and reactivity to external stimulation after TBI across time points. We restricted our perturbative protocol to the introduction of random values in all nodes to modify the local bifurcation parameter of the model. That being said, it would be interesting to implement an excitatory protocol for individual nodes or a subset of relevant nodes in order to identify brain targets that can be stimulated with techniques such as transcranial magnetic stimulation (TMS) (Demirtas-Tatlidede et al., 2012). Notably, some promising results indicate that TMS and neurorehabilitation may act synergistically to improve functional outcomes after TBI (Nardone et al., 2020). Yet, we leave the implementation of such perturbative protocol for future work.

It is also worthwhile to acknowledge some limitations in the current study. Although this dataset represents a relatively homogenous sample in terms of age and time since injury, it remains of small size and a replication with a larger cohort of patients is needed to assess the generalizability of these results. In addition, the resting-state sequence comprised 150 scans with a temporal resolution of 2 s that were acquired over a 5-min period. While it has been shown that functional correlation of RSNs only stabilizes after 5 min with TR = 2.3 s (Van Dijk et al., 2010), recent evidence indicates that metastability is less affected by temporal resolution and shows greater reliability in short-term scans (Yang et al., 2022). Even if our simulated attack approach allowed us to take into account the effect of focal lesions on the SC, it leaves aside interindividual differences in SC, and future studies with personalized models are warranted to validate our findings. Finally, although we were able to partially replicate our results using the 400-node Schaefer parcellation, it remains unknown how they may be affected when using a different atlas.

In conclusion, we have shown that whole-brain turbulent dynamics complement previous findings on the temporal evolution of functional connectivity after TBI in an open-access dataset of healthy controls and patients with data collected after 3, 6, and 12 months post-injury. Specifically, our empirical results suggest the presence of maladaptive neuroplasticity at 6 months post-injury as manifested in the reduction of global amplitude turbulence at long distances. Thus, we demonstrated that the framework of turbulent dynamics may facilitate our understanding of how the brain responds to traumatic injury. Importantly, our computational approach suggests that in this cohort of moderate to severe TBI patients, their brain is the least reactive to external perturbations at 6 months post-injury during the recovery trajectory.

We would like to thank Professor Frank Hillary and Andrew Cwiek for kindly sharing the lesion maps and neuropsychological assessment from the TBI patients included in this OpenNeuro dataset.

All code written in support of this study is publicly available on https://github.com/noechan/TBI_ON_turbu_Hopf_NC (Martínez-Molina, 2023).

Supporting information for this article is available at https://doi.org/10.1162/netn_a_00346.

Noelia Martínez-Molina: Conceptualization; Data curation; Formal analysis; Funding acquisition; Project administration; Visualization; Writing – original draft. Anira Escrichs: Formal analysis; Software; Writing – review & editing. Yonatan Sanz-Perl: Software; Validation; Writing – review & editing. Aleksi J. Sihvonen: Writing – review & editing. Teppo Särkämö: Writing – review & editing. Morten L. Kringelbach: Methodology; Software; Visualization; Writing – review & editing. Gustavo Deco: Methodology; Software; Supervision; Writing – review & editing.

Noelia Martínez-Molina, BdP Programme, Award ID: 2019-BP-00032. Yonatan Sanz-Perl, H2020 Marie Skłodowska-Curie Actions (https://dx.doi.org/10.13039/100010665), Award ID: 896354. Anira Escrichs, Human Brain Mapping Project, Award ID: 945539. Aleksi J. Sihvonen, Suomen Kulttuurirahasto (https://dx.doi.org/10.13039/501100003125), Award ID: 191230. Aleksi J. Sihvonen, Orionin Tutkimussäätiö (https://dx.doi.org/10.13039/501100007083). Aleksi J. Sihvonen, Signe ja Ane Gyllenbergin Säätiö (https://dx.doi.org/10.13039/501100004325). Teppo Särkämö, Academy of Finland (https://dx.doi.org/10.13039/501100002341), Award ID: 338448 & 346211. Teppo Särkämö, H2020 European Research Council (https://dx.doi.org/10.13039/100010663), Award ID: 803466. Morten L. Kringelbach, Danish National Research Foundation, Award ID: DNRF117. Morten L. Kringelbach is the founder of the Centre for Eudaimonia and Human Flourishing at Linacre College, funded by the Pettit and Carlsberg Foundations. Gustavo Deco, Ministerio de Ciencia e Innovación (https://dx.doi.org/10.13039/501100004837), Award ID: PID2019-105772GB-I00 MCIU AEI.

Turbulence:

In coupled oscillators, turbulence is characterized by the high variability across space and time of the Kuramoto local order parameter capturing the local synchronization, so in this sense it is a spatiotemporal extension of the level of metastability.

Small world:

A class of networks that show high clustering, much like a lattice, but with a short characteristic path length, much like a random graph.

Diffuse axonal injury:

Damage caused to axons and axoplasmic membranes by shearing forces produced at the time of brain injury. This impairs axonal transport and disrupts brain function and structure.

Metastability:

Fundamental concept describing complex systems’ behavior in terms of the variability in the global state of synchronization as a function of time.

Tractography:

A class of methods typically applied to diffusion MRI data to trace the trajectories of axonal fiber bundles based on preferred directions of water diffusion in the brain.

Information cascade flow:

A measure that estimates the stream of information between a given scale and a subsequent lower scale in consecutive time steps.

Information cascade:

A measure that captures the entire efficient information-processing behavior across scales.

Information transfer:

A measure that characterizes how the information propagates across space at a particular scale.

Bassett
,
D. S.
, &
Bullmore
,
E.
(
2006
).
Small-world brain networks
.
Neuroscientist
,
12
(
6
),
512
523
. ,
[PubMed]
Benjamini
,
Y.
, &
Hochberg
,
Y.
(
1995
).
Controlling the false discovery rate: A practical and powerful approach to multiple testing
.
Journal of the Royal Statistical Society: Series B (Methodological)
,
57
(
1
),
289
300
.
Bonnelle
,
V.
,
Leech
,
R.
,
Kinnunen
,
K. M.
,
Ham
,
T. E.
,
Beckmann
,
C. F.
,
De Boissezon
,
X.
,
Greenwood
,
R. J.
, &
Sharp
,
D. J.
(
2011
).
Default mode network connectivity predicts sustained attention deficits after traumatic brain injury
.
Journal of Neuroscience
,
31
(
38
),
13442
13451
. ,
[PubMed]
Breakspear
,
M.
(
2017
).
Dynamic models of large-scale brain activity
.
Nature Neuroscience
,
20
(
3
),
340
352
. ,
[PubMed]
Cabral
,
J.
,
Luckhoo
,
H.
,
Woolrich
,
M.
,
Joensson
,
M.
,
Mohseni
,
H.
,
Baker
,
A.
,
Kringelbach
,
M. L.
, &
Deco
,
G.
(
2014
).
Exploring mechanisms of spontaneous functional connectivity in MEG: How delayed network interactions lead to structured amplitude envelopes of band-pass filtered oscillations
.
NeuroImage
,
90
,
423
435
. ,
[PubMed]
Caeyenberghs
,
K.
,
Leemans
,
A.
,
Leunissen
,
I.
,
Gooijers
,
J.
,
Michiels
,
K.
,
Sunaert
,
S.
, &
Swinnen
,
S. P.
(
2014
).
Altered structural networks and executive deficits in traumatic brain injury patients
.
Brain Structure and Function
,
219
(
1
),
193
209
. ,
[PubMed]
Cruzat
,
J.
,
Perl
,
Y. S.
,
Escrichs
,
A.
,
Vohryzek
,
J.
,
Timmermann
,
C.
,
Roseman
,
L.
,
Luppi
,
A. I.
,
Ibañez
,
A.
,
Nutt
,
D.
,
Carhart-Harris
,
R.
,
Tagliazucchi
,
E.
,
Deco
,
G.
, &
Kringelbach
,
M. L.
(
2022
).
Effects of classic psychedelic drugs on turbulent signatures in brain dynamics
.
Network Neuroscience
,
6
(
4
),
1104
1124
.
Deco
,
G.
,
Cabral
,
J.
,
Saenger
,
V. M.
,
Boly
,
M.
,
Tagliazucchi
,
E.
,
Laufs
,
H.
,
Van Someren
,
E.
,
Jobst
,
B.
,
Stevner
,
A.
, &
Kringelbach
,
M. L.
(
2018
).
Perturbation of whole-brain dynamics in silico reveals mechanistic differences between brain states
.
NeuroImage
,
169
,
46
56
. ,
[PubMed]
Deco
,
G.
,
Cruzat
,
J.
,
Cabral
,
J.
,
Tagliazucchi
,
E.
,
Laufs
,
H.
,
Logothetis
,
N. K.
, &
Kringelbach
,
M. L.
(
2019
).
Awakening: Predicting external stimulation to force transitions between different brain states
.
Proceedings of the National Academy of Sciences
,
116
(
36
),
18088
18097
. ,
[PubMed]
Deco
,
G.
, &
Kringelbach
,
M. L.
(
2014
).
Great expectations: Using whole-brain computational connectomics for understanding neuropsychiatric disorders
.
Neuron
,
84
(
5
),
892
905
. ,
[PubMed]
Deco
,
G.
, &
Kringelbach
,
M. L.
(
2020
).
Turbulent-like dynamics in the human brain
.
Cell Reports
,
33
(
10
),
108471
. ,
[PubMed]
Deco
,
G.
,
Kringelbach
,
M. L.
,
Jirsa
,
V. K.
, &
Ritter
,
P.
(
2017
).
The dynamics of resting fluctuations in the brain: Metastability and its dynamical cortical core
.
Scientific Reports
,
7
(
1
),
3095
. ,
[PubMed]
Deco
,
G.
,
Sanz Perl
,
Y.
,
Vuust
,
P.
,
Tagliazucchi
,
E.
,
Kennedy
,
H.
, &
Kringelbach
,
M. L.
(
2021
).
Rare long-range cortical connections enhance human information processing
.
Current Biology
,
31
(
20
),
4436
4448
. ,
[PubMed]
Deco
,
G.
,
Tagliazucchi
,
E.
,
Laufs
,
H.
,
Sanjuán
,
A.
, &
Kringelbach
,
M. L.
(
2017
).
Novel intrinsic ignition method measuring local-global integration characterizes wakefulness and deep sleep
.
eNeuro
,
4
(
5
),
ENEURO.0106-17.2017
. ,
[PubMed]
Deco
,
G.
,
Tononi
,
G.
,
Boly
,
M.
, &
Kringelbach
,
M. L.
(
2015
).
Rethinking segregation and integration: Contributions of whole-brain modelling
.
Nature Reviews Neuroscience
,
16
(
7
),
430
439
. ,
[PubMed]
De Filippi
,
E.
,
Uribe
,
C.
,
Avila-Varela
,
D. S.
,
Martínez-Molina
,
N.
,
Gashaj
,
V.
,
Pritschet
,
L.
,
Santander
,
T.
,
Jacobs
,
E. G.
,
Kringelbach
,
M. L.
,
Sanz Perl
,
Y.
,
Deco
,
G.
, &
Escrichs
,
A.
(
2021
).
The menstrual cycle modulates whole-brain turbulent dynamics
.
Frontiers in Neuroscience
,
15
,
753820
. ,
[PubMed]
Demertzi
,
A.
,
Antonopoulos
,
G.
,
Heine
,
L.
,
Voss
,
H. U.
,
Crone
,
J. S.
,
de Los Angeles
,
C.
,
Bahri
,
M. A.
,
Di Perri
,
C.
,
Vanhaudenhuyse
,
A.
,
Charland-Verville
,
V.
,
Kronbichler
,
M.
,
Trinka
,
E.
,
Phillips
,
C.
,
Gomez
,
F.
,
Tshibanda
,
L.
,
Soddu
,
A.
,
Schiff
,
N. D.
,
Whitfield-Gabrieli
,
S.
, &
Laureys
,
S.
(
2015
).
Intrinsic functional connectivity differentiates minimally conscious from unresponsive patients
.
Brain
,
138
(
9
),
2619
2631
. ,
[PubMed]
Demirtas-Tatlidede
,
A.
,
Vahabzadeh-Hagh
,
A. M.
,
Bernabeu
,
M.
,
Tormos
,
J. M.
, &
Pascual-Leone
,
A.
(
2012
).
Noninvasive brain stimulation in traumatic brain injury
.
Journal of Head Trauma Rehabilitation
,
27
(
4
),
274
292
. ,
[PubMed]
Escrichs
,
A.
,
Perl
,
Y. S.
,
Uribe
,
C.
,
Camara
,
E.
,
Türker
,
B.
,
Pyatigorskaya
,
N.
,
López-González
,
A.
,
Pallavicini
,
C.
,
Panda
,
R.
,
Annen
,
J.
,
Gosseries
,
O.
,
Laureys
,
S.
,
Naccache
,
L.
,
Sitt
,
J. D.
,
Laufs
,
H.
,
Tagliazucchi
,
E.
,
Kringelbach
,
M. L.
, &
Deco
,
G.
(
2022
).
Unifying turbulent dynamics framework distinguishes different brain states
.
Communications Biology
,
5
(
1
),
638
. ,
[PubMed]
Fox
,
M. D.
,
Snyder
,
A. Z.
,
Vincent
,
J. L.
,
Corbetta
,
M.
,
Van Essen
,
D. C.
, &
Raichle
,
M. E.
(
2005
).
The human brain is intrinsically organized into dynamic, anticorrelated functional networks
.
Proceedings of the National Academy of Sciences
,
102
(
27
),
9673
9678
. ,
[PubMed]
Friston
,
K. J.
(
1997
).
Transients, metastability, and neuronal dynamics
.
NeuroImage
,
5
(
2
),
164
171
. ,
[PubMed]
Glerean
,
E.
,
Salmi
,
J.
,
Lahnakoski
,
J. M.
,
Jaaskelainen
,
I. P.
, &
Sams
,
M.
(
2012
).
Functional magnetic resonance imaging phase synchronization as a measure of dynamic functional connectivity
.
Brain Connectivity
,
2
(
2
),
91
101
. ,
[PubMed]
Hellyer
,
P. J.
,
Leech
,
R.
,
Ham
,
T. E.
,
Bonnelle
,
V.
, &
Sharp
,
D. J.
(
2013
).
Individual prediction of white matter injury following traumatic brain injury
.
Annals of Neurology
,
73
(
4
),
489
499
. ,
[PubMed]
Hellyer
,
P. J.
,
Scott
,
G.
,
Shanahan
,
M.
,
Sharp
,
D. J.
, &
Leech
,
R.
(
2015
).
Cognitive flexibility through metastable neural dynamics is disrupted by damage to the structural connectome
.
Journal of Neuroscience
,
35
(
24
),
9050
9063
. ,
[PubMed]
Horn
,
A.
, &
Blankenburg
,
F.
(
2016
).
Toward a standardized structural-functional group connectome in MNI space
.
NeuroImage
,
124
(
Pt. A
),
310
322
. ,
[PubMed]
Horn
,
A.
,
Neumann
,
W. J.
,
Degen
,
K.
,
Schneider
,
G. H.
, &
Kühn
,
A. A.
(
2017
).
Toward an electrophysiological “sweet spot” for deep brain stimulation in the subthalamic nucleus
.
Human Brain Mapping
,
38
(
7
),
3377
3390
. ,
[PubMed]
Jilka
,
S. R.
,
Scott
,
G.
,
Ham
,
T.
,
Pickering
,
A.
,
Bonnelle
,
V.
,
Braga
,
R. M.
,
Leech
,
R.
, &
Sharp
,
D. J.
(
2014
).
Damage to the salience network and interactions with the default mode network
.
Journal of Neuroscience
,
34
(
33
),
10798
10807
. ,
[PubMed]
Jolly
,
A. E.
,
Bălăeţ
,
M.
,
Azor
,
A.
,
Friedland
,
D.
,
Sandrone
,
S.
,
Graham
,
N. S. N.
,
Zimmerman
,
K.
, &
Sharp
,
D. J.
(
2021
).
Detecting axonal injury in individual patients after traumatic brain injury
.
Brain
,
144
(
1
),
92
113
. ,
[PubMed]
Kawamura
,
Y.
,
Nakao
,
H.
, &
Kuramoto
,
Y.
(
2007
).
Noise-induced turbulence in nonlocally coupled oscillators
.
Physical Review E: Statistical, Nonlinear, and Soft Matter Physics
,
75
(
3 Pt. 2
),
036209
. ,
[PubMed]
Kinnunen
,
K. M.
,
Greenwood
,
R.
,
Powell
,
J. H.
,
Leech
,
R.
,
Hawkins
,
P. C.
,
Bonnelle
,
V.
,
Patel
,
M. C.
,
Counsell
,
S. J.
, &
Sharp
,
D. J.
(
2011
).
White matter damage and cognitive impairment after traumatic brain injury
.
Brain
,
134
(
2
),
449
463
. ,
[PubMed]
Kringelbach
,
M. L.
, &
Deco
,
G.
(
2020
).
Brain states and transitions: Insights from computational neuroscience
.
Cell Reports
,
32
(
10
),
108128
. ,
[PubMed]
Kolmogorov
,
A. N.
(
1941
).
Dissipation of energy in the locally isotropic turbulence
.
Doklady Akademii Nauk SSSR A
,
32
,
16
18
.
Kuramoto
,
Y.
(
1984
).
Chemical turbulence
. In
Chemical oscillations, waves, and turbulence
(pp.
111
140
).
Berlin
:
Springer
.
López-González
,
A.
,
Panda
,
R.
,
Ponce-Alvarez
,
A.
,
Zamora-López
,
G.
,
Escrichs
,
A.
,
Martial
,
C.
,
Thibaut
,
A.
,
Gosseries
,
O.
,
Kringelbach
,
M. L.
,
Annen
,
J.
,
Laureys
,
S.
, &
Deco
,
G.
(
2021
).
Loss of consciousness reduces the stability of brain hubs and the heterogeneity of brain dynamics
.
Communications Biology
,
4
(
1
),
1037
. ,
[PubMed]
Martínez-Molina
,
N.
(
2023
).
Whole-brain turbulent-like dynamics and Hopf computational model, GitHub
, https://github.com/noechan/TBI_ON_turbu_Hopf
Martínez-Molina
,
N.
,
Siponkoski
,
S.-T.
,
Kuusela
,
L.
,
Laitinen
,
S.
,
Holma
,
M.
,
Ahlfors
,
M.
,
Jordan-Kilkki
,
P.
,
Ala-Kauhaluoma
,
K.
,
Melkas
,
S.
,
Pekkola
,
J.
,
Rodriguez-Fornells
,
A.
,
Laine
,
M.
,
Ylinen
,
A.
,
Rantanen
,
P.
,
Koskinen
,
S.
,
Cowley
,
B. U.
, &
Särkämo
,
T.
(
2021
).
Resting-state network plasticity induced by music therapy after traumatic brain injury
.
Neural Plasticity
,
2021
,
6682471
. ,
[PubMed]
McAllister
,
T. W.
,
Flashman
,
L. A.
,
Sparling
,
M. B.
, &
Saykin
,
A. J.
(
2004
).
Working memory deficits after traumatic brain injury: Catecholaminergic mechanisms and prospects for treatment—A review
.
Brain Injury
,
18
(
4
),
331
350
. ,
[PubMed]
Medaglia
,
J. D.
,
Erickson
,
B. A.
,
Pustina
,
D.
,
Kelkar
,
A. S.
,
DeMarco
,
A. T.
,
Dickens
,
J. V.
, &
Turkeltaub
,
P. E.
(
2022
).
Simulated attack reveals how lesions affect network properties in poststroke aphasia
.
Journal of Neuroscience
,
42
(
24
),
4913
4926
. ,
[PubMed]
Momi
,
D.
,
Wang
,
Z.
, &
Griffiths
,
J. D.
(
2023
).
TMS-evoked responses are driven by recurrent large-scale network dynamics
.
eLife
,
12
,
e83232
. ,
[PubMed]
Nardone
,
R.
,
Sebastianelli
,
L.
,
Versace
,
V.
,
Brigo
,
F.
,
Golaszewski
,
S.
,
Manganotti
,
P.
,
Saltuari
,
L.
, &
Trinka
,
E.
(
2020
).
Repetitive transcranial magnetic stimulation in traumatic brain injury: Evidence from animal and human studies
.
Brain Research Bulletin
,
159
,
44
52
. ,
[PubMed]
National Academies of Sciences, Engineering, and Medicine
. (
2022
).
Traumatic brain injury: A roadmap for accelerating progress
.
Washington, DC
:
National Academies Press
.
Nieto-Castanon
,
A.
(
2020
).
Handbook of functional connectivity magnetic resonance imaging methods in CONN
.
Boston
:
Hilbert Press
.
Ponsford
,
J. L.
,
Downing
,
M. G.
,
Olver
,
J.
,
Ponsford
,
M.
,
Acher
,
R.
,
Carty
,
M.
, &
Spitz
,
G.
(
2014
).
Longitudinal follow-up of patients with traumatic brain injury: Outcome at two, five, and ten years post-injury
.
Journal of Neurotrauma
,
31
(
1
),
64
77
. ,
[PubMed]
Roy
,
A.
,
Bernier
,
R. A.
,
Wang
,
J.
,
Benson
,
M.
,
French
,
J. J.
, Jr.
,
Good
,
D. C.
, &
Hillary
,
F. G.
(
2017
).
The evolution of cost-efficiency in neural networks during recovery from traumatic brain injury
.
PLOS ONE
,
12
(
4
),
e0170541
. ,
[PubMed]
Schaefer
,
A.
,
Kong
,
R.
,
Gordon
,
E. M.
,
Laumann
,
T. O.
,
Zuo
,
X. N.
,
Holmes
,
A. J.
,
Eickhoff
,
S. B.
, &
Yeo
,
B. T. T.
(
2018
).
Local-global parcellation of the human cerebral cortex from intrinsic functional connectivity MRI
.
Cerebral Cortex
,
28
(
9
),
3095
3114
. ,
[PubMed]
Setsompop
,
K.
,
Kimmlingen
,
R.
,
Eberlein
,
E.
,
Witzel
,
T.
,
Cohen-Adad
,
J.
,
McNab
,
J. A.
,
Keil
,
B.
,
Tisdall
,
M. D.
,
Hoecht
,
P.
,
Dietz
,
P.
,
Cauley
,
S. F.
,
Tountcheva
,
V.
,
Matschl
,
V.
,
Lenz
,
V. H.
,
Heberlein
,
K.
,
Potthast
,
A.
,
Thein
,
H.
,
Van Horn
,
J.
,
Toga
,
A.
, …
Wald
,
L. L.
(
2013
).
Pushing the limits of in vivo diffusion MRI for the Human Connectome Project
.
NeuroImage
,
80
,
220
233
. ,
[PubMed]
Shanahan
,
M.
(
2010
).
Metastable chimera states in community-structured oscillator networks
.
Chaos
,
20
(
1
),
013108
. ,
[PubMed]
Sharp
,
D. J.
,
Scott
,
G.
, &
Leech
,
R.
(
2014
).
Network dysfunction after traumatic brain injury
.
Nature Reviews Neurology
,
10
(
3
),
156
166
. ,
[PubMed]
Tognoli
,
E.
, &
Kelso
,
J. A.
(
2014
).
The metastable brain
.
Neuron
,
81
(
1
),
35
48
. ,
[PubMed]
Van Dijk
,
K. R. A.
,
Hedden
,
T.
,
Venkataraman
,
A.
,
Evans
,
K. C.
,
Lazar
,
S. W.
, &
Buckner
,
R. L.
(
2010
).
Intrinsic functional connectivity as a tool for human connectomics: Theory, properties, and optimization
.
Journal of Neurophysiology
,
103
(
1
),
297
321
. ,
[PubMed]
Yang
,
L.
,
Wei
,
J.
,
Li
,
Y.
,
Wang
,
B.
,
Guo
,
H.
,
Yang
,
Y.
, &
Xiang
,
J.
(
2022
).
Test–retest reliability of synchrony and metastability in resting state fMRI
.
Brain Sciences
,
12
(
1
),
66
. ,
[PubMed]

Competing Interests

Competing Interests: The authors have declared that no competing interests exist.

Author notes

Handling Editor: Bratislav Misic

This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. For a full description of the license, please visit https://creativecommons.org/licenses/by/4.0/legalcode.

Supplementary data