QUANTITATIVE ASSESSMENT OF ELECTROENCEPHALOGRAM REACTIVITY IN COMATOSE PATIENTS ON EXTRACORPOREAL MEMBRANE OXYGENATION

2.1. Subjects and study protocol

This study was undertaken with the approval of the Institutional Review Board of Johns Hopkins University. Subjects were enrolled in a prospective observational study of ECMO patients in whom continuous EEG was conducted as part of standard clinical care¹⁷. This study was a multidisciplinary endeavor coordinating efforts from the cardiac intensive care, cardiovascular surgical intensive care, medical intensive, and neurocritical care units, and the neurophysiology department. Adult ECMO patients (age >18 years) with Glasgow Coma Scale score (GCS) of 3–8 were enrolled from November 2017 to April 2020. As part of routine multidisciplinary clinical care, all patients were followed by a neurocritical care physician from time of ECMO cannulation until hospital discharge or death. Baseline neurological examination of each patient was assessed on day one of ECMO cannulation. Cessation or reduction of sedative medications (“sedation holidays”) were routinely given for each patient as deemed safe by the treating physician(s). Sedation holidays were usually undertaken between days three to five of ECMO support and EEG was placed to allow an EEG assessment while sedation was minimized or absent. GCS was collected at time of routine clinical exams when off sedation. Continuous EEG was performed as a clinical standard of care for patients with DOC and collected in a clinical cEEG database. EEG data were de-identified and transferred for further analysis.

2.2. EEG data collection

A montage of 19 electrodes, placed according to the international 10–20 system, of continuous EEG was used to record at a sampling frequency of 200 Hz. Studies were performed using Nihon Kohden digital EEG systems (Nihon Kohden, Tokyo, Japan). The collected EEG signals were de-identified and stored in the European Data Format (EDF). VAR and qEEG analysis were then conducted independently.

2.3. Stimulation trial

The stimulation protocol was conducted by an experienced clinical team of fellows and EEG technicians as part of clinical care. Three types of stimulation were administered to patients: auditory, peripheral, and nostril. A set of stimulation trials was performed in the same order each time (with increasing intensity of stimulation):

• Auditory Stimulation (AUD): with increasing loudness, calling of the patient’s name followed by clapping near their ear

• Peripheral Stimulation (PER): pinching the finger and toe nailbeds bilaterally

• Nostril Stimulation (NOS): inserting a cotton swab deep into the nasal cavity

Each individual stimulation trial was performed at least once and was followed by a rest period to allow for separation of the effects from each individual stimulation. Each subject had at least one set of stimulation trials using all three stimuli, and some subjects had more than one. Those patients that had more than one set of stimulation trials had at least one day in between each set. As EEG features and patterns may change due to the dynamic nature of acute critical illness^18,19, each stimulation trial was classified as a separate EEG sample, and samples from the same patient were treated as independent for purposes of analysis.

When stimulation was performed, annotations and timing marks were added to the EEG signal to synchronize the onset of stimulation and recorded EEG. Subjects’ responses were observed and noted by the clinical team.

2.4. Visual analysis of EEG-R

All EEG data were reviewed retrospectively by one neurologist (S.-M.C.) and one neurophysiology specialist (E.K.R.). EEG data were de-identified during VAR, such that clinicians were blinded to which of the patients’ EEG they were examining, nor were they aware of the survival or current state of the patient whose EEG they were examining. Each sample of EEG data was classified, following standard clinical practice of consensus discussion, with “presence of visual reactivity” (R+) or “absence of visual reactivity” (R−) based on previously established criteria²³. R+ indicates that a subject demonstrated EEG-R in response to at least one type of stimulation, while R− indicates that a patient did not show a response to any stimulation. To prevent biasing of patient outcomes, determination from visual analysis by these clinicians was not used to guide patient treatment, nor were the clinicians informed of the progression of the patient’s condition.

2.5. Artifact detection and data preprocessing

