Neural oscillations track recovery of consciousness in acute traumatic brain injury patients

Abstract Electroencephalography (EEG), easily deployed at the bedside, is an attractive modality for deriving quantitative biomarkers of prognosis and differential diagnosis in severe brain injury and disorders of consciousness (DOC). Prior work by Schiff has identified four dynamic regimes of progressive recovery of consciousness defined by the presence or absence of thalamically‐driven EEG oscillations. These four predefined categories (ABCD model) relate, on a theoretical level, to thalamocortical integrity and, on an empirical level, to behavioral outcome in patients with cardiac arrest coma etiologies. However, whether this theory‐based stratification of patients might be useful as a diagnostic biomarker in DOC and measurably linked to thalamocortical dysfunction remains unknown. In this work, we relate the reemergence of thalamically‐driven EEG oscillations to behavioral recovery from traumatic brain injury (TBI) in a cohort of N = 38 acute patients with moderate‐to‐severe TBI and an average of 1 week of EEG recorded per patient. We analyzed an average of 3.4 hr of EEG per patient, sampled to coincide with 30‐min periods of maximal behavioral arousal. Our work tests and supports the ABCD model, showing that it outperforms a data‐driven clustering approach and may perform equally well compared to a more parsimonious categorization. Additionally, in a subset of patients (N = 11), we correlated EEG findings with functional magnetic resonance imaging (fMRI) connectivity between nodes in the mesocircuit—which has been theoretically implicated by Schiff in DOC—and report a trend‐level relationship that warrants further investigation in larger studies.

