Alterations in topology, cost, and dynamics of gamma-band EEG functional networks in a preclinical model of traumatic brain injury

TL;DR

This study explores how traumatic brain injury alters brain network structures and may lead to epilepsy, using an animal model and gamma-band EEG data.

Contribution

The paper introduces a novel computational normalization method for network cost and identifies gamma-band network reorganization in TBI for the first time.

Findings

TBI alters small-world propensity and rich-club structure in gamma-band functional networks.

Post-injury hyperexcitability is observed as increased cortical sound-evoked response amplitudes.

Altered network topology is linked to potential mechanisms of posttraumatic epilepsy.

Abstract

Traumatic brain injury (TBI) is a major cause of disability leading to multiple sequelae in cognitive, sensory, and physical domains, including posttraumatic epilepsy. Despite extensive research, our understanding of its impact on macroscopic brain circuitry remains incomplete. We analyzed electrophysiological functional connectomes in the gamma band from an animal model of blast-induced TBI over multiple time points after injury. We revealed differences in small-world propensity and rich-club structure compared with age-matched controls, indicating functional reorganization following injury. We further investigated cost-efficiency trade-offs, propose a computationally efficient normalization procedure for quantifying the cost of spatially embedded networks that controls for connectivity strength differences, and observed dynamic changes across the injury timeline. To explore potential links between altered network topology and epileptic activity, we employed a brain-wide computational model of seizure dynamics and attribute brain reorganization to a homeostatic mechanism of activity regulation with the potential unintended consequence of driving generalized seizures. Finally, we demonstrated post-injury hyperexcitability that manifests as an increase in sound-evoked response amplitudes at the cortical level. Our work characterizes, for the first time, gamma-band functional network reorganization in a model of brain injury and proposes potential causes of these changes, thus identifying targets for future therapeutic interventions.

Traumatic brain injury (TBI) is a prevalent neurological disorder caused by factors such as accidents, contact sports, or military conflict. While animal studies have provided detailed insights into the molecular and cellular mechanisms underlying TBI, and clinical research has revealed its effects on large-scale brain network organization, these findings often lack integration. Our research seeks to address this gap by examining large-scale brain activity in an animal model of TBI. By analyzing brain network structures, we identify reorganization, highlighting novel connections to weight gain and posttraumatic epilepsy. Further exploration using this model may deepen our understanding of cross-scale mechanisms and inform the development of therapeutic interventions.

Figures

Top view schematic of the shock tube, including subject positioning, and of the electrode array. (A) The subject is positioned perpendicular to the incoming pressure shockwave, while the rest of the body besides the head is protected. The arrow illustrates that the subject’s head is allowed to move freely. The rest of the body besides the head is protected by a plastic tube with internal padding. (B) Macroscopic neural activity was recorded via epicranial EEG electrodes covering a large area of the skull. Adapted from Costa, G. (2020). Rat from the top. Zenodo. https://doi.org/10.5281/zenodo.3926343, licensed under CC BY 4.0.

Illustration of our cost normalization procedure. For this example, we consider a simple, fully connected, three-node network with node names 1, 2, and 3. Furthermore, we consider this network to be spatially embedded, that is, nodes have physical distances from one another, that are of lengths a, b, and c, with a < b < c. We also assume an actual (observed) network where the three edges have weights reflected on edge width (leftmost network). We define the nonnormalized network cost as the sum of products of distance and weight over all network edges. Our normalization method pertains to finding permutations of the edge weights such that the permuted networks have the maximum or minimum cost attainable given the edge distances and weights (middle and rightmost networks). For this toy network, the maximum cost is achieved through a clockwise permutation of all edges (middle network), while the minimum cost is achieved by swapping edges {1,2} and {2,3}. In the general case, finding these global optima is performed analytically and, thus, comprises a computationally efficient operation. Finally, the normalized cost of the observed network is defined as the fractional deviation of the observed cost from the minimum cost, that is, the deviation from minimum cost divided by the range of attainable costs (bottom schematic). Thus, our normalized measure is bounded in the interval [0, 1].

Decreased global centrality for the gamma connectivity networks chronically after injury. (A) Global betweenness centrality for all groups. (B) Global eigenvector centrality for all groups. All panels: Permutation-based two-way ANOVA (Group × Time Postblast with interaction, 10,000 permutations). Symbols denote p values for main effect of group. #: p < 0.05 (blast, 1 month: n = 13; blast, 3 months: n = 20; sham, 1 month: n = 8; sham, 3 months: n = 9).

Decreased deviation from regular network clustering for the gamma connectivity networks after injury. (A) SWP for all groups. (B) Fractional deviation from the mean clustering coefficient of an equivalent regular network for all groups. (C) Fractional deviation from the mean path length of an equivalent random network for all groups. All panels: permutation-based two-way ANOVA (Group × Time Postblast with interaction, 10,000 permutations). Symbols denote p values for main effect of group. #: p < 0.05 (blast, 1 month: n = 13; blast, 3 months: n = 20; sham, 1 month: n = 8; sham, 3 months: n = 9).

Loss of rich-club structure for the gamma connectivity networks chronically after injury. Permutation-based two-way ANOVA (Group × Time Post-Blast with interaction, 10,000 permutations) with FDR correction for multiple comparisons (n = 31 values of the richness parameter). Symbols denote corrected p values for main effect of group. *: p < 0.05 (blast, 1 month: n = 13; blast, 3 months: n = 20; sham, 1 month: n = 8; sham, 3 months: n = 9).