Data were processed in MATLAB (version 2020a on macOS Big Sur, MathWorks, Natick, MA) using both scripts developed by the authors and scripts included in the publicly available software package EEGLAB²⁰.

A finite impulse response band-pass filter (0.1 to 50 Hz) was used to remove ambient electrical noise. An independent component analysis was performed to identify common artifacts, such as muscle activity, eye movements, or body movements, and data containing any of these artifacts were excluded from analysis. To remove eye movement (electro-ocular, EOG) and muscle and motion (electro-muscular, EMG) artifacts, the EEG data were processed by a second-order blind identification algorithm from the Automatic Artifact Removal (AAR) toolbox in EEGLAB^21,22.

An epoch consisting of two consecutive ten-second, periods of EEG recording corresponding to the time before and after a stimulation was extracted from each trial. Each ten-second period was additionally divided into five two-second periods (sub-epochs) for further processing. Each stimulation period lasted from 3–5 seconds; these were not included in the record. Figure 1 depicts the stimulation protocol and data processing methodology. Following EEG data preprocessing, the spectral power and permutation entropy of the frequency bands [δ (0.5–4 Hz), θ (4–8 Hz), α (8–13 Hz), β (13–30 Hz), γ (30–50 Hz)] were computed, resulting in ten metrics.

2.6. Quantitative EEG measures

2.7. Computation of reactivity score

To reduce the potential effect of noise from and the dimensionality of each individual channel, each metric (the spectral power or permutation entropy in a given frequency band) was aggregated across all 19 EEG channels; the default aggregation method used was to calculate the mean of the channels, however other methods were also applied (see Algorithm Permutations). The computation of reactivity score was then based on a standard score (SS):

where SS_metric,post is the standard score of metric in a particular post-stimulation sub-epoch post;; metric_post is the aggregate value of metric (either spectral power or permutation entropy) across all EEG channels at time post; mean_pre and std_pre are the mean and standard deviation of the aggregate value of metric over the pre-stimulation period. This standardizes the time series of each metric in the post-stimulation period against the pre-stimulation period. Thus, a positive SS_metric,post in the post-stimulation period indicates an increase from the pre-stimulation level while a negative SS_metric,post indicates a decrease. The value of a metric was considered to be indicative of the presence of EEG-R if it had an absolute normalized value higher than the threshold of 1.96, chosen as the approximate value of the 97.5 percentile point of the standard normal distribution. Alternative threshold values for this criterion were also tested (see Algorithm Permutations).

In order to aggregate the values of SS across the five-time segments of the post-stimulation period, the mean value of the difference between the criterion and any metric in any time segment that met the criterion were computed (see Algorithm Permutations). The resulting Reactivity Score (R-score) then becomes:

where SS_metric,post is any standard score in the time series of a metric that meets the criterion for EEG-R; R_metric is the R-score for a specific metric (i.e. sub-R-score). This results in ten sub-R-scores - one for each of the metrics calculated for the sub-epoch. The R-score for an epoch is then the sum of all of the sub-R-scores.

2.8. Algorithm permutations

A subset of parameters of the algorithm used to calculate the R-scores from the EEG were identified for analysis by permutation. Parameter permutations were included based solely on mathematical validity (rather than underlying physiology) to limit bias regarding potential neurological mechanisms of effectiveness:

1. Baseline high end limit of frequencies

2. Aggregation across multiple EEG channels

3. Threshold of standard deviations from baseline

4. Aggregation of peak values post-stimulation

5. Aggregation of multiple SS readings

All combinations of these permutations of the algorithm were used to calculate R-scores that were fed into the modeling stage. The goal of this permutational approach was to blindly assess which combination provided the best approximation of the conclusions drawn by the expert neurologists, and thus reverse engineer through algorithmic means the methods being applied by clinicians.

2.9. Baseline comparison