. Biomarkers of recovery are greatly needed in DOC to inform prognosis and differential diagnosis, as diagnostic error rates in DOC are as high as 40%, largely owing to the challenges posed by behavioral assessments of level of consciousness (Andrews, Murphy, Munday, & Littlewood, 1996;Childs, Mercer, & Childs, 1993;Schnakers et al., 2009). Diagnostic biomarkers are needed in DOC to identify instances of covert consciousness, that is, consciousness that occurs in the absence of behavioral responsiveness (Huang et al., 2018;Monti et al., 2010). Furthermore, patients' prognoses are also frequently inaccurate, with many families often advised to consider withdraw of care despite over two thirds of patients recovering consciousness when DOC results from traumatic brain injury (TBI) (Giacino et al., 2014;Peberdy et al., 2003;Turgeon et al., 2011). However, predicting which patients are likely to recover consciousness in the absence of prognostic biomarkers remains challenging, thus underscoring the need for such biomarkers (Provencio et al., 2020).
EEG is an attractive modality for biomarkers of postinjury recovery. As a direct readout of cortical activity, EEG is inexpensive, portable, and easily deployed at the bedside. In particular, EEG may be well suited to test hypotheses concerning the role of functional reafferentation of the cortex during coma recovery, that is, restoration of thalamocortical integrity. One such hypothesis, the mesocircuit model, theorizes that the globus pallidus interna (GPi) is disinhibited following diffuse brain injury, and thus, silences the central thalamus, resulting in functional deafferentation of the cortex (Schiff, 2010(Schiff, , 2016. Thus, recovery from DOC requires restoration of striatal functioning and thalamocortical integrity. The latter may be inferred from EEG, given that cortical oscillations such as theta and alpha are thought to be driven by the thalamus (Hughes & Crunelli, 2005;Lindgren et al., 1999;Liu et al., 2012;Sarnthein & Jeanmonod, 2007; Sarnthein, Morel, Von Stein, & Jeanmonod, 2003;Schreckenberger et al., 2004). As such, the loss and recovery of thalamocortical integrity is visible in noninvasive recordings and motivates the "ABCD" model by Schiff. Based on the assumption that specific cortical oscillations indicate varying levels of thalamocortical integrity, Schiff (2016) has defined four dynamic regimes that build on the mesocircuit model, each detectable with EEG and corresponding to a thalamocortical state that indicates progressive circuit recovery. In particular, this model emphasizes thalamic projections to frontal cortical areas, given the privileged role of central thalamic nuclei in anterior forebrain arousal (Schiff, 2020). These EEG types, labeled A-D (hence, ABCD model) are summarized in Table 1. Later types (C, D) denote more progressive recovery (i.e., are "better") than earlier types, (A, B), which correspond to a quiescent thalamic state. Specifically, A-type EEG spectra (featuring no or only low frequency oscillations) are thought to indicate complete cortical deafferentation on a circuit level and a vegetative state on a behavioral level, whereas B-type spectra (featuring theta oscillations) indicate severe deafferentation and a minimally conscious state.
Next, C-type spectra (featuring theta and beta oscillations) may occur when thalamic nuclei fire in burst mode, corresponding to less severe deafferentation and emergence from the minimally conscious state.
Finally, D-type spectra (featuring alpha and beta oscillations) indicate an approximately normal EEG, corresponding to tonic firing of thalamic nuclei and a normal capacity for wakeful consciousness. During the progression from A to D, excitatory synaptic background activity, as well as metabolic rates in cortical, central thalamic, and pallidal tissues, are increasingly restored to normal levels (Comanducci et al., 2020).
Two recent studies have tested the mesocircuit hypothesis in the context of the ABCD model using EEG. Forgacs et al. (2017) found that EEGs from 44 patients who had lost consciousness following cardiac arrest displayed a progression of EEG patterns consistent with the ABCD model. In particular, EEG patterns indicative of greater circuit-level recovery correlated with better outcomes at hospital discharge. More recently, Alkhachroum et al. (2020) used the ABCD model to examine recovery in DOC patients with largely anoxic etiologies treated with amantadine. The best recorded ABCD type increased (A-D) linearly with the percentage of patients who recovered the ability to follow commands.
The foregoing studies offer first evidence for the ABCD model. However, it is unknown how EEG dynamic regimes relate empirically to recovery of the mesocircuit in DOC and whether these EEG types can be used as biomarkers. Deploying the ABCD model as a clinical biomarker will depend crucially on whether the different types can be detected at an acute stage and related to recovery of the mesocircuit.
With the exception of three patients (Alkhachroum et al., 2020), the T A B L E 1 EEG types, criteria, and descriptive statistics Note: Theoretical meanings are based on the original ABCD model by (Schiff, 2016). Some EEG observations (10.6%) could not be classified according to the ABCD model due to peak combinations that did not fit any ABCD type (e.g., a peak in beta without an accompanying peak in theta or alpha). ABCD model has never been applied to patients with severe TBI in the acute stage.

| EEG classification
To classify EEG observations according to ABCD type, we splineinterpolated channel-level power spectral densities (PSDs) to achieve 100 frequency bins per octave. Next, we identified local maxima in PSDs [MATLAB function: findpeaks, min width = 0.1 log 2 (Hz), min prominence = 0.001 log 10 (μV 2 /Hz)] and determined the ABCD type using criteria in Table 1 unclassifiable. We also classified EEG observations according to whether they contained a theta and/or alpha peak (henceforth: θα type) using the mode across channels, given that the long oscillatory periods (>83 ms) of these EEG rhythms are compatible with physiological conduction delays between thalamus and cortex (Swadlow & Waxman, 2012; i.e., thalamocortical communication occurs fully within the excitable phase of the oscillation or one half of the oscillatory period [Fries, 2005]) and are thus especially valuable for inferring thalamocortical integrity.

| fMRI data processing
Patients in our cohort are prone to high motion so we have implemented several preprocessing measures and exclusion criteria to F I G U R E 1 EEG types and behavioral trajectories. (a) Thirteen EEG channels common to all patients were imported for analysis. Actual channel positions varied from patient to patient to accommodate bone flaps and injuries. (b) EEG observations (green circles) were sampled at timepoints corresponding to local maxima of the Glasgow Coma Scale (GCS, black trace), with a 12-hr buffer in between observations. Purple highlights show times with available EEG data. Time is referenced to the patient's earliest GCS score. (c) ABCD type (color) and θα type (

| Brain parcellation using independent component analysis
We performed ICA using the GIFT toolbox (http://mialab.mrn.org/ software/gift/index.html) to parcellate the brain as implemented in Allen et al. (2014). We chose this parcellation approach (Crone et al., 2015;Crone, Lutkenhoff, Vespa, & Monti, 2020) as opposed to standard atlases, because the brains we investigated were severely injured. Thus, it is problematic to assume that parcellation resulting from the average of young and healthy brains adequately represents function in a brain that has been subject to reorganization due to TBI. We defined the cortical regions of interest (ROIs) at an individual level based on the individual functional covariance using groupICA. See Supporting Information Materials and Methods for further details of the parcellation and seed-based connectivity analysis.

| Principal components space clustering
Having classified EEG observations according to the presence or absence of power spectral peaks, we next asked how the foregoing approach would compare with a data-driven approach. To implement such an approach, we began by averaging power (unnormalized) across channels for each patient and log-scaling. Due to the large number of frequency bins (44), we next applied PCA to this feature space and retained the top two PCs according to variance explained.
We then used k-means clustering to identify two clusters in PC space.
Only one EEG observation was used per patient for clustering; EEG observations were selected according to the procedure described below under statistical analysis.

| Statistical analysis
To determine the extent to which EEG spectral peaks predicted patients' behavior, we related acute, longitudinal EEG data to both acute, longitudinal behavioral data (i.e., daily maximum GCS scores) and a single, chronic behavioral datum (i.e., chronic GOSe score). For the former, we used linear mixed models (LMMs) with random intercepts, that is, varying-intercept models. We allowed intercepts but not slopes to vary between patients as we expected patients to have different baselines but did not expect predictors to exert differing levels of influence across patients. Furthermore, given that all models had at least five predictor variables, modeling random slopes for each predictor would result in a cumbersome level of model complexity relative to our sample size. This LMM had the formula where GCS is the daily maximum GCS score, EEG is either an ordinal Note: Our heavily male sample reflects higher risk for TBI in males. "Time to follow up" gives the number of days postinjury when GOSe assessments were performed. Note that GOSe = 1 indicates that the patient was deceased (thus, for Patient 32, the time to follow up was not applicable). Comments give additional details including which analyses, if any, patients were excluded from and why. Abbreviation: N/A, not applicable. of medications. LMMs were fit in MATLAB using the function fitlme.
Instances of unclassifiable ABCD type were treated as missing data.
Patients with fewer than three usable EEG observations were excluded from this analysis due to an insufficient number of observations for inclusion in LMMs.
Additionally, we utilized GCS subscales to infer consciousness in patients and related the same predictors to a binary variable denoting conscious state. Specifically, we inferred the presence of consciousness from a GCS motor score ≥5 or a GCS verbal score ≥4 (see Crone et al., 2020 for further details). We then tested the relationship between predictors and conscious state using generalized LMMs (GLMMs; MATLAB function fitglme) with the logit link function and the formula: where CONSCIOUS is a binary variable denoting the presence or absence of consciousness. For each of the above LMMs and GLMMs, we performed an F-test to evaluate the EEG term.
An alternative model that predicts GCS scores and conscious state based on EEG spectral band power, rather than peak combinations, was also assessed. Separate models were fit for the absolute spectral power (unnormalized) and relative power (normalized by the total 1-45 Hz power). Specifically, we fit LMMs with the formula For the second aim (i.e., relating multiple longitudinal EEG observations per patient to a single chronic outcome), we were unable to fit LMMs or GLMMs because they are not compatible with an unbalanced design featuring dynamic/longitudinal predictors with a static outcome variable. Rather than arbitrarily choosing one EEG timepoint for each patient to create a balanced design, we utilized multiple linear regression models with a resampling approach that randomly sampled one EEG observation per patient with replacement for each of 9,999 resamples and constructed an empirical distribution of test statistics.
Resamples that yielded invalid combinations (i.e., a rank deficient regression design matrix) were discarded and replaced with a new sampling. We reported the results of the resample that yielded the median t-statistic for the variables EEG, DELTA, THETA, ALPHA, or BETA; for this reason, the number of resamples N was chosen as an odd number such that the middle-most test statistic would be defined.
Our initial multiple linear regression model was specified as where GOSe is the chronic GOSe score. For alternative models using spectral power as predictors, we used multiple linear regression models with the formula Next, we compared GCS and GOSe scores between clusters identified using k-means clustering. This clustering approach requires only one EEG observation per patient. Because including all data from each patient in our clustering would have biased clusters toward patients with more observations, we once again utilized a resampling approach with replacement using 9,999 resamples to include only one EEG observation per patient per resample. Invalid resamples were discarded and replaced as described above. For each random resampling, we constructed a feature space using channel-averaged spectral power across all frequency bins, applied PCA, retained the two PCs that accounted for the highest proportion of variance, and identified two clusters in PC space using k-means clustering. For consistent labeling of clusters across resamples, we used one label for all clusters that had a centroid with a PC 1 -coordinate ≥0 and another label for those with a PC 1 -coordinate <0. We then performed multiple linear regression using the model where SCORE is the behavioral measure (GCS, conscious state, or GOSe) and CLUSTER is a binary variable denoting the patient's power space cluster membership. For predicting conscious state (a binary variable), we utilized logistic multiple regression. GCS scores were selected corresponding to the day/time of each patient's randomly resampled EEG observation. The resample yielding the median tstatistic for CLUSTER was then selected for reporting. Note that this resampling procedure was performed separately for GCS and GOSe.
Finally, to relate EEG to fMRI connectivity, we correlated the mean proportion of channels with a θα peak with the z-scored fMRI connectivity for each of four fMRI ROI pairings: thalamus-striatum, striatum-globus pallidus, thalamus-prefrontal cortex (PFC), and thalamus-posterior cingulate cortex (PCC). We used θα type rather than ABCD type due to the greater variance and weaker skew of the former (see Section 3). Furthermore, to reduce the number of data points at floor (i.e., all θαÀ) or ceiling (i.e., all θα+), we choose to examine the proportion of EEG channels with a θα peak, rather than the modal (and thus binary yes/no) θα type of each EEG observation. To derive noise-robust estimates for each patient, we averaged this proportion across all EEG observations within 48 hr of fMRI acquisition.
Correlations were derived using the Pearson coefficient. We did not include covariates in this analysis for two reasons: (a) EEG observations from both before and after fMRI acquisition were included, and so temporal precedence of predictors could not always be established, and (b) because EEG and fMRI are both intimately related measures of brain activity, their relationship is not confounded by the same variables that were covaried for in other analyses. Outliers were identified and excluded from correlation analysis using a threshold of three scaled median absolute deviations from the median of either variable (peak proportion or fMRI connectivity) using the MATLAB function isoutlier with default parameters. To account for multiple testing, we applied false discover rates (FDR, Benjamini-Hochberg) to correct pvalues (Benjamini & Hochberg, 1995).
In all models, we corrected for the four hypothesis tests outlined in the introduction using a Bonferroni correction, yielding a test-wise criterion of α = .0125 (for fMRI connectivity, this was performed in addition to the FDR correction).

| RESULTS
Following preprocessing and quality control, we retained 320 EEG observations across 38 patients (31 male) with ages ranging from 19 to 84 years (40 ± 17 years, mean ± SD). Patient demographics are summarized in Table 2 and Figure S1. The number of usable EEG observations per patient ranged from 1 to 16 (8.4 ± 4.3, mean ± SD), with 1.8-30 min of usable data per observation (24 ± 7.2, mean ± SD).
Thus, despite our modest patient sample size, we analyzed an average of 3.4 hr of data per patient.
EEG observations were most commonly categorized as A-type (71.6%), with those remaining categorized as B-type (15.0%), C-type (1.88%), D-type (0.94%), or unclassifiable (10.6%). Separately, 30.9% of observations exhibited a peak in the theta-alpha band in the majority of channels (θα+). Distributions of GCS scores for each type are described in Table 1 and Figure 2. See Figure 3 for examples of each ABCD type and Figure S3 for behavioral trajectories and EEG variables for all patients.

| Thalamically-driven oscillations track behavioral recovery
Given the very small proportion of EEG observations that were categorized as C-type or D-type, we opted to create a binary variable by grouping types B, C, and D together, thus avoiding outlier effects. To investigate ABCD type in relation to acute behavioral state, we used an LMM with GCS as the dependent variable. We excluded four patients from the model due to having fewer than three observations, plus an additional patient with missing medication data, yielding n = 33. We found that ABCD type (i.e., A vs. B, C, D) significantly predicted GCS (p < 0.0001, F (1,268) = 17.0), with more progressive states (B, C, D) associated with higher GCS scores (t = 4.13); see Table S1 for F-statistics and p-values of the covariates and intercept in this and other models). However, unclassifiable EEG observations, which were omitted from the above model, were significantly more likely to correspond with higher GCS scores than classifiable EEG observations (t = À2.59) as revealed using an LMM with ABCD classifiability (yes/no) as the EEG predictor (p = .010, F (1,296) = 6.71).
ABCD type also significantly predicted the presence/absence of consciousness in patients, as inferred from GCS subscales (p = .0037, Next, we fit the same models substituting θα type for ABCD type. Similar to ABCD type, θα type was significantly predictive of both GCS (p = <.0001, F (1,296) = 16.2) and conscious state (p = <.0001, F (1,296) = 15.6), where θα + patients were more likely to have higher GCS scores (t = 4.03) and to be conscious (t = 3.95). Adding θα type as a predictor to the LMM fit using ABCD type only improved prediction of GCS scores before adjusting our α level for the additional hypotheses outlined in the introduction (p = .032, likelihood ratio stat = 4.59); thus, a larger sample might demonstrate an added value for including both predictors in the model. On the other hand, adding ABCD type as a predictor to the GLMM fit using θα type unambiguously improved prediction of conscious state (p <.0001, likelihood ratio stat = 75.4).
As a benchmark to compare the above models against, we also predicted GCS and conscious state using spectral power in the delta, theta, alpha, and beta frequency bands. Absolute alpha (p = .0037, F (1,293) = 8.56) and beta (p = 2.9 Â 10 À5 , F (1,293) = 18.1) power significantly predicted GCS, with lower alpha power (t = À2.93) and higher beta power (4.25) corresponding to higher GCS scores after accounting for covariates. Note, however, that alpha power was positively related to GCS in raw correlations that did not control for other predictors (absolute power, r = .29; relative power, r = .065). Adding absolute alpha power to LMMs that predicted GCS using ABCD type (p = .19, log likelihood stat = 1.70) or θα type (p = .23, log likelihood stat = 1.41) did not significantly improve model fit in either case.
However, adding absolute beta power to the LMM that predicted GCS using θα type significantly improved model fit (p = .0028, likelihood ratio stat = 8.93), and a trend level improvement was observed when added to the LMM that predicted GCS using ABCD type (p = .019, likelihood ratio stat = 5.52). Using relative power, we again found that alpha (p = .0022, F (1,293) = 9.58) and beta (p = .0011, F (1,293) = 10.9) power significantly related to lower and higher GCS scores, respectively. Model fits were significantly improved when relative alpha power was added to LMMs predicting GCS (ABCD type: p = 3.2 Â 10 À4 ; θα type: p = 5.1 Â 10 À5 ), but relative beta power did not significantly improve model fits (ABCD type: p = .079, log likelihood stat = 3.08; θα type: p = .32, log likelihood stat = 0.98).
Having examined spectral power as a predictor of GCS, we next examined it as a predictor of conscious state using GLMMs. No absolute power features significantly predicted conscious state, though we observed a trend for beta power (p = .058, F (1,293) = 3.63), with higher power corresponding to consciousness (t = 1.91); however, relative delta power did significantly predict conscious state (p = .0069, F (1,293) = 7.41), with lower power corresponding to consciousness (t = À2.72). Nonetheless, adding relative delta power to GLMMs that predicted conscious state from ABCD type or θα type did not improve model fit (the original models lacking relative delta power featured greater maximized log-likelihoods and thus no test was performed).
Next, we used multiple linear regression with resampling to investigate whether ABCD type predicts chronic outcomes as measured with GOSe (again, types B, C, and D were grouped together to avoid outlier effects). Four patients missing GOSe score were excluded and four previously excluded patients with fewer than three EEG observations were reincluded here, see Table 2. The number of postinjury days to follow up did not correlate with GOSe scores (r = .05, p = .76) and thus was not considered to be a confound. For all hypotheses tested, the median t-statistic across resamples stabilized after~1,000 resamples ( Figure S4), that is, well within the number of resamples we performed (N = 9,999). ABCD type did not significantly predict GOSe (p = .37, t = 0.92). We then substituted θα type for ABCD type in the original multiple linear regression model and repeated the resampling procedure. As with ABCD type, θα type did not significantly predict GOSe (p = .52, t = 0.65). Neither absolute nor relative power in any frequency band significantly predicted GOSe.
Finally, we used a data-driven clustering approach to determine whether clusters based on EEG spectral power (unnormalized) would predict GCS (n = 37) and/or GOSe (n = 33). Once again, we utilized N = 9,999 resamples to choose randomly one EEG observation per patient per resample. Two PCs were retained that, for all resamples, explained at least 80% of the variance in spectral power (resamples that did not explain this proportion of variance in the first two PCs were discarded and not counted toward N). K-means clustering was used to identify clusters in PC space (see Figure 4). Each resample yielded one cluster whose centroid had a positive PC 1 -coordinate and one cluster whose centroid had a negative PC 1 -coordinate (Figure 4d), and clusters were therefore labeled accordingly. Cluster labeling did not significantly predict GCS (p = .11, t = 1.64), conscious state (p = .32, t = 1.00), or GOSe (p = .68, t = 0.42). See Tables S2-S4 for p-values, t-statistics, and F-statistics from resampled tests.
3.2 | Relating thalamically-driven oscillations to mesocircuit recovery in a small sample  (Table 3). We therefore correlated the proportion of channels with a θα peak (peak proportion), averaged across EEG observations within 48 hr of scanning, with fMRI connectivity.
F I G U R E 3 Examples of ABCD model types. Power spectral densities (PSDs) from each channel are color-coded according to their ABCD type; unclassifiable channels are colored black. Peaks detected for classification are indicated with circles. Channels with a θα peak are dashed. Shaded areas are colored according to frequency band. (a) A-type EEG corresponding to a GCS score of 3. All channels were categorized as Atype, and peaks were only present in the delta band. (b) B-type EEG corresponding to a GCS score of 9. In total, one channel was categorized as A-type, nine channels as B-type, two channels as C-type, and one channel as unclassifiable (due to the presence of a beta peak without an accompanying theta or alpha peak). Eleven channels showed a θα peak. (c) C-type EEG corresponding to a GCS score of 7. In total, two channels were categorized as C-type, and the remaining nine channels were uncategorizable. Seven channels showed θα peaks. (d) D-type EEG corresponding to a GCS score of 14. In total, one channel was categorized as A-type, one channel as D-type, and the remaining nine channels were uncategorizable. Because D-type is a more progressive type than A-type, the tie between A-type and D-type (one channel each) is broken by D-type. Twelve channels showed a θα peak We found a trend relating peak proportion to fMRI connectivity between thalamus and PFC (r = .68, Pearson coefficient, p = .08, FDR corrected). Trend-level relationships were also observed between peak proportion and thalamo-striatal (r = .63, p = .10, FDR corrected) and striatal-pallidal BOLD signal coupling (r = .59, p = .10, FDR corrected). No relationship was observed between peak proportion and BOLD signal coupling between thalamus and PCC (r = .01, p = .97, FDR corrected). See Figure 5 for scatter plots and correlations. To ensure that these correlations were not overly sensitive to the 48-hr time limit used to select EEG observations, we also performed a sensitivity analysis and computed Pearson coefficients for time windows ranging 6-72 hr in length. Correlations appeared stable and maximized usable data when EEG observations were included within 13-51 hr of fMRI scanning ( Figure S5). Thus, our 48-hr time limit appears to maximize the amount of available data without including EEG observations outside of the observed window of stability. Note: Forty-four EEG observations from 11 patients were included in the correlation of EEG with fMRI connectivity measures. EEG observations were considered if they fell within 48 hr of MRI. Type proportion gives the proportion of EEG channels with a type progressed beyond the A-type (B, C, or D). Peak proportion gives the proportion of EEG channels with a θα peak. EEG time post-fMRI gives the number of hours before that the EEG observation took place after the patient's MRI scan (negative values indicate that the EEG observation preceded the MRI scan). Given the more favorable statistics of peak proportion (less skew and greater variance), this measure was correlated with fMRI connectivity ( Figure 5) after averaging across all available observations. also predict these variables, only relative alpha power significantly improves model fit when added as a predictor alongside ABCD type.

| DISCUSSION
Given its ability to predict conscious state, the ABCD model may have applications as a diagnostic biomarker, for example, to detect instances of covert consciousness (Huang et al., 2018;Monti et al., 2010) in patients lacking behavioral responsiveness. However, our findings suggest that, in acute patients, more progressive types (C and D) are rare, as patients may generally achieve this level of recovery after leaving the ICU. Accordingly, a more parsimonious categorization of EEG type based on the presence or absence of θα (4-12 Hz) peaks may be equally useful, as it concentrates specifically on thalamically-entrained cortical rhythms that form the core of the ABCD model (see Table 1). Furthermore, this approach has several practical advantages over the ABCD model: (a) all usable data are classifiable, (b) one classification will capture the majority of channels for all EEGs with an odd number of channels, and (c) because all data are classifiable, this approach does not introduce the sampling bias that occurs using the ABCD model. Finally, our study is the first to offer provisional evidence of the relationship between EEG and mesocircuit recovery, as measured with fMRI connectivity. Importantly, as predicted by the mesocircuit hypothesis, which emphasizes central thalamic projections to frontal cortex (Schiff, 2016), we found a stronger relationship between θα peak proportion and coupling between thalamus and PFC, versus coupling between thalamus and PCC.

| EEG oscillations as a readout of mesocircuit recovery
The mesocircuit model explains postinjury forebrain dysfunction in terms of inactive striatal medium spiny neurons (MSNs; Schiff, 2010Schiff, , 2016). MSN firing rates may be particularly sensitive to diffuse injury, as they require high levels of both dopaminergic neuromodulation and spontaneous background corticostriatal and thalamostriatal synaptic F I G U R E 5 Correlations of EEG (proportion of channels with θα peak or "peak proportion") with fMRI connectivity. Peak proportion was averaged across all usable EEG observations within 48 hr of MRI. Data points are sized proportionally to the number of averaged EEG observations and colored according to the mean hours elapsed between EEG and MRI (absolute value, warmer colors indicate greater mean time elapsed). Outliers (>3 scaled median absolute deviations from the median of either variable) are indicted with red circles around data points and were excluded from analysis. Dotted black lines represent the least-squares linear fit for included data points. We corrected for testing across four region of interest (ROI) pairings using false discovery rates (FDR). Correlations are reported as Pearson coefficients. (a) Peak proportion versus z-scored thalamo-striatal connectivity: r = .63, p = .10, FDR corrected (uncorrected: p = .07). Note the presence of two excluded outliers. (b) Peak proportion versus z-scored striatal-pallidal connectivity: r = .59, p = .10, FDR corrected (uncorrected: p = .07). Note the presence of one excluded outlier. (c) Peak proportion versus z-scored thalamo-prefrontal cortical (PFC) connectivity: r = .68, p = .08, FDR corrected (uncorrected: p = .02). (d) Peak proportion versus z-scored thalamo-posterior cingulate cortical (PCC) connectivity: r = .01, p = .97, FDR corrected (uncorrected: p = .97). Note the presence of two excluded outliers input to reach firing threshold (Grillner, Hellgren, Menard, Saitoh, & Wikström, 2005). When these necessary conditions are disrupted by diffuse, multifocal brain-injury, GABAergic MSNs go offline, and the GPi is released from striatal inhibition (Schiff, 2010). As a result, GPi powerfully inhibits central thalamus, leading to functional deafferentation of cortex ( Figure S6). Furthermore, the striatum receives weaker cortical and thalamic drive as a result of GPi inhibiting thalamus, leading to a positive feedback loop in which GPi is further disinhibited and central thalamus further inhibited.
As originally proposed by Schiff, the ABCD model is a readout of mesocircuit recovery, with thalamically-driven EEG oscillations indicating progressive levels of thalamocortical integrity (Schiff, 2016).
Consistent with this hypothesis, we found that ABCD type predicts conscious state, suggesting clinical application as a diagnostic biomarker in DOC. Because frontal cortical areas receive denser projections from central thalamus than other cortical areas (Deschenes, Bourassa, & Parent, 1996;Morel, Liu, Wannier, Jeanmonod, & Rouiller, 2005), frontal cortex is proposed to be disproportionately affected by mesocircuit dysfunction (Schiff, 2010(Schiff, , 2016. Accordingly, we found that the proportion of EEG channels displaying a θα peak is correlated with thalamo-prefrontal connectivity (r = .68), but not with thalamo-posterior cingulate connectivity (r = .01), in a small subsample of patients. Although neither relationship was significant after correcting for multiple comparisons, the former correlation appears promising and therefore warrants further investigation in larger samples. Furthermore, we also observed provisional, trend-level evidence relating EEG peak proportion to subcortical mesocircuit connections (thalamo-striatal and striatal-pallidal connectivity). Taken together, our findings in this small subset of patients are a promising, albeit inconclusive, indicator that EEG oscillations in the 4-12 Hz range may be reflective of mesocircuit recovery.

| Improving the ABCD model
While alternative models that predicted acute variables based on spectral power-used as benchmarks for the ABCD model-found that both absolute and relative alpha and beta power predicted GCS and that relative delta power predicted conscious state, only relative alpha power significantly improved the fit of any model with ABCD type as a predictor, as determined using log likelihood ratio tests. Adding absolute beta power, rather than relative alpha power, to the same model resulted in a trend-level improvement, but significantly improved model fit in the same LMM with θα type substituted for ABCD type. These results suggest that incorporating both information regarding peaks in frequency bands and the area under the curve in one or more frequency bands might improve the utility of the ABCD model as a diagnostic biomarker in DOC. Surprisingly, we found that even though the presence of alpha peaks in the ABCD model was associated with higher GCS scores, alpha power itself was associated with lower GCS scores in our models, though only after accounting for beta power and other covariates, as the raw correlation between alpha power and GCS was positive for both absolute and relative power when not controlling for other predictors. This multifaceted relationship between alpha oscillations and behavioral recovery underscores the potential importance of considering both oscillatory power and peaks.
Furthermore, as the resonant frequency of alpha oscillations changes with development and aging (Chiang, Rennie, Robinson, Van Albada, & Kerr, 2011;Donoghue et al., 2020), the area under the curve in the alpha frequency range might be more sensitive to alpha activity than the presence or absence of a local maximum in the alpha range. For example, consider a slow 7 Hz posterior alpha rhythm peak that would still "leak" energy into the Finally, given that delta oscillations are widely regarded as indicators of cortical down states and unconsciousness (Buzsaki, 2006;Koch, Massimini, Boly, & Tononi, 2016;Massimini, Ferrarelli, Sarasso, & Tononi, 2012), one might be surprised by our finding that relative delta power, while predictive of conscious state, did not improve the fit of any peak-based model, especially given that these models did not already consider any peaks <4 Hz. While delta oscillations are clearly seen in states of unconsciousness including the slow wave sleep (Brown, Lydic, & Schiff, 2010;Franks, 2008;Murphy et al., 2011), anesthesia (Franks, 2008;Murphy et al., 2011;Purdon et al., 2013;Supp, Siegel, Hipp, & Engel, 2011), and DOC (Hussain et al., 2019;Kaplan, 2004;Sutter & Kaplan, 2012) including coma and the vegetative state (Sutter & Kaplan, 2012), high amplitude delta oscillations can nonetheless be observed in a variety of circumstances in which individuals are fully conscious and responsive, such as Angelman syndrome (Frohlich et al., 2020); for a comprehensive review of this and other cases, see Frohlich, Toker, and Monti (2021). Based on these other findings and our results herein, the consideration of parameters in the delta band may not be necessary for improving the ABCD model. Although the ABCD model outperformed a data-driven clustering approach based on PCA and K-means clustering, it is possible that further refinement of the clustering model or a more sophisticated unsupervised learning approach could yield competitive results. Furthermore, while the clustering approach assigned all EEG observations to a cluster, the ABCD model was unable to classify some EEG observations, resulting in a sampling bias revealed in our study using an LMM, where EEG observations corresponding to higher GCS scores were less likely to be classified. This is likely due to the fact that in conscious states, there are a greater number of "illegal" peak combinations that cannot be classified (e.g., a beta peak without any accompanying theta or alpha peak). This bias can be addressed by replacing ABCD type with θα type, which allows all EEG observations to be classified.
We also observed that dynamic regimes from the ABCD model reflecting high levels of thalamocortical integrity (C and D types) were rare in acute patients. This is perhaps unsurprising, as more recovered patients would exhibit C and D types, but they are typically moved from the ICU to less intense care and, thus, not part of our acute sample. We must therefore consider whether circuit-level recovery precedes behavioral recovery early enough to be useful for predicting outcomes, as an ideal biomarker should predict a patient's pending recovery before the patient is moved out of the ICU. It is possible that the mere reemergence of theta oscillations (B-type), prior to that of other oscillations, already heralds recovery. However, we did not observe a predictive relationship between ABCD type (after grouping types B, C, and D together to avoid outlier effects) and chronic outcomes as measured with GOSe. This might be due to the differing sensitivities of the GCS and GOSe assessments, where the former is better suited for detecting the presence/absence of consciousness, whereas the latter is better suited for detecting recovery of daily living skills. Because ABCD types of B or higher were rarely observed in conscious states, ABCD type relates very closely to GCS but less well to GOSe. It is possible that had our sample included more instances of C-types and D-types, which might relate to recovery of daily living skills, a relationship with GOSe would have also been detected.
Finally, the mean GCS score corresponding to EEG observations that were categorized as C-type or D-type did not exceed that of those categorized as B-type (Table 1, Figure 2c), suggesting that Ctype and D-type do not indicate progressive recovery beyond B-type.
However, given the small number of EEG observations that were categorized as C-type or D-type, it is possible that we did not have sufficient data to observe the progressive hierarchy predicted by the ABCD model. GCS ranges for A-type (GCS = 3-15) and B-type (GCS = 4-15) were both broad, suggesting low specificity. A-type was observed at ceiling (GCS = 15), though conversely, B-type was never observed at floor (GCS = 3). Given these data, EEGs that have progressed beyond A-type likely indicate that patients have already begun recovering behavioral responsiveness.
Besides relating EEG to behavioral data, an important aim of our study was to validate an automated procedure for categorizing EEG observations according to ABCD type using quantitative objective criteria. We found that our automated EEG classification corresponded significantly with contemporaneous behavior, suggesting that the ABCD model can be implemented computationally, without the need for manual scoring. However, we found that 10.6% of EEG observations could not be classified according to the ABCD model due to peak combination that are not defined by the model. Additionally, since each EEG channel had five possible classifications (A, B, C, D, or unclassifiable), one classification did not capture the majority of EEG channels in all instances. The foregoing issues are solved using a more parsimonious classification between two types, θαand θα+, with one type capturing the majority of EEG channels in all instances.
We found that classification based on θα type performed equally well as ABCD classification in predicting GCS (both were strongly predictive) and GOSe (neither was predictive). The main advantage of the ABCD model was in improving prediction of conscious state when added to the GLMM with θα type as a predictor. Our findings suggest that the ABCD model performs similarly to a parsimonious and computationally simpler categorization scheme that infers thalamocortical integrity based on the presence of theta and/or alpha EEG oscillations in acute TBI patients.

| Limitations and future directions
Our study has a number of methodological limitations: (a) Because EEG data were collected in a clinical setting, EEG channel placement was variable and the number of channels common to all patients was low. This precluded the use of EEG source localization or spatial filtering (e.g., surface Laplacian) in our analysis. Furthermore, although patients with more progressive EEG types (C and D) were largely absent from our study due to the fact that EEG data were only collected in the ICU, we believe that our findings have greater translatability since the need for diagnostic and prognostic biomarkers is most prevalent in the acute stage. (b) Despite collecting a large amount of EEG data (multiple days) from each patient, our patient sample size was relatively small (38 patients analyzed total, 33 patients in the main analysis), and so caution should be used in generalizing these findings to larger patient populations. (c) Patients were administered a large number of different medications, the influences of which could not be entirely controlled in our study. However, by utilizing logistic PCA, we covaried for two PCs that explained roughly two thirds of the variance in medication data. (d) The GOSe, our measure of chronic (~6 month) outcome, was conducted in some cases as a phone interview and thus an indirect assessment. The indirect nature of the assessment may have added substantial noise that could reduce correlations with EEG measures, possibly explaining our absence of findings when relating EEG to GOSe. (e) Our k-means clustering approach to predicting GCS did not perfectly mirror that based on EEG spectral peak classification (ABCD type or θα type), as the latter approach utilized LMMs and thus discarded four patients with insufficient longitudinal data who were included in the former approach.
However, we believe that our results unambiguously show no relationship between data-driven clusters and GCS; thus, the asymmetry in our two analyses is likely inconsequential, though we acknowledge that a more sophisticated or better developed clustering or unsupervised learning approach could perhaps yield competitive results. (f) Our correlation of EEG with fMRI measures was limited by the fact that many patients were not scanned until leaving the ICU, at which point EEG was no longer acquired. Thus, only 11 patients had at least one EEG observation within 48 hr of fMRI. We were therefore likely underpowered to detect significant relationships between neural oscillations (EEG) and mesocircuit integrity (fMRI). However, after removing outliers that would easily influence our small sample, we observed statistical trends.
Two future directions are strongly encouraged by our results.
Firstly, given the paucity of C and D types found in our data from acute patients in the ICU, future work should apply our automated peak fitting approach to data from patients moved to less intensive care, for whom more progressive types (C and D) may be observed.
Such an investigation might be more successful in relating ABCD type to long-term outcome, given the potentially larger spread in ABCD type after sampling more recovered patients. This sample might also help determine if a more parsimonious approach based on θα type is indeed preferable over the ABCD model, even when patients are more likely to exhibit a large spread in ABCD types. Secondly, we advocate for further work relating the proportion of EEG channels displaying θα peaks to mesocircuit functional connectivity in larger samples. Given the large correlation we observed between EEG peak proportion and thalamo-prefrontal connectivity (r = .68), only three additional patients (n = 14) are needed to achieve >80% statistical power. Thus, even moderately larger samples may be beneficial in detecting significant correlations. Additionally, dynamic causal modeling (Friston, Harrison, & Penny, 2003) applied to larger datasets may be useful for inferring the directionality of interactions between mesocircuit nodes, for example, to show that EEG peak proportion correlates with increased thalamic drive to cortex or decreased pallidal inhibition of central thalamus.

| CONCLUSIONS
Noninvasive readouts of mesocircuit recovery are desirable for diagnosing DOC in TBI patients. Because circuit-level recovery should precede behavioral recovery, such a readout may also be more useful than behavioral scales (e.g., GCS) in predicting which patients will recover consciousness. Our findings show that the reemergence of neural oscillations tracks recovery in acute TBI patients, as relevant EEG measures correspond strongly and significantly with behavioral responsiveness and conscious state in the ICU while also exhibiting a trend of greater spatial extent (i.e., peak proportion) with increasing thalamocortical connectivity. These results suggest that thalamicallydriven EEG oscillations may serve as diagnostic biomarkers in TBI and DOC, although the failure to detect a relationship between EEG and chronic outcome in our patient cohort may be only due to the low occurrences of C and D types in our sample, and thus, may still have prognostic value when investigated at a less acute timepoint.