Dynamic changes in cost of gamma functional networks chronically after injury. (A) Nonnormalized network cost, similar to Roy et al. (2017). (B) Normalized network cost introduced here. All panels: permutation-based two-way ANOVA (Group × Time Postblast with interaction, 10,000 permutations). Symbols denote p values for main effect of group (#), interaction (£), or two-sample t test (*). £: p < 0.05; ##: p < 0.01; *: p < 0.001 (blast, 1 month: n = 13; blast, 3 months: n = 20; sham, 1 month: n = 8; sham, 3 months: n = 9).

Opposite correlation signs between normalized network cost and normalized weight gain between sham and injury groups in the chronic phase. Spearman correlation values and associated p values are reported. Lines correspond to best fit lines using least squares (red line: blast groups; blue line: sham groups). Blast, 1 month: n = 13; blast, 3 months: n = 20; sham, 1 month: n = 8; sham, 3 months: n = 9.

Dynamic changes in cost-efficiency of gamma functional networks chronically after injury. (A) Segragative cost-efficiency. (B) Integrative cost-efficiency. All panels: permutation-based two-way ANOVA (Group × Time Postblast with interaction, 10,000 permutations). Symbols denote p values for interaction effect (£) or two-sample t test (*). £, *: p < 0.05; ££, : p < 0.01; *: p < 0.001 (blast, 1 month: n = 13; blast, 3 months: n = 20; sham, 1 month: n = 8; sham, 3 months: n = 9).

AEP waveforms for all groups—broadband clicks. Grand average waveforms for the 1-month (top) and 3-month groups (bottom). Shaded areas denote the standard deviation across the group. Blast, 1 month: n = 14; blast, 3 months: n = 20; sham, 1 month: n = 8; sham, 3 months: n = 10.

Increase in AEP N1 magnitudes 1 month after blast in response to broadband clicks. Permutation-based two-way ANOVA (Group × Time Postblast with interaction, 10,000 permutations). Symbols denote p values for an interaction effect (£) or two-sample t test (*). : p < 0.01; £££: p < 0.001 (blast, 1 month: n = 14; blast, 3 months: n = 20; sham, 1 month: n = 8; sham, 3 months: n = 10).

Increased network ictogenicity and seizure probability for the gamma functional networks 3 months post-injury. (A) BNI. Permutation-based two-way ANOVA (Group × Time Postblast with interaction, 10,000 permutations). Symbols denote p values for interaction effect (£) or two-sample t test (*). £, *: p < 0.05 (blast, 1 month: n = 13; blast, 3 months: n = 20; sham, 1 month: n = 8; sham, 3 months: n = 9). (B) Seizure probability (Psz) as a function of excitability (I0).

Decreased NI for the gamma-band functional networks chronically after injury. (A) Topoplots depicting mean group NI. (B) Mean NI. (C) Standard deviation of NI. (B, C) Permutation-based two-way ANOVA (Group × Time Postblast with interaction, 10,000 permutations). Symbols denote p values for the main effect of group. #: p < 0.05 (blast, 1 month: n = 13; blast, 3 months: n = 20; sham, 1 month: n = 8; sham, 3 months: n = 9).

Tables

Table 1.. Parameters of the computational model of brain ictogenicity

Parameter	Meaning	Value
I 0	Node excitability	[−1.7, −0.5]
K	Global scaling of coupling strength	10
A	Connectivity matrix	dwPLI values scaled to [0, 1]
N	Network size	32
σ	Noise level	6
T	Simulation steps	4 ⋅ 106
δt	Simulation time step	10−2

Full text

INTRODUCTION

Traumatic brain injury (TBI) affects millions annually (Maas et al., 2017) with post-injury phenotypes spanning sensory (Lew, Cifu, Crowder, & Grimes, 2012), physical (Golub & Reddy, 2022; Reid, Bragin, Giza, Staba, & Engel, 2016), cognitive (Adhikari et al., 2023; Silver, McAllister, & Arciniegas, 2009), and psychological (Silver et al., 2009) domains. Blast-induced TBI (bTBI) is prevalent in war zones (Hoge et al., 2008; Rona et al., 2012) and has been coined the signature injury of warfare (Martin, Lu, Helmick, French, & Warden, 2008), affecting both military personnel and civilians (Sattin, Sasser, Sullivent, & Coronado, 2008; Wade et al., 2008).

To date, a comprehensive understanding of how macroscopic TBI phenotypes emerge from the underlying pathophysiology remains elusive (Carron, Alwis, & Rajan, 2016). At the microscopic scale, Cantu et al. (2015) revealed signatures of an imbalance of excitation and inhibition (EI) as a precursor of posttraumatic epilepsy, a documented sequela of injury (Verellen & Cavazos, 2010). Indeed, many studies have highlighted the vulnerability of inhibitory neurons to TBI (Cantu et al., 2015; Frankowski, Kim, & Hunt, 2019; Hunt, Scheff, & Smith, 2011; Lowenstein, Thomas, Smith, & McIntosh, 1992; Nichols, Bjorklund, Newbern, & Anderson, 2018; Toth, Hollrigel, Gorcs, & Soltesz, 1997). Crucially, EI imbalance can lead to maladaptive plasticity as loss of inhibition and excessive firing can induce synaptic strengthening, and if uncontrolled, this process can lead circuits into a pathological regime (Chen, Epstein, & Stern, 2010). Recent studies have begun to map the multiscale reorganization of cortical circuits (Chen et al., 2010; Frankowski et al., 2022; Hunt et al., 2011; Vascak, Jin, Jacobs, & Povlishock, 2018), while Frankowski et al. (2022) identified a widespread restructuring of inhibitory networks following focal TBI.

Yet, another consequence of EI imbalance and excessive neuronal firing following TBI is metabolic stress. Restoring altered ionic gradients requires sodium-potassium pump activity, leading to elevated adenosine triphosphate (ATP) consumption, and this heightened ATP demand, coupled with mitochondrial dysfunction (Arun et al., 2013; Xiong, Gu, Peterson, Muizelaar, & Lee, 1997), can result in an energy crisis within neural tissue. This may disrupt the brain’s balance between energy cost and efficiency (Bullmore & Sporns, 2012), typically associated with “small-world” topology in brain networks. This topology combines densely interconnected clusters (low energetic costs) and short paths linking any two nodes (geodesic distance), ensuring efficient information transfer (Watts & Strogatz, 1998) as metabolic demands shape brain networks to balance efficient information transfer and minimal wiring cost (Bullmore & Sporns, 2012). Metabolic stress following TBI is expected to disrupt this balance, and studies have shown deviations from small-world topology in this context (Caeyenberghs et al., 2014; Cao & Slobounov, 2010; Pandit et al., 2013). However, an underexplored avenue involves conceptualizing functional brain networks as spatially embedded, considering physical rather than geodesic node distances and calculating network cost, assuming that strong or long-distance coupling incurs higher costs (Roy et al., 2017; van den Heuvel, Kahn, Goñi, & Sporns, 2012). This approach characterizes the cost-efficiency of network organization, aligning with theories on the emergence of small-world organization in biological networks (Bullmore & Sporns, 2012).

Another critical feature of brain networks is the presence of a highly interconnected core of hub nodes, often referred to as the “rich-club” structure, which supports computations and synchronization across the entire network (Kim & Min, 2020; van den Heuvel et al., 2012). Due to its strong coupling density, this functional backbone incurs elevated costs (van den Heuvel et al., 2012) and could also become disrupted in the context of a TBI-induced metabolic crisis. Crucially, insights from computational modeling suggest that this topology can influence the propensity of a network to generate seizures, while the contribution of individual nodes to this behavior is influenced by their centrality, that is, their hub-like characteristics (Lopes et al., 2017). These network measures might therefore also be disrupted following TBI, in the context of EI imbalance.

Despite progress in TBI research, powerful analytical tools like network theory (Bassett & Sporns, 2017; Brynildsen, Rajan, Henderson, & Bassett, 2023; Farahani, Karwowski, & Lighthall, 2019; Heiney et al., 2021; Sharp, Scott, & Leech, 2014) and dynamic brain modeling (Pathak, Roy, & Banerjee, 2022; Sharp et al., 2014), along with translationally relevant recording modalities, have been underappreciated in preclinical models. In this work, we address this gap by further analyzing data from a previous study using a rat model of bTBI (May et al., 2024). We experimentally simulated blast pressure (Figure 1A) and recorded brain activity under anesthesia at 1 and 3 months post-injury using on-skull, multielectrode EEG (Figure 1B). Consistent with previous findings of hyperexcitability and altered connectivity post-injury, we observed elevated broadband spectral power and functional hyperconnectivity in the gamma frequency band (25–80 Hz), suggesting potential reorganization of distributed brain networks, associated with decreased inhibitory interneuron density. In the current work, we first analyze small-world propensity (SWP), a key topological property disrupted by TBI (Pandit et al., 2013), and hypothesize an energetic driver of this disruption reflected in network cost (Bullmore & Sporns, 2012) and reduction of costly topological features like rich-club structure. We introduce a novel metric for quantifying network cost, normalized to equivalent null models. SWP and rich-club structure were disrupted post-bTBI, while network cost and its efficiency trade-off showed a biphasic response: low cost at 1 month post-injury, high cost at 3 months, undetectable with the nonnormalized metric (Figure 2). This cost negatively correlated with weight gain, used as a proxy for the post-injury metabolic state. We also provide evidence for hyperexcitability in our model, reflected by increased auditory evoked potential amplitudes. Finally, we investigate the altered network topology’s impact on whole brain dynamics through computational modeling, observing increased seizure risk at the later timepoint, driven by most network nodes.

Top view schematic of the shock tube, including subject positioning, and of the electrode array. (A) The subject is positioned perpendicular to the incoming pressure shockwave, while the rest of the body besides the head is protected. The arrow illustrates that the subject’s head is allowed to move freely. The rest of the body besides the head is protected by a plastic tube with internal padding. (B) Macroscopic neural activity was recorded via epicranial EEG electrodes covering a large area of the skull. Adapted from Costa, G. (2020). Rat from the top. Zenodo. https://doi.org/10.5281/zenodo.3926343, licensed under CC BY 4.0.

Illustration of our cost normalization procedure. For this example, we consider a simple, fully connected, three-node network with node names 1, 2, and 3. Furthermore, we consider this network to be spatially embedded, that is, nodes have physical distances from one another, that are of lengths a, b, and c, with a < b < c. We also assume an actual (observed) network where the three edges have weights reflected on edge width (leftmost network). We define the nonnormalized network cost as the sum of products of distance and weight over all network edges. Our normalization method pertains to finding permutations of the edge weights such that the permuted networks have the maximum or minimum cost attainable given the edge distances and weights (middle and rightmost networks). For this toy network, the maximum cost is achieved through a clockwise permutation of all edges (middle network), while the minimum cost is achieved by swapping edges {1,2} and {2,3}. In the general case, finding these global optima is performed analytically and, thus, comprises a computationally efficient operation. Finally, the normalized cost of the observed network is defined as the fractional deviation of the observed cost from the minimum cost, that is, the deviation from minimum cost divided by the range of attainable costs (bottom schematic). Thus, our normalized measure is bounded in the interval [0, 1].

RESULTS

Blast-Induced Changes in the Topology of Resting Functional Networks

Betweenness centrality (number of shortest paths in a network that traverse that node) was reduced at the global level after blast (Figure 3A; F(1, 46) = 4.89, p < 0.05, main effect of group). Eigenvector centrality (the centrality of a node as a function of the centrality of its neighbors) showed a trend towards increase for the injury group (Figure 3B; p > 0.05). Similarly, normalized betweenness centrality was not different between groups (Supporting Information Figure S1). No timepoint or interaction effects were detected for these measures.

Decreased global centrality for the gamma connectivity networks chronically after injury. (A) Global betweenness centrality for all groups. (B) Global eigenvector centrality for all groups. All panels: Permutation-based two-way ANOVA (Group × Time Postblast with interaction, 10,000 permutations). Symbols denote p values for main effect of group. #: p < 0.05 (blast, 1 month: n = 13; blast, 3 months: n = 20; sham, 1 month: n = 8; sham, 3 months: n = 9).

In terms of whole-network measures, we observed a trend toward increased SWP in the blast group compared with the sham group (Figure 4A). The small-world topology combines the high clustering coefficients of regular and the short path lengths of random networks (Muldoon, Bridgeford, & Bassett, 2016; Watts & Strogatz, 1998). Hence, we dissected SWP into its constituents by analyzing the deviation of the mean clustering coefficient from that of an equivalent regular network and the deviation of the mean path length from that of an equivalent random network (ΔC and Δ_L, _ respectively; see the Materials and Methods section for more details) and revealed a reduced ΔC for the blast group at 1 month post-injury. This was accompanied by a larger deviation from the path length of an equivalent random network compared with the other experimental groups (Figure 4B–C). This trend was confirmed for ΔC (F(1, 46) = 4.20, p < 0.05, main effect of injury group) but did not reach significance for ΔL (p > 0.05). The presence of low values of ΔL across most of our cohort indicates that path lengths were identical to those of equivalent random networks. This result was robust to binarizing and thresholding networks across a wide range of thresholds (Supporting Information Figure S2). No timepoint or interaction effects were observed for these metrics.

Decreased deviation from regular network clustering for the gamma connectivity networks after injury. (A) SWP for all groups. (B) Fractional deviation from the mean clustering coefficient of an equivalent regular network for all groups. (C) Fractional deviation from the mean path length of an equivalent random network for all groups. All panels: permutation-based two-way ANOVA (Group × Time Postblast with interaction, 10,000 permutations). Symbols denote p values for main effect of group. #: p < 0.05 (blast, 1 month: n = 13; blast, 3 months: n = 20; sham, 1 month: n = 8; sham, 3 months: n = 9).

Finally, the presence of a rich-club structure was assessed using the normalized rich-club coefficient (Figure 5). Over a range of richness parameter values, a significant main effect of injury group was found with blast groups having diminished rich-club structure compared with sham groups.

*Loss of rich-club structure for the gamma connectivity networks chronically after injury. Permutation-based two-way ANOVA (Group × Time Post-Blast with interaction, 10,000 permutations) with FDR correction for multiple comparisons (n = 31 values of the richness parameter). Symbols denote corrected p values for main effect of group. : p < 0.05 (blast, 1 month: n = 13; blast, 3 months: n = 20; sham, 1 month: n = 8; sham, 3 months: n = 9).

Blast-Induced Changes in Network Cost

The small-world topology is believed to achieve a cost-efficiency trade-off (Bullmore & Sporns, 2012), while the rich-club structure incurs high costs to a network (van den Heuvel et al., 2012). To identify putative drivers of the observed group differences in these features, we turn to network cost. Cost was first calculated similarly to Roy et al. (2017) as the inner product of euclidean distances and connectivity values across node pairs. However, this metric could be sensitive to group differences in connectivity. Hence, we propose a normalization procedure based on equivalent null models (see the Materials and Methods section). Briefly, we calculate lower and upper bounds of network cost for two inputs: a set of euclidean distances and a set of connectivity values. The normalized cost is then defined as the fractional deviation of the empirically observed cost from the lower cost bound. This measure is robust to differences in connectivity strength and can be interpreted as reflecting a “strategy” of resource (connectivity strengths) allocation to a system of prespecified properties (node distances), a potential plasticity-based phenomenon in the post-TBI brain. Our measure is computationally tractable, involving only a sorting operation of the distance and connectivity sets.

Network cost was elevated for the blast groups (F(1, 46) = 7.35, p < 0.01, main effect of group; Figure 6A) when using the nonnormalized cost metric. However, when considering our normalized cost measure, we observed a significant interaction effect (F(1, 46) = 16.85, p < 0.05; Figure 6B). Post hoc pairwise comparisons revealed a decreased normalized cost for the blast group at 1 month compared with the sham group at the same timepoint (t(30.37) = −4.60, p < 0.001, 95% CI [−0.23, −0.09]), and an elevated normalized cost for the blast group at 3 months compared with the blast group at 1 month (t(18.32) = −5.28, p < 0.001, 95% CI [−0.21, −0.09]). The blast group at 3 months was not found to be different to the sham group at 3 months.

Dynamic changes in cost of gamma functional networks chronically after injury. (A) Nonnormalized network cost, similar to Roy et al. (2017). (B) Normalized network cost introduced here. All panels: permutation-based two-way ANOVA (Group × Time Postblast with interaction, 10,000 permutations). Symbols denote p values for main effect of group (#), interaction (£), or two-sample t test (). £: p < 0.05; ##: p < 0.01; **: p < 0.001 (blast, 1 month: n = 13; blast, 3 months: n = 20; sham, 1 month: n = 8; sham, 3 months: n = 9).

Interestingly, our normalized cost metric was found to be significantly negatively correlated with normalized weight gain for the blast groups (Figure 7, red hue symbols; Spearman ρ = −0.52, 95% CI [−0.76, −0.17], p < 0.01), but positively correlated for the sham groups (Figure 7, blue hue symbols; Spearman ρ = 0.71, 95% CI [0.31,0.89], p < 0.01). This pattern was preserved when correlations were examined in all four subgroups independently except for the blast group measured at 1 month (Supporting Information Figure S3).

Opposite correlation signs between normalized network cost and normalized weight gain between sham and injury groups in the chronic phase. Spearman correlation values and associated p values are reported. Lines correspond to best fit lines using least squares (red line: blast groups; blue line: sham groups). Blast, 1 month: n = 13; blast, 3 months: n = 20; sham, 1 month: n = 8; sham, 3 months: n = 9.

Alterations in network cost following TBI could reflect broader changes in cost-efficiency trade-offs (Bullmore & Sporns, 2012). Hence, we further integrate our cost measure in cost-efficiency metrics using the fractional deviations of SWP (see the Materials and Methods section) as deviations from efficiency. Specifically, we interpret clustering deviation from the equivalent regular network as segregative inefficiency and path length deviation from the equivalent random network as integrative inefficiency. We define high cost-efficiency as having low cost and efficiency deviations from their optimal values (see the Materials and Methods section). By virtue of their constituent quantities being normalized to appropriate null models, these measures are also robust to differences in overall density and connectivity strength.

An interaction effect was found for both segregative (Figure 8A; F(1, 46) = 6.04, p < 0.05) and integrative (Figure 8B; F(1, 46) = 11.86, p < 0.01) cost-efficiency, and post hoc pairwise comparisons revealed increased cost-efficiencies for the blast group at 1 month compared to the sham group at the same timepoint (segregative: t(13.07) = 4.13, p < 0.001, 95% CI [0.11, 0.34], Figure 8A; integrative: t(17.18) = 3.05, p < 0.05, 95% CI [0.02, 0.12], Figure 8B), as well as compared with the blast group at 3 months (segregative: t(31) = 3.69, p < 0.01, 95% CI [0.08, 0.28], Figure 8A; integrative: t(30.25) = 3.49, p < 0.01, 95% CI [0.03,0.13], Figure 8B). Decreased integrative cost-efficiency was found for the blast group at 3 months, compared with the sham group at the same timepoint (t(22.47) = −2.56, p < 0.05, 95% CI [−0.11, −0.01], Figure 8B).

Dynamic changes in cost-efficiency of gamma functional networks chronically after injury. (A) Segragative cost-efficiency. (B) Integrative cost-efficiency. All panels: permutation-based two-way ANOVA (Group × Time Postblast with interaction, 10,000 permutations). Symbols denote p values for interaction effect (£) or two-sample t test (). £, *: p < 0.05; ££, **: p < 0.01; **: p < 0.001 (blast, 1 month: n = 13; blast, 3 months: n = 20; sham, 1 month: n = 8; sham, 3 months: n = 9).

Blast-Induced Changes in Auditory Evoked Potentials

As supporting evidence for the presence of a disrupted EI balance in our injury model, we probed the excitability of central auditory regions using evoked potentials. For broadband clicks, grand averages of evoked potential waveforms indicated increased amplitudes for the blast group, particularly at the 1 month timepoint (Figure 9). Analysis of the first negative component (N1) amplitude showed a significant interaction (F(1, 48) = 11.93, p < 0.001; Figure 10), with post hoc pairwise comparisons revealing a significant increase of magnitude for the blast group at 1 month compared with the sham group at the same timepoint (t(15.64) = −4.05, p < 0.01, 95% CI [−10.22, −3.19]) and to the blast group at 3 months (t(16.57) = −3.82, p < 0.01, 95% CI [−9.93, −2.85]). Similarly, increased amplitudes for the blast group were found for most other stimuli employed (Supporting Information Figures S6–S8).

AEP waveforms for all groups—broadband clicks. Grand average waveforms for the 1-month (top) and 3-month groups (bottom). Shaded areas denote the standard deviation across the group. Blast, 1 month: n = 14; blast, 3 months: n = 20; sham, 1 month: n = 8; sham, 3 months: n = 10.

Increase in AEP N1 magnitudes 1 month after blast in response to broadband clicks. Permutation-based two-way ANOVA (Group × Time Postblast with interaction, 10,000 permutations). Symbols denote p values for an interaction effect (£) or two-sample t test (). *: p < 0.01; £££: p < 0.001 (blast, 1 month: n = 14; blast, 3 months: n = 20; sham, 1 month: n = 8; sham, 3 months: n = 10).

Blast-Induced Changes in Ictogenicity and the Influence of Topology

To decipher the role of altered network topology on seizure dynamics in the presence of hyperexcitability, we employed the brain ictogenicity framework (Lopes et al., 2017, 2019; Tait et al., 2021). Briefly, the functional connectivity values for each pair of nodes in sensor space are used as edge weights to define networks as for the graph-theoretical analysis above, while now, each node is imbued with dynamics that unfold over time. These dynamics are governed by a phase oscillator model operating in either a resting or a rotating dynamical regime, interpreted as resting or seizure dynamics, respectively (see the Materials and Methods section). Stronger coupling would result in higher ictogenicity and vice versa (e.g., see Figure 4a in Lopes et al., 2017). A control parameter I0 tunes the baseline excitability of individual nodes, on top of which, random noise and input currents from other nodes are summed to form the total current received by a node (Equation 12). This value is assigned to be equal across all nodes following previous work (Lopes et al., 2017; Tait et al., 2021). Subsequently, the network activity is simulated while parametrically varying this value across a range of negative values that do not drive the network into persistent epileptic dynamics (Tait et al., 2021). See the Connectivity-Informed Whole-Brain Computational Modeling section and Table 1 for more details on the simulation study. Outcome measures include seizure probability (Psz), brain network ictogenicity (BNI), and node ictogenicity (NI), the latter quantifying each node’s contribution to seizure probability in the whole network.

Globally, an interaction effect was found for BNI (Figure 11A; F(1, 46) = 4.47, p < 0.05). Post hoc pairwise comparisons revealed an increased BNI for the blast group at 3 months post-injury compared with the corresponding sham group (Figure 11A; t(15.84) = 2.78, p < 0.05, 95% CI [0.03, 0.23]). The plot of seizure probability (Psz) over excitability (I0) values revealed that the increased BNI for the blast group at 3 months was due to the increased seizure probability at intermediate to high excitability values (Figure 11B).

Increased network ictogenicity and seizure probability for the gamma functional networks 3 months post-injury. (A) BNI. Permutation-based two-way ANOVA (Group × Time Postblast with interaction, 10,000 permutations). Symbols denote p values for interaction effect (£) or two-sample t test (). £, : p < 0.05 (blast, 1 month: n = 13; blast, 3 months: n = 20; sham, 1 month: n = 8; sham, 3 months: n = 9). (B) Seizure probability (Psz) as a function of excitability (I0).

NI quantifies the contribution of individual nodes to ictogenicity (Figure 12A). A significant effect of injury group was revealed for mean NI over nodes with decreased mean NI for the blast groups compared with the shams (Figure 12B; main effect of group, F(1, 46) = 4.40, p < 0.05). The standard deviation of NI over nodes showed a trend toward reduction for the blast groups across timepoints (Figure 12C).

Decreased NI for the gamma-band functional networks chronically after injury. (A) Topoplots depicting mean group NI. (B) Mean NI. (C) Standard deviation of NI. (B, C) Permutation-based two-way ANOVA (Group × Time Postblast with interaction, 10,000 permutations). Symbols denote p values for the main effect of group. #: p < 0.05 (blast, 1 month: n = 13; blast, 3 months: n = 20; sham, 1 month: n = 8; sham, 3 months: n = 9).

As expected by the hub-like characteristics of a node as a proxy for its influence on network dynamics (Lopes et al., 2017), NI of individual nodes strongly correlated with centrality measures (Supporting Information Figure S4A, Spearman ρ = 0.6, 95% CI [0.57, 0.64], p < 0.001; Supporting Information Figure S4B, Spearman ρ = 0.85, 95% CI [0.83, 0.87], p < 0.001). Furthermore, mean NI correlated with the maximum value of the rich-club coefficient across richness parameter values (Supporting Information Figure S5; β = 0.02, 95% CI [0.015, 0.023], t(48) = 8.74, p < 0.001, adjusted R^2^ = 0.611).

DISCUSSION

In this study, we observed changes in gamma-band functional network topology post-injury. Locally, betweenness centrality decreased. Globally, networks showed increased SWP due to higher clustering, and a pronounced loss of rich-club structure. Network costs decreased at 1 month post-injury, recovering by the third month, a pattern only discernible with our normalized metric. Network cost correlated with post-injury weight gain, positively in sham and negatively in blast cohorts. We found evidence of post-injury hyperexcitability through increased auditory evoked potential amplitudes. Computational modeling revealed increased network ictogenicity at 3 months post-injury and a reduced role of individual nodes in epileptogenesis, significantly correlated with centrality and rich-club structure.

Our findings align with prior research. Altered small-worldness has been reported in EEG-based TBI studies (Cao & Slobounov, 2010), and MRI-based studies have observed increased path lengths (Caeyenberghs et al., 2014; Pandit et al., 2013). Changes within the gamma frequency range suggests a plastic reorganization of inhibitory networks (Hunt et al., 2011; Vascak et al., 2018), as these neurons shape activity in this band (Cardin et al., 2009) and are impacted by injury in our model (May et al., 2024). Supporting this, Frankowski et al. (2022) documented increased local and reduced long-range input to inhibitory interneurons, even in distant areas from the injury locus. We observed this effect along with a disruption in the rich-club structure, consistent with a diffusion imaging study (Verhelst, Vander Linden, De Pauw, Vingerhoets, & Caeyenberghs, 2018). Magnetoencephalography (MEG) studies have shown changes in rich-club hub representation (Antonakakis, Dimitriadis, Zervakis, Papanicolaou, & Zouridakis, 2017a) and hyperactivation of rich-club networks (Antonakakis, Dimitriadis, Zervakis, Papanicolaou, & Zouridakis, 2017b), echoing fMRI observations (Hillary et al., 2014; Li et al., 2022).

To identify drivers of altered small-worldness, we investigated the cost-efficiency trade-off thought to be afforded by the small-world topology. This trade-off is usually quantified using cost metrics incorporating connection distance and strength (Roy et al., 2017; van den Heuvel et al., 2012). A recent modeling study (Achterberg, Akarca, Strouse, Duncan, & Astle, 2023) suggests that artificial neural networks optimizing performance and cost converge to a small-world architecture. However, topological measures should be compared with appropriate null models, as confounding factors such as weight distribution may spuriously drive group differences in topological measures (Váša & Mišić, 2022), an important confound in TBI-induced connectivity strength differences. We integrated these arguments by introducing a novel, computationally efficient normalization procedure for calculating cost in weighted networks based on equivalent null models. This revealed a complex post-TBI network cost trajectory, indiscernible with the standard metric. The trajectory of diminished cost in early chronic stages followed by an increase in later stages might suggest a temporally evolving metabolic perturbation post-injury. Literature connects TBI with such aberrations (Bergsneider et al., 2000; Kawamata, Katayama, Hovda, Yoshino, & Becker, 1992; McGill et al., 2022; Shiga et al., 2006; Sullivan, Keller, Mattson, & Scheff, 1998; Vallez Garcia, Otte, Dierckx, & Doorduin, 2016; Xiong et al., 1997), noting a multiphasic response typified by initial hypermetabolism followed by hypometabolism and potential recovery (Koenig & Dulla, 2018). Alternatively, the increased cost for the blast group at 3 months relative to the blast group at 1 month might reflect a shift toward more efficient information processing by allocating high-weight edges to distant node pairs. This could reflect a transition away from a low-cost, low-efficiency topology that does not satisfy the behavioral needs of the organism. Importantly, post-TBI cerebral glucose metabolism measured by positron emission tomography (PET) has been positively correlated with full-scale intelligence quotient scores (Kato et al., 2007), corroborating the existence of trade-offs between energy conversion and cognition.

Interestingly, there was a reversal in the correlation direction between our proposed network cost metric and weight gain across groups. Our metric, interpreted as a proxy for resource allocation strategy (lower cost reflecting a lower part of the cost range given network distances and connectivity strengths), suggests potential differences in energy intake due to TBI. In a typical brain, increased energy supply might allow elevated network costs, promoting enhanced computational capabilities. In contrast, putative TBI-induced metabolic challenges could drive increased food intake and employ a more frugal resource allocation strategy, reducing network costs. This would make weight gain inversely associated with network cost, as shown here. The restored costs at 3 months post-injury could reflect a compensatory response. Specifically, the reduced cost at 1 month might curtail computational capabilities enough to disrupt normal function, necessitating a shift to a more resource-intensive configuration to offset these challenges. These findings on cost dynamics align with cost-efficiency analyses, supporting this interpretation.

The localization of these results in the gamma frequency band could also provide converging evidence for a metabolic driver in the observed effects, as fast-spiking inhibitory interneurons’ activity is energetically demanding (Kann, Papageorgiou, & Draguhn, 2014). Kann (2016) hypothesized that these neural subtypes might be preferentially affected by TBI through an energetic insult. Our present and previous (May et al., 2024) findings support this hypothesis and extend it to the macroscopic spatial scale and reorganization of brain-wide networks as measured by EEG electrodes.

The rich-club structure is postulated to function as a densely interconnected core that integrates information processed by peripheral components (Kim & Min, 2020; van den Heuvel et al., 2012). Hence, this functional backbone incurs a high cost (van den Heuvel et al., 2012). Consequently, a metabolic perspective could also explain the observed post-injury disruption of this topological feature observed in this study. However, a complementary perspective from our computational modeling emphasizes neural dynamics. The rich-club’s capacity to facilitate activity propagation can be a double-edged sword. While beneficial for information transfer in a healthy brain, this could, in the backdrop of TBI and its associated hyperexcitability, induce uncontrolled activation cascades, possibly culminating in seizures. Lopes et al. (2017) demonstrated that targeted elimination of rich-club nodes curtails network ictogenicity in silico, and this intervention corresponded with favorable outcomes in postresective surgery of epileptic patients. Thus, the loss of rich-club organization post-TBI might signify an adaptive strategy to counteract excessive activity spread and mitigate seizure risks, especially amid EI imbalances. Importantly, our interpretation rests on the assumption that changes in edge-level features such as the topological properties investigated here occur in response to node-level changes such as disrupted EI balance, in order to mitigate their effects on the rest of the network. Our findings indicate an opposing effect of injury on global and local ictogenicity with interesting implications. On one hand, network-wide seizure probability remains unchanged at 1 month but is increased by the third month post-injury, compared with shams at the same timepoint. In parallel, the mean contribution of each node as well as the heterogeneity in these contributions decreases for the blast group, suggesting that all nodes equally drive the network toward this more prominent epileptic regime. This could reflect the interplay of local and global topological changes both between injury groups and across the injury timeline that are shared (like rich-club organization) or unique (small-worldness or cost). Finally, the switch-like behavior observed at the later timepoint, combined with the decreased contribution of individual nodes to seizures throughout the whole network, could suggest a potential substrate for generalized seizures rather than focal seizures, as observed in a study using data from clinical populations with these two types of epilepsy (Lopes et al., 2019). The reduction in individual node contribution, as well as its correlation with topological measures, lends support to the interpretation of an adaptive loss of rich-club structure. Intriguingly, a similar pattern in global ictogenicity and node-level contributions has been identified in research focusing on Alzheimer’s disease (Tait et al., 2021), suggesting a potentially general adaptive mechanism applicable across diverse neurological conditions.

The influence of injury on evoked potentials was prominent, pointing to hyperexcitability (Bruyns-Haylett et al., 2017) in cortical auditory regions and echoing increased amplitudes in pediatric epileptic patients (Zgorzalewicz, 2006). This pattern, combined with disruptions to inhibitory interneurons across multiple brain regions in the same model (May et al., 2024), reinforces our interpretation of the observed network differences as indicative of a disrupted EI balance (Carron et al., 2016). Although we cannot rule out a hearing loss–related driver of this phenotype (Popelar, Grecova, Rybalko, & Syka, 2008; Rumschlag & Razak, 2021), recent research in a Fragile X syndrome model (Jonak, Lovelace, Ethell, Razak, & Binder, 2020) noted a similar trend of increased auditory evoked potential amplitudes. Given that this syndrome shares a mechanistic substrate of hyperexcitability, it is likely that a significant portion of the effects observed here stem from bTBI. Future studies could further explore the interplay between network topology and node excitability on ictogenicity by incorporating heterogeneous excitability across network nodes, for example, by using a proxy measure such as the power spectral exponent (Gao, Peterson, & Voytek, 2017).

Overall, the graph theoretical characterization of functional networks has yielded a nuanced view of injury consequences. Supplemented by a new methodological contribution in calculating network cost, which shows promise of increased sensitivity to experimental effects while controlling for potential confounds, this approach has integrated electrophysiology with a metabolic proxy post-TBI. Computational modeling results suggest that these alterations might be the brain’s response as part of a homeostatic mechanism to maintain dynamics post-injury. The study also reiterates the dynamic nature of injury progression and the potential role of the inhibitory plasticitome (McFarlan et al., 2023), aligning with the hypothesis that TBI reopens the early developmental critical period where the EI balance has not yet stabilized, and plastic changes occur more frequently and intensely (Masri, 2020).

Recording brain activity during an anesthetized state is a major limitation of the current work, and future investigations should record activity in awake behaving rodents (Choi, Koch, Poppendieck, Lee, & Shin, 2010; Jonak, Lovelace, Ethell, Razak, & Binder, 2018; Jonak et al., 2020). We do, however, expect some of our results to transfer to the awake state, based on evidence on the robustness of gamma-band connectivity differences to anesthetic state when using urethane (Mondino et al., 2024) as well as the similarity of fMRI-based functional connectivity patterns to those observed during wakefulness (Paasonen, Stenroos, Salo, Kiviniemi, & Gröhn, 2018). Importantly, unlike other anesthetics that target specific receptors (e.g., isoflurane that targets GABAergic receptors (Jones, Brooks, & Harrison, 1992) and influences gamma-band oscillations (Hudetz, Vizuete, & Pillay, 2011)), urethane has a nonspecific mechanism of action (Hara & Harris, 2002), making it an appropriate choice for studies investigating EI balance. However, this anesthetic has been shown to modulate the excitability (Schumacher, Schneider, & Woolley, 2011), geometry of population activity (Kobak, Pardo-Vazquez, Valente, Machens, & Renart, 2019), and spike rate and representational capacity (Pachitariu, Lyamzin, Sahani, & Lesica, 2015) but not the activity structure (Harris et al., 2011) of auditory neurons or inhibitory gating in nonauditory regions (Wang, Ma, Wang, & Qin, 2018). Incorporating the collection of behavioral data into future study designs would also help the speculated functional significance of the electrophysiological changes observed in the current study to be elucidated. Moreover, the suspected metabolic disruptions in this model and their association with the proposed network cost metric should be validated by future studies using additional techniques able to directly interrogate brain metabolism such as PET scanning. Finally, using a longitudinal within-subjects design could improve sensitivity to experimental effects that dynamically unfold over the post-injury period. By benchmarking our cost metric using larger datasets from clinical cohorts, its applicability can be gauged against other energetically relevant neural measures (Chennu et al., 2017; Smailovic et al., 2020).

MATERIALS AND METHODS

Experimental Procedures

All animal experiments were conducted in compliance with the Home Office License ([Scientific Procedures] Act 1986 and European Union legislation). Briefly, male Sprague Dawley rats were used, housed in standard conditions. We employed a 2 × 2 between-subjects factorial design based on injury group (sham or blast) and post-injury timepoint (1 or 3 months). Blast exposure was conducted using a shock tube (Campos-Pires et al., 2023; Eftaxiopoulou et al., 2016; Kazezian, Yu, Ramette, Macdonald, & Bull, 2021; T.-T. Nguyen, 2016; T.-T. Nguyen et al., 2019; T.-T. N. Nguyen, Sory, Amin, Rankin, & Proud, 2018; T.-T. N. Nguyen, Wilgeroth, & Proud, 2014), and the injury was unilaterally delivered on the right side of isoflurane-anesthetized animals, isolated to the head, which was allowed to move freely (Figure 1A). Animals in the blast group received a pressure wave with peak overpressure of 232.5826(6.4486) kPa and positive phase duration of 1.3870(0.0427) ms (both values M(SD)), while the sham animals received an identical procedure minus the pressure wave. Electrophysiological recordings involved terminal anesthesia under urethane, and 32-channel flexible EEG multielectrode arrays and a SmartBox Pro/Radiens Allego system (all electrophysiological equipment by Neuronexus, Michigan, USA). Electrodes were implanted on the dorsal skull surface and span 11.7 mm on the rostrocaudal and 9.6 mm on the left-to-right axis, with a minimum interelectrode distance of 1.2 mm (Figure 1B; Neuronexus, Rat EEG Functional). Signals were referenced to the nuchal musculature and sampled at 30 kHz. More details can be found in the Supporting Information and in our previous work (May et al., 2024).

Analysis of Functional Connectomes

Data preprocessing and connectome construction.

Data preprocessing was performed using functions from the EEGLAB Toolbox (Delorme & Makeig, 2004) and custom scripts in MATLAB (Version 2019b). Briefly, signals were downsampled to 1 kHz and line noise was removed using the CleanLine plugin (Bigdely-Shamlo, Mullen, Kothe, Su, & Robbins, 2015; Mullen, 2012). This was followed by automatic artifact rejection (EEGLAB function clean_artifacts; less than three electrodes were rejected per recording and all accepted ones were used for further analysis) and an independent-component decomposition using the Infomax algorithm (function pop_runica) to separate artifactual components. Components visually inspected those containing primarily respiratory-, cardiac-, or muscle-related artifacts were conservatively removed from the data. Three animals were excluded due to noisy data. More details on preprocessing can be found in the Supporting Information.

Functional connectivity was assessed using the debiased weighted phase-lag index (dwPLI; Vinck, Oostenveld, van Wingerden, Battaglia, & Pennartz, 2011), using MATLAB and FieldTrip functions. dwPLI is a frequency-resolved, phase-based metric quantifying the consistency of phase differences between two time series, interpreted as the time delay for neural activity to reach the target region (Cohen, 2014). Importantly, phase differences are weighted by magnitude such that zero-phase connectivity (which might be attributable to volume-conduction in sensor-level analysis) is discarded (Vinck et al., 2011). Signals were segmented into 3-s windows and the cross-spectral density for each electrode pair was calculated using a 1.5-s sliding window and 50% overlap (MATLAB function cpsd). dwPLI values were calculated using the FieldTrip function ft_connectivity_wpli. Connectivity values for a frequency band are determined by taking the mean across the frequency bins specific to the gamma band (25–80 Hz).

Local topological measures.

Topological analysis of connectomes was conducted using custom scripts, MATLAB functions, functions from the Brain Connectivity Toolbox (BCT; Rubinov, Kötter, Hagmann, & Sporns, 2009), and publicly available code (Muldoon et al., 2016).

To characterize local structure, we quantified the importance of network nodes via centrality measures. Betweenness centrality (Matlab function centrality) is defined as the number of shortest paths in the network that transverse that node, while eigenvector centrality (BCT function eigenvector_centrality_und) measures the centrality (importance) of a node as a function of the centrality of its neighbors (Lohmann et al., 2010).

Global variants of node-specific measures were calculated as the mean across the network. Since betweenness centrality can be influenced by additive shifts of the weight distribution, a normalized measure was also computed by dividing the global centrality of the observed networks by the mean global centrality of 100 equivalent random networks.

Global topological measures.

To characterize global topology, we employed SWP (Muldoon et al., 2016), which quantifies the small-world characteristics of a network by positioning it along the continuum from regular to random (as in the original formulation from Watts & Strogatz, 1998) using appropriate null models. This measure quantifies for the observed network, the deviation of its clustering coefficient ΔC, from an equivalent regular network, and the deviation of its path length ΔL, from an equivalent random network and is computed according to the following equation:

$$[eqn]$$

The deviations are calculated as follows:

$$[eqn]$$ $$[eqn]$$

where obs, reg, rand refer to the observed, equivalent regular, and equivalent random networks, respectively.

To calculate the clustering coefficients, the measure by Onnela, Saramäki, Kertész, and Kaski (2005) is used:

$$[eqn]$$

where Ci is the clustering coefficient of node i, ki the number of edges adjacent to node i, and $[eqn]$ the weight of the edge between nodes i and j divided by the maximum network weight.

For path lengths, edge weights is defined as the functional coupling strength between adjacent nodes; hence, nodes connected by a strong edge can be interpreted as being in close proximity from a neural communication perspective. Hence, we interpret inverse edge weight as distance (Muldoon et al., 2016). Shortest paths are calculated using Dijkstra’s algorithm and is used to characterize the whole network:

$$[eqn]$$

where N is the number of nodes and dij the length of the shortest path connecting nodes i and j.

To generate equivalent random networks, observed network edges are randomly reassigned. For equivalent regular networks, edges are organized in descending order and are distributed along the subdiagonals of the adjacency matrix. These equivalent networks preserve the weight distribution, making them appropriate in the context of TBI.

Rich-club structure was assessed using the approach by van den Heuvel et al. (2013). Specifically, the rich-club coefficient was calculated for each value of the richness parameter defined as the node degree (BCT function rich_club_wu), and the same calculation was performed for 100 equivalent random networks (BCT function randmio_und_connected, which ensures that the randomized network maintains connectedness). Intuitively, the rich-club structure can be understood as the presence of strong coupling between the hub (i.e., rich in connections) nodes of a network. Finally, the normalized rich-club coefficient was calculated as the ratio of the observed coefficient to the mean of the coefficients of the random networks, for each value of the richness parameter.

A novel normalization procedure for cost quantification.

The concept of network cost has been particularly illuminating in the context of TBI. A notable study (Roy et al., 2017) determined network cost by computing the dot product of the functional connectivity matrix and the physical distance matrix, based on the premise that connections that are strongly coupled or spatially distant demand higher energetic expenditure.

In the context of TBI, as with the network topology measures previously discussed, variations in cost—given a network’s spatial configuration (physical distances between nodes)—might trivially arise simply due to disparities in overall connectivity strength. Consequently, to draw valid comparisons between networks, it is necessary to employ a suitable normalization method or a null model. To this end, we introduce the concept of total effective wiring cost (TEWC). This metric is derived by initially determining the cost of the observed network, following the methodology of Roy et al. (2017), where the cost is the dot product of edge weights and corresponding Euclidean distances between the adjacent nodes:

$$[eqn]$$

where di, j denoted the node distances and wi, j the corresponding edge weights.

Following this, two surrogate null networks are constructed via reassignment of edge weights, representing the extremes of the cost spectrum for the given configuration of the network: the minimum cost equivalent network and the maximum cost equivalent network. Constructed to incur either the minimal or maximal cost, these networks maintain the sets of physical distances and connection strengths between nodes, and their costs can be calculated using Equation 6. Subsequently, the TEWC is defined as the fractional deviation from the minimum cost, relative to the maximum cost attainable:

$$[eqn]$$

where the subscripts obs, min, max refer to the observed, maximum cost equivalent, and minimum cost equivalent networks, respectively.

This formulation can be interpreted in terms of a strategy of resource allocation (i.e., weights) across a set of nodes with fixed physical distances. Thus, it is suitable to quantify potential network reorganization schemes that may be implemented after injury by plasticity mechanisms, for example, as an adaptive response to a post-TBI energy crisis (Bergsneider et al., 2000; Giza & Hovda, 2014; Kawamata et al., 1992; McGuire et al., 2019; Shiga et al., 2006; Sullivan et al., 1998; Vallez Garcia et al., 2016; Xiong et al., 1997).

This reduction of the problem to one of edge allocation naturally leads to an efficient implementation through a closed-form solution, making use of the rearrangement inequality (Hardy, Littlewood, & Pólya, 1952). Given two equally sized sets of real numbers x1, x2, . . , xn and y1, y2, …, yn sorted in ascending order, the following inequality holds:

$$[eqn]$$

for every permutation σ of the indices 1, 2, …, n.

It follows that the cost calculation for the equivalent maximum and minimum cost networks described above reduces to sorting the weights of the observed network in descending and ascending order and taking their dot product with the sorted distance vector according to ascending order, yielding the cost of the minimum- and maximum-cost equivalent networks, respectively. Having a closed form, this metric is computationally efficient and reduces to one sorting operation per network. Our implementation makes use of the sort function of the core MATLAB package, which, in turn, employs an implementation of the Quick Sort algorithm (Hoare, 1961). This algorithm has an average performance of $[eqn]$ , pointing to the scalability of TEWC for parcellations involving a large number of brain regions.

Brain organization does not strictly minimize cost (Bullmore & Sporns, 2012), instead striking a balance between efficient organization and information processing capacity. Consequently, cost metrics are often evaluated within the framework of cost-efficiency (Bassett et al., 2009), quantifying this aforementioned equilibrium.

Here, we introduce two metrics of cost-efficiency by combining the newly introduced TEWC with the ΔC and ΔL deviations, yielding segregative and integrative cost-efficiency (SCE and ICE, respectively) as described below:

$$[eqn]$$ $$[eqn]$$

Analogously to how SWP gauges the coexistence of short path lengths with pronounced clustering, these cost-efficiency metrics evaluate the presence of low cost alongside computational capacity. Here, higher computational capacity is quantified as smaller deviation from the null equivalent models encompassed by the SWP framework, namely, the equivalent regular and random networks. This, in turn, allows for a systematic decomposition of cost-efficiency into components associated with segregation and integration.

Evoked Potentials Analysis

Auditory stimuli were 0.1-ms-long broadband clicks, as well as 25-ms-long tone pips at 10, 37.5, and 65 kHz, all at an intensity of 80 dB SPL. Stimuli were repeated 1,000 times with an interstimulus interval of 500 ms using ultrasonic dynamic speakers (Vifa, Avisoft). Electrophysiological data under auditory stimulation were downsampled to 2 kHz, and line noise was removed as described above. Voltage excursions exceeding 400 μV were also removed and data filtered between 3 and 300 Hz using the EEGLAB function pop_eegfiltnew. Finally, data were segmented into epochs and meaned over trials to extract evoked potential waveforms. More details on stimuli and data preprocessing can be found in the Supporting Information.

For the analysis of auditory evoked potentials, the approach and terminology used by Jonak et al. (2020) are adopted. For each electrode, the N1 component was automatically identified using the MATLAB function findpeaks within a predetermined time window (15 to 100 ms poststimulus). Subsequently, the peak amplitude was extracted and stored for group-level analysis.

Connectivity-Informed Whole-Brain Computational Modeling

Simulations using a phase-oscillator model.

Networks of nodes were defined based on the functional connectomes, with nodes corresponding to electrodes (32 nodes for a full dataset) and edge weights set based on functional connectivity (dwPLI) values. Node dynamics were simulated using the following theta model (Ermentrout & Kopell, 1986):

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where θ_j_ is the phase of node j, Ij(t) is the input current received by node j at time t, $[eqn]$ is the excitability of node j, ξ^(j)^(t) are noisy inputs to node j from distant regions, α_ij_ is the coupling strength between nodes i and j, K is a global scaling constant, N is the number of nodes, and $[eqn]$ is the stable phase of the oscillators in the absence of coupling and noise at which the values of θ were initialized for all nodes. Simulations were run for a total of 6 million timesteps and repeated five times, while the values of the excitability parameter I0 were varied.

Values for K and I0 were selected similarly to Tait et al. (2021). Specifically, I0 was varied from −1.7 to −0.5 in 40 equally spaced values across simulations for a functional connectome. Connectivity matrices were scaled to [0,1] by division with the maximum value in a static fashion, that is, based on the whole recording and were kept fixed for each simulation for that subject while I0 is varied. The specific parameter values used are detailed in Table 1.

Code accompanying the publication by Tait et al. (2021) available from GitHub (https://github.com/lukewtait/AlzheimersBNI.git) was used for simulations.

Outcome measures of ictogenicity.

Following the simulations, two summary measures of the dynamics were calculated. Globally, BNI is defined as the integral of the seizure probability over the excitability parameter values (Lopes et al., 2017):

$$[eqn]$$

where the seizure probability is defined as the proportion of simulation time during which the system is engaged in epileptic dynamics:

$$[eqn]$$

At the node level, NI quantifies the contribution of each node to the whole network’s BNI:

$$[eqn]$$

where BNIpre and BNIpost are the BNI of the full network and the BNI following the removal of node i, respectively. Importantly, NI stands as a normalized metric, ensuring that potential group effects are not spuriously driven by differences in global connectivity strength.

Statistical Analysis

Permutation-based two-way analysis of variance (ANOVA) tests (Anderson, 2001; Anderson & Braak, 2003) were conducted, including terms for the main effect of group, main effect of timepoint, and an interaction effect. This procedure is robust to deviations from normality and outlier values but is not unaffected by violations of homoscedasticity. It can be argued, however, that variance changes in the outcome measures of an experimental group (particularly the injury group) is indeed an effect of interest and sensitivity to it is desirable.

In the case of a significant interaction effect, post hoc pairwise comparisons were performed for the contrasts of interest (blast at 1 month vs. blast at 3 months, blast at 1 month vs. sham at 1 month, blast at 3 months vs. sham at 3 months), using permutation-based two-sample t tests (MATLAB function ttest2 with unequal variances from the statistics and machine learning toolbox), from which the confidence interval for the difference of population means was also extracted.

Correlations between continuous variables were assessed using Spearman correlation (MATLAB function corr). Ninety-five percent confidence intervals were calculated by bootstrapping (10,000 samples). In the presence of outlier observations, a robust regression approach was followed (MATLAB function fitlm with default parameters).

For tests examining group pairs for post hoc pairwise analyses, correction for multiple comparisons was performed using the Bonferroni correction method (VanderWeele & Mathur, 2019). When considering values within a parameterized space (such as the richness or excitability parameters), the Benjamini-Hochberg false discovery rate method was used (Benjamini & Hochberg, 1995). The significance threshold was set at 0.05 for adjusted p values in all tests.

In all figures, “#” denotes a main effect of injury group, while “%” denotes an interaction effect. “*” denotes a significant effect in two-sample t tests. The number of symbols has the usual interpretation (one symbol: p < 0.05, two symbols: p < 0.01, three symbols: p < 0.001).

Acknowledgments

We acknowledge Gil Costa for the use of the scientific figure from SciDraw (DOI: doi.org/10.5281/zenodo.3926343), which was adapted for this work. We also thank SciDraw for providing open-access scientific illustrations that enhance scientific communication. Additionally, we acknowledge Servier Medical Art (licensed under CC BY 4.0) for their medical illustrations.

Author Contributions

Konstantinos Tsikonofilos: Conceptualization; Data curation; Formal analysis; Investigation; Methodology; Software; Visualization; Writing – original draft; Writing – review & editing. Michael Bruyns-Haylett: Conceptualization; Data curation; Investigation; Methodology; Project administration; Supervision; Writing – review & editing. Hazel G. May: Investigation; Writing – review & editing. Cornelius K. Donat: Investigation; Writing – review & editing. Andriy S. Kozlov: Conceptualization; Funding acquisition; Methodology; Project administration; Supervision; Writing – review & editing.

Funding Information

Andriy Kozlov, Royal British Legion (https://dx.doi.org/10.13039/100008631). Andriy Kozlov, Army Research Office (https://dx.doi.org/10.13039/100000183).

Supplementary Material