An EEG epoch during which no stimulation was performed (baseline) was taken from each patient and the reactivity scores for each of these epochs were computed. The baseline scores were divided into two cohorts: those judged by the expert clinicians to contain signs of reactivity (R+) and those that were not (R−). The groups were then compared by a Cramer–von Mises two-sample test at α = 0.01 under the null hypothesis that the baseline EEG-R scores for the groups were the same. Verifying that the baseline values are of equivalent underlying distribution is key to validating the algorithm.

2.10. Statistical modeling

Statistical analyses were performed using MATLAB (version 2021a on 64-bit Windows 10, MathWorks, Natick, MA). In order to avoid bias, different investigators performed the reactivity score algorithm permutations (Y.Z.) and the statistical modeling (A.W.). Details of permutations used were anonymized prior to being transmitted between investigators and thus the statistical analysis was blinded to which permutations were used in each set of reactivity scores.

In order to convert from the metric space defined by the algorithm for scoring the presence of reactivity in the filtered EEG into a probability space that could be used for prediction of patient outcome in a graded and intuitive manner, each of the five sub-bands of spectral power and permutation entropy, as well as the total values thereof, were analyzed in a paradigm of univariate and bivariate stepwise logistic regression models. Models were generated for four binary outcome classes: one for each class of stimulation (AUD, PER, NOS), and one in which all stimulation classes were merged (MER). Models were provided a set of twelve potential predictors (five sub-band R-scores for spectral power, one total R-score for spectral power, five sub-band R-scores for permutation, one total R-score for permutation entropy) and one response variable (the VAR decision). Models were considered successful if both the overall model and all coefficients fulfilled a P < 0.10 threshold. Generated models were then pruned using χ2 goodness of fit test against a baseline (predictor-less) model; deviance value comparison for selection between nested or branched models; and magnitude and significance of penultimate step coefficients to limit overfitting. In order to account for the restricted sample size and correct for potential population bias in our use of a sample of convenience, a King-Zeng Prior Correction (KZPC) factor²⁶ was implemented during the model generation process, selecting τ = 0.48 based on the clinical population occurrence rate of patient EEG-R reported by two recent multi-center studies^11,12. This logistic regression model produces the target graded probability value of the clinician judgement for the presence of EEG-R for each epoch. In order to provide a baseline measure of differentiation to illustrate the potential capabilities of the approach, the output of the models were dichotomized with a 50% threshold. These dichotomized example models were ranked by sensitivity and specificity against the outcome criterion and only those models performing at greater than fifty percent on both were kept.

3.1. Baseline comparison

The cohort sizes for the Cramer-von Mises two-sample test were N_Baseline,R+=5 and N_{Baseline,R−}=21, with a computed test statistic of 0.1848. Given the critical value for these cohort sizes at α = 0.01 is 0.7322, we thus fail to reject the null hypothesis that the baseline values are similarly distributed. This indicates that the algorithm is identifying qualities particular to the post-stimulation period rather than qualities intrinsic to the EEG of the reactive and unreactive cohorts.

3.2. Stimulation Prediction Outcomes

Of the 20, 87, 22, and 129 epochs analyzed for AUD, PER, NOS, and MER, 13 (65%), 74 (85%), 16 (73%), and 103 (80%) were from patients rated R− on visual analysis of the EEG.

Due to the large number of permutations tested on the algorithm used to generate the reactivity, a similarly large number of statistical models needed to be generated to cover all of the scores. Those example models that met the earlier stated threshold for p-value and pruning criteria were then ranked, separately and equally weighted, on sensitivity (percent of true positives correctly predicted) and specificity (percent of true negatives correctly predicted); positive and negative in this experiment correspond to whether the clinicians identified the patient as R+ or R−. As the outcome in this case was binary, both sensitivity and specificity needed to be higher than 50% to outperform chance. Specifically, 8400 regression models were generated for each stimulation type; the vast majority of these did not reach statistical significance or were otherwise pruned. Out of the total 8400 potential models generated for each of AUD, PER, NOS, and MER stimulation, only 65 (0.77%), 272 (3.2%), 55 (0.65%), and 994 (11.8%), respectively, were judged as acceptable by the inclusion criteria.

For AUD example models (Table 2.A), nine achieved a sensitivity higher than 50%, and six of these achieved a specificity of 100%; all but one of them identified the spectral power in the beta band as the sole predictor, while the last identified total spectral power.

For PER example models (Table 2.B.i), six achieved a sensitivity of 100% with a specificity higher than 90%. All but one of these models identified the total value of spectral power as its predictor, while the last identified the permutation entropy of the beta band. Conversely, seven models (Table 2.B.ii) achieved 100% specificity with a sensitivity of 80%; all but one of which identified the total value of spectral power, while the last identified the permutation entropy in the alpha band.

For NOS example models (Table 2.C), four achieved a sensitivity higher than 50%, with three of these achieving a specificity of 100%. Two of them identified the permutation entropy in the gamma band as the sole predictor, while the other identified total spectral power. While the gamma band has the most potential for contamination from EMG artifacts, we are confident that our noise and artifact reduction protocol removed unrelated activity from the gamma band based on the existing evidence of these techniques in other applications^21,22.The R-scores for reactive and non-reactive patients differed substantially in the gamma band, many of the latter showing minimal gamma band activity at all.

Finally, for MER example models (Table 2.D.i), seven models performed well on both measures, achieving higher than 75% sensitivity and higher than 90% specificity. Two of these models used the combination of the total spectral power and total permutation entropy, with a stronger weight placed on the latter, while the other five used solely spectral power in the gamma band. Alternatively, four models (Table 2.D.ii) achieved higher than 84% on both sensitivity and specificity, however, each model identified different predictors.

3.3. Successful Permutations

When evaluating the breakdown of permutations by examining their prevalence in distinct categories of stimulation example models (Figure 2), one can clearly see that the majority of these used the full range of frequencies, as is to be expected given the prevalence of beta and gamma frequency ranges and total sums in the identified predictors. There is a similar pattern for the AUD and NOS categories on the standard deviation threshold permutation; these permutations involve changing the initial threshold (1.96) at which a post-stimulation baseline deviation is recognized as displaying reactivity to be more (1.64) or less (2.58) permissive. Both AUD and NOS evaluation showed increased accuracy when assigned the latter permutation, indicating their epochs to have a lower signal-to-noise ratio in comparison to the PER epochs, which would lead to higher false positive deviations. Conversely., the PER and merged categories included more permissive models with lower thresholds.

Regarding aggregation of EEG channels, models for AUD and PER stimulation found measures of dispersion (StdDev, IQR, RMS) to be more useful than measures of central tendency (Mean, Median) or position (Max, Upper Quartile) in determining presence of reactivity; further, all such models found a positive association between the predictor and the response. As broader dispersion indicates greater differences between input factors, this suggests that reactivity for these two types of stimulation may suggest a differential effect on different brain regions. Conversely, example models for NOS used measures of central tendency, indicating an overall strong response from all regions. The merged category showed, as intuition would suggest, a balance of both groups. None of the example models found the maximum value of the EEG channels to be useful in predicting decision of reactivity, providing support for the need to examine multiple electrode channels.

Regarding aggregation when examining post-stimulation peaks, which used a two-factor matrix for its permutations (mean vs max, first peak vs all peaks), there was a more mixed response. The example models for the merged class identified the mean of the first peak as the predominantly useful method in making prediction, whereas the AUD class did not use this in any of its example models, and the PER class was split between the mean of the first peak and the mean of all peaks. The AUD and PER classes found similar levels of utility in the application of extrema, to either set of peaks.

Finally, for aggregation of multiple successful readings, the PER category predominantly applies the mean of the readings and thus provides a more stable reading, whereas the AUD and MER classes found the measures of extrema to be more accurate. This is, however, expected given that the PER category had more trials.