Brain functional connectivity during the first day of coma reflects long-term outcome.

OBJECTIVE: In patients with disorders of consciousness (DOC), properties of functional brain networks at rest are informative of the degree of consciousness impairment and of long-term outcome. Here we investigate whether connectivity differences between patients with favorable and unfavorable outcome are already present within 24 h of coma onset.
METHODS: We prospectively recorded 63-channel electroencephalography (EEG) at rest during the first day of coma after cardiac arrest. We analyzed 98 adults, of whom 57 survived beyond unresponsive wakefulness. Functional connectivity was estimated by computing the 'debiased weighted phase lag index' over epochs of five seconds duration. We evaluated the network's topological features, including clustering coefficient, path length, modularity and participation coefficient and computed their variance over time. Finally, we estimated the predictive value of these topological features for patients' outcomes by splitting the patient sample in training and test datasets.
RESULTS: Group-level analysis revealed lower clustering coefficient, higher modularity and path length variance in patients with favorable compared to those with unfavorable outcomes (p < 0.01). Within all features, the path length variance in the network provided the best positive predictive value (PPV) for favorable outcome and specificity for unfavorable outcome in the test dataset (PPV: 0.83, p < 0.01; specificity: 0.86, p < 0.01) with above-chance negative predictive value and accuracy. Of note, the exclusion of patients with epileptiform activity (20 in total) eliminates all false positive predictions (n = 6) for path length variance.
INTERPRETATION: Topological features of functional connectivity differ as a function of long-term outcome in patients on the first day of coma. These differences are not interpretable in terms of consciousness levels as all patients were in a deep unconscious state. The time variance of path length is informative of comatose patients' outcome, as patients with favorable outcome exhibit a richer repertoire of path length than those with unfavorable outcomes.

Introduction

Converging evidence on brain functional connectivity in humans has highlighted the organizational principles underlying normal cognition during conscious wakefulness (Avena-Koenigsberger et al., 2017). Functional studies describe this brain network architecture as a modular structure organized around few – highly interconnected – hubs playing a central role in the information transfer across the whole network (Meunier et al., 2010). During loss of consciousness, these patterns of functional connectivity typically break down into loosely connected units of weakened long-range spatial and temporal connections in accordance with current influential theories of consciousness (Mashour and Hudetz, 2018). On these bases, prognostic and diagnostic markers for disorders of consciousness (DOC) patients have benefited from analyzing electrophysiological and functional magnetic resonance imaging (fMRI) data from a network perspective. In DOC patients, several weeks to months after injury, functional connectivity features stratify patients as a function of clinical severity (Chennu et al., 2017, Engemann et al., 2018), provide predictive information for long-term outcome (Chennu et al., 2017, Demertzi et al., 2019) and show complementarity with other metrics of consciousness evaluation (Chennu et al., 2017, Engemann et al., 2018). During early coma, a state of unarousable unresponsiveness with absence of awareness of both self and the environment, neuroimaging and electrophysiological studies of functional connectivity patterns and their relation to patients’ outcome are sparser compared to DOC patients. A seminal functional neuroimaging study uncovered a radical reorganization of the anatomical distribution of main ‘hubs’ in functional networks of comatose patients in comparison to healthy controls and irrespective of long-term outcome (Achard et al., 2012). In another study, comatose patients with poor outcome exhibited a disruption of the connectivity underlying the so-called ‘default mode network’ (Norton et al., 2012), a set of regions encompassing posterior cingulate cortex, precuneus, medial prefrontal cortex, and bilateral temporoparietal junctions. More recently, the functional connectivity strength between posterior cingulate cortex and medial prefrontal cortex has been identified as a discriminative feature between comatose patients who went on to recover (favorable outcome – FO) and those with unfavorable outcomes (UO) three months after brain injury (Malagurski et al., 2017).

This body of functional neuroimaging evidence encourages further investigation of the predictive power of the functional connectivity patterns during coma based on electrophysiological recordings which up to now have been only reported in cohorts with few patients and low density montages (Beudel et al., 2014). EEG can be used at bedside, and – due to its superior temporal resolution – presents complementary investigative power in comparison to fMRI. In this study, we focus in particular on the first day of coma for two reasons: from the clinical perspective, the availability of quantitative predictive markers at this early stage can complement current neurological tests carried out during the first hours of coma for which the overall predictive performance remains limited (Rossetti et al., 2016), especially for favorable outcome. Understanding whether discriminative connectivity features between patients with good or bad prognosis arise already within hours of brain injury – while patients remain deeply unconscious – can help to identify precursors of consciousness recovery at early stages after brain injury.

With these two motivations, we recorded resting 63-channel EEG during the first 24 h of coma after cardiac arrest (CA) while patients were mildly sedated and treated with targeted temperature management. In each recording, we derived functional connectivity matrices based on established synchronization measures across the EEG electrodes (Vinck et al., 2011). In order to derive global metrics of network organization, we extracted their topological features, which were previously shown to co-vary with consciousness levels in disorders of consciousness patients (Lehembre et al., 2012, Chennu et al., 2014). In addition, we estimated the temporal variability of these topological measures within each recording in order to investigate the richness in the repertoire of functional network configurations. Based on both the topological features and on their temporal variability we aimed at: first, evaluating whether functional network features vary with outcome at group level; second, estimating the performance of these properties in predicting patients’ outcome during the first day after coma onset.

Materials and methods

Post-anoxic comatose patients

Electroencephalography was prospectively recorded from 138 comatose patients after CA on the first and – when possible – on the second day of coma at the intensive care units of the University Hospitals of Lausanne (n = 69), Sion (n = 4), Fribourg (n = 2) and Bern (n = 63) between July 2014 and January 2018. The study was approved by the ethics committees of the respective institutions. Prior to any recording, informed written consent was signed by a family member, legal representative, or treating clinician not involved in this study. Patients in the study routinely received targeted temperature treatment at either 33 °C (n = 18), or 36 °C (n = 115), or no intervention (spontaneous temperature; n = 5) during the first 24 h of coma. Ice packs, cooling pads or intravenous ice-cold fluids with a feedback controlled cooling device (Arctic Sun System, medivance, Louisville or thermogard XP, ZOLL Medical, Zug, Switzerland) were used for temperature control. Midazolam (0.1 mg/kg/h), propofol (2–3 mg/kg/h) and fentanyl (1.5 μg/kg/h) were administered for analgesia-sedation and neuromuscular blocking agents vecuronium, rocuronium, or atracurium were given in order to control shivering when needed. After 24 h the administration of sedation, muscle relaxants and temperature treatment was weaned off. All recordings included in the present study were acquired during the first day of coma.

Patients functional and neurological state was assessed using the Full Outline of UnResponsiveness (FOUR (Wijdicks et al., 2005)), repeated at least twice within the 72 h following CA.

Individual patient’s outcome was defined according to the best functional outcome within three months post CA. We considered both the Cerebral Performance Categories (CPC (Booth et al., 2004)), a semi-structured phone interview at three months after CA, and the neurological state during hospitalization. The CPC characterizes the patient’s capacity to conduct independent daily activities at work and at home. Full recovery and conscious with moderate disability correspond to CPC of 1 and 2 respectively; a CPC 3 indicates conscious with severe disability; CPC 4 corresponds to coma or persistent vegetative state, and CPC 5 to death. Patients scoring a CPC of 1 to 2 at any point in time during the first three months after cardiac arrest were considered to have favorable outcomes, while those who never awoke and ultimately died within this time span were considered to have unfavorable outcomes (our cohort contained no patients with CPC 4 at three months, i.e. no patients remained in a vegetative state at three months).

Clinical assessments

A certified neurologist assessed pupillary, oculocephalic, corneal reflexes and motor reactivity to pain stimulation after our EEG recordings on the second and third day of coma. Based on bedside EEG recordings, background reactivity was evaluated by certified neurologists (Tsetsou et al., 2015). Bilateral median nerve somatosensory evoked potentials (SSEP) evaluation was carried out within 24–48 h after CA. Withdrawal of care was decided based on a multidisciplinary approach, requiring at least two of the following criteria (Rossetti et al., 2010): a) unreactive EEG background after targeted temperature treatment and off sedation, b) treatment-resistant epileptiform EEG and/or myoclonus, c) bilateral absence of N20 in SSEP, and d) incomplete return of brainstem reflexes. One or two clinical EEG recordings were assessed, the first 6–24 h after CA, the second at 36–48 h after CA after withdrawal of targeted temperature treatment. In some cases, recording EEG was only possible on the first day of coma due to medical interventions, patients awakening or life support being withdrawn before a second recording could be performed.

EEG methods

Acquisition and preprocessing

We recorded continuous resting EEG for 8 to 20 min at 1200 Hz using a 63 active ring electrode array (g.HIamp, g.tec medical engineering, Graz, Austria) in the 10–10 system, referencing to the right ear lobe. Additional electrodes were attached to the patient’s chest to record electrocardiography (ECG). Impedances of all active electrodes were kept below 50 kΩ and all data were recorded with online 0.1–100 Hz band-pass and 48–52 Hz notch filters. Patients also underwent an auditory protocol after completion of resting EEG recordings (Tzovara et al., 2013).

The continuous EEG data were first segmented into epochs of five seconds, downsampled to 500 Hz and filtered using a windowed finite-impulse band-stop filter for ± 5 Hz at 50 Hz, 100 Hz, and 150 Hz to minimize line noise.

Based on visual inspection, we discarded epochs and electrodes containing artifacts. Patients with high amounts of muscle noise (e.g. due to shivering) and machine noise were excluded from further analysis.

We carried out an independent component analysis (Jung et al., 2001) to remove components containing identifiable ECG artifacts. In a next step, the discarded electrodes were interpolated based on a weighted average of neighboring channels before re-referencing the EEG to the average offline reference.

EEG data analysis

Functional connectivity

We estimated undirected connectivity from the EEG electrodes in FieldTrip (Oostenveld et al., 2011) using a debiased version of the weighted phase lag index (Vinck et al., 2011) (hereafter referred to as dwPLI) at 10 Hz with a slepian tapering sequence resulting in 1 Hz smoothing, following previously observed spectral differences in partially overlapping data, centered around the alpha band (Kustermann et al., 2019).

For each recording, based on the absolute magnitude of the connections, 50% to 90% of connectivity matrices were pruned in steps of 2.5%. We binarized the resulting matrices to account for inter-subject differences based on overall connectivity strength, thus easing results interpretation. We computed the topological features from the pruned and binarized connectivity matrices.

In order to estimate the variability of topological features over time, we first extracted the dwPLI over non-overlapping sliding windows of 5 epochs each. Binarization and pruning as described above were applied for each connectivity matrix estimation. Based on this set of connectivity matrices, we assessed the variability of topological features over time, within each recording, by computing the time variance of the topological properties extracted from each matrix.

Topological features

We extracted topological features of the functional connectivity using the brain connectivity toolbox (BCT 23) based on the dwPLI at 10 Hz. Graph visualization was based on Gephi (Bastian et al., 2009). We estimated the clustering coefficient, characteristic path length, modularity and participation coefficient. These measures are derived from graph theory (Bullmore and Sporns, 2009), where electrodes and synchronization level represent nodes and connections of the graph, respectively. All topological features were computed based on the connectivity matrix from the whole, epoched continuous resting EEG recording and – separately - on connectivity matrices based on the 5 consecutive epochs (corresponding to 25 s) for the time-variance.

The clustering coefficient reflects the degree of clustering within a given network, as it computes the degree of connectivity among nodes that connect to a given node (Watts and Strogatz, 1998). Overall clustering is then estimated by averaging its values across nodes.

The characteristic path length (Watts and Strogatz, 1998) (hereafter referred to as ‘path length’) is a measure of distance from a given node to another node. In a binarized matrix this will reflect the number of steps from one electrode to another via interconnected nodes. Averaging yields an aggregate measure of the distance between any two given nodes in the network.

The modularity (Girvan and Newman, 2002, Newman, 2006) reflects the degree of sub-division of the network into highly connected 'sub-communities’ (or groups) with limited size. This is computed by maximizing within group connections and minimizing between-group connections, using the Louvain algorithm (Fortunato, 2010). The visualization of modularity (Fig. 2) is carried out via the Yifan-Hu algorithm (Hu, 2005).

Fig. 2

Exemplary graphs of one patient with FO and one patient with UO at 82.5% pruning. Colors reflect modules as identified by the Louvain algorithm 29. Spatial configuration of sub-communities as calculated through the Yifan-Hu algorithm 30. Each inset represents the projection of the community affiliation on the EEG electrode array. One can note that communities are more spatially segregated in the patient with FO than in the patient with UO reflecting higher modularity.

The participation coefficient (Guimera and Amaral, 2005) computes the degree to which nodes within a sub-community retain connections to other sub-communities. Thus, a high participation coefficient indicates a high degree of diversity – of connections within and outside the sub-community.

Group-level statistical analysis

Group-level statistical analyses of topological features and of their time-variances were based on their summed values across pruning levels. This summation has the advantage of avoiding the multiple comparison problem that would stem from comparing topological values at each pruning level. In addition, sum-based differences arise from general topological properties of the network irrespective of specific pruning level and are thus easier to interpret. We conducted independent sample t-tests for all summed topological feature and considered significant any comparison yielding p < 0.013, i.e. p < 0.05 corrected for multiple comparisons.

EEG-based outcome prediction

The second main analysis of this study was the prediction of patient outcome based on EEG topological features of 10 Hz dwPLI and their temporal variability. Preprocessing, connectivity matrices and topological features extraction were identical to those described for the group-level analysis. All available recordings were split in training and test datasets using about 50% of the data in each split. We used the summed values across pruning levels for outcome prediction. The training dataset was used to derive the area under the Receiver Operating Characteristic and for optimizing the threshold of each topological value in order to maximize the positive predictive value (PPV; number of correctly predicted favorable outcomes relative to all predicted favorable outcomes) with the constraint of predicting at least ten patients with FO. The reason for optimizing the prediction for favorable outcome is that these results could be used in the clinical routine to avoid inappropriate end-of-life decisions and are less affected by the ‘self-fulfilling prophecy’ which confounds the interpretability of unfavorable outcome predictions. The early prediction of favorable outcome is mostly valuable in clinical contexts with high prevalence of withdrawal of life sustaining treatment, a rate that can vary widely across countries (Mashour and Hudetz, 2018).

We used the test dataset for an unbiased selection of the best predictive performance in an independent dataset. We evaluated its significance by estimating confidence intervals based on a permutation distribution (n = 100) of outcome labels of the training dataset where for each permutation we estimated an optimal threshold as to maximize the PPV, and this new threshold was used to quantify the prediction performance on the test dataset. Prediction results were considered significant if different from the permuted values based on a Wilcoxon signed-rank test (p < 0.001, Table 1).

Table 1

Prediction results based on the time-variance of clustering coefficient, path length, modularity and participation coefficient in training and test sets derived from the functional connectivity matrices on the first day of coma. Statistically significant predictions are highlighted in bold. PPV = Positive Predictive Value; NPV = Negative Predictive Value.

Condition	PPV	NPV	Sensitivity	Specificity	Accuracy
Clustering coefficient
Training	0.71	0.47	0.37	0.79	0.54
Test	0.69	0.46	0.30	0.82	0.52
Path length
Training	0.79	0.50	0.41	0.84	0.59
Test	0.83	0.56	0.50	0.86	0.65
Modularity
Training	0.65	0.83	0.96	0.26	0.67
Test	0.67	1.00	1.00	0.32	0.71
Participation coefficient
Training	0.67	0.46	0.44	0.68	0.54
Test	0.68	0.48	0.43	0.73	0.56

Results

Preprocessing and included patients

After preprocessing, because of diffuse artifacts (either muscle or machine noise) we had to exclude EEG recordings of 40 patients, 23 of which had FO, (Supplemental Material 1 for an exemplar case). The total number of patients included in the analysis was therefore 98, 57 had FO.

Within the included patients, the final number of epochs for FO and UO after preprocessing did not differ significantly (mean and standard deviation indicated throughout text, x̅ = 176.68 ± 5.03, x̅ = 174.61 ± 5.48, t(96) = 0.28, p = 0.78; removed artifacts amounted to 8.93 ± 12.70, UO 9.83 ± 19.27 trials for FO and UO, respectively; t(96) = −0.28, p = 0.78). The number of electrodes interpolated did not differ between FO and UO (x̅ = 0.79 ± 0.20, x̅ = 1.41 ± 0.34, t(96) = −1.66, p = 0.10). As a result of the ICA analysis, the number of removed components did not significantly differ between FO and UO patients (number of removed ICA components: 0.84 ± 0.82, 1.15 ± 1.13 for FO and UO respectively, t(96) = -1.55, p = 0.13).

Within the included patients, the average FOUR score of FO was 4.84 ± 2.31 whereas UO scored 2.39 ± 1.96. In all patients for which individual sub-scores could be retrieved (34 in total), the eye and motor responses were lower than 2 and 4 respectively, indicating a deep unconscious state on the first day.

Group level analyses of topological features

Clustering coefficient differed significantly by outcome (t(96) = −4.29, p < 0.01; Fig. 1 and Fig. 2) as well as modularity (t(96) = 2.98, p < 0.01). There were no significant differences observed between FO and UO for the path length (t(96) = 1.99, p = 0.05) and for the participation coefficient (t(96) = -1.11, p = 0.27).

Fig. 1

Group-level analysis of the 10 Hz dwPLI derived topological features for FO (n = 57 in blue) and UO (n = 41 in red). The left column shows the values as a function of percentage of connections pruned. The right column shows the probability densities of the summed values and their individual scores (lower half of each subplot). Boxplots reflect the median and the interquartile range (IQR) via the box and the 1.5 IQR via the whiskers. Displayed significant p-values result from t-test comparison between FO and UO, corrected for multiple comparisons. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Variability of topological features

For the time-variance of the clustering coefficient sum there was no difference between FO and UO (t(96) = 0.44, p = 0.665, Fig. 3). Patients with FO exhibited higher path length time-variance compared to patients with UO, t(96) = 3.67, p < 0.001. Neither modularity time-variance (t(96) = 1.81, p = 0.073) nor the participation coefficient time-variance (t(96) = 0.83, p = 0.411) differed between patients with FO and UO (Fig. 4).

Fig. 3

Group-level analysis of the 10 Hz dwPLI derived time-variance of the topological features for FO (n = 57 in blue) and UO (n = 41 in red). The left column shows the values as a function of percentage of connections pruned. The right column shows probability densities of summed values and their individual scores (lower half of each subplot). Boxplots reflect the median and IQR via the box and the 1.5 IQR via the whiskers. Displayed significant p-values result from t-test comparison between FO and UO, corrected for multiple comparisons. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4

Outcome prediction for FO (n = 57) and UO (n = 41) based on the time-variance of the topological features. Each dot represents the sum of topological features across pruning levels for an individual patient. Blue and red dots refer to patients included in the training and test datasets respectively. The horizontal dashed line indicates the classification cut-off for outcome prediction. Values above the threshold are predictive to have FO. Classification cut-offs were chosen to maximize the positive predictive value in the training dataset. The selected parameters were then validated using an independent test set. * = significant at p < 0.01 level. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

In all the topological features and their time-variance, we evaluated the effect of the removed artifacts by comparing the results before and after artifacts removal. Results did not significantly differ at group-level. At single-subject level we observed an over estimation at high pruning level of the majority of the considered metrics when considering data without removing artifacts.

Prediction analysis

Topological features

We estimated the area under the Receiver Operating Characteristics (AUC) values for all topological features as well as the prediction results after optimizing for the threshold in order to maximize the PPV. For clustering coefficient the AUC was 0.77, modularity 0.69, path length 0.72 and for the participation coefficient 0.59.

We found significant PPV and prediction results in the test dataset for all metrics but the participation coefficient, although none of these results provided very high scores (<=0.70, Supplemental Materials 2 and 3).

The AUC for the variability of the topological features were 0.64 for clustering coefficient, 0.50 for modularity, 0.56 for participation coefficient and 0.68 for path length. Hereafter, we present threshold specific predictions tuned for PPV in the training set. Using path length time-variance we obtained the highest PPV in the test set (Table 1). We correctly predicted 15 of 18 patients to have FO (PPV: 0.83, z = 8.68, p < 0.001), whereby all false positives had epileptiform activity (n = 6; Table 2). NPV - the fraction of correctly predicted UO to patients predicted as UO - was 0.48. Sensitivity was 0.50 - the ratio of correctly identified FO relative to all patients with FO - with 15 out of 30 correct FO predictions and specificity (0.86, z = 8.69, p < 0.001) - the ratio of correctly identified UO relative to all patients with UO - were significant with 19 out of 22 UO patients correctly predicted. Overall, accuracy was significant as well (0.62, p < 0.001) with 29 out of 52 patients correctly predicted.

Table 2

Demographic and clinical comparisons of patients above and below the path length time-variance threshold, separately for FO and UO; the number of patients for which values were not available are shown in parentheses. Statistical results are shown based on t-statistic or Fisher’s exact test. ROSC = Return of spontaneous circulation; FOUR = Full Outline of UnResponsiveness; SSEP = Somatosensory evoked potentials.

	Path length time-variance
	Above threshold	Below threshold	p-value
Favorable outcome
N	26	31
Time to ROSC (min), M ± SD	19.11 ± 14.28	22.25 ± 13.66	0.41
Age (y), M ± SEM	57.76 ± 15.87	65.90 ± 13.70	0.05
Gender (male), N (Nmissing)	22 (0)	26(0)	1
Time to 1st EEG (h), M ± SEM	18.85 ± 6.92	21.76 ± 5.34	0.10
Time from 1st to 2nd EEG (h), M ± SD	23.84 ± 3.02	22.51 ± 2.77	0.16
Temperature (°C), M ± SEM	35.61 ± 0.85	35.64 ± 1.11	0.92
Propofol (mg/kg/h), M ± SD (n)	2.46 ± 1.53 (21)	2.54 ± 0.74 (23)	0.85
Midazolam (mg/kg/h), M ± SD (n)	0.13 ± 0.05 (12)	0.11 ± 0.05 (8)	0.46
Fentanyl (μg/kh/h), M ± SD (n)	0.89 ± 0.46 (20)	1.55 ± 1.49 (22)	0.07
FOUR score	4.83 ± 2.76	4.85 ± 1.79	0.99
Eye	0 ± 0	0 ± 0
Motor	0.86 ± 1.36	0.20 ± 0.60	0.30
Brainstem	2.57 ± 0.90	3.30 ± 1.10	0.18
Respiratory	0.57 ± 0.49	0.60 ± 0.49	0.91
EEG reactivity 1st day, N (Nmissing)	5 (6)	5 (3)	0.72
EEG discontinuity 1st day, N (Nmissing)	4 (1)	14 (1)	0.02
Electrographic epileptic activity 1st day, N(Nmissing)	0 (1)	1 (1)	1
Pulmonary etiology of CA, N (Nmissing)	6 (0)	5 (0)	0.52
Pupillary reflexes, N (Nmissing)	20 (6)	25 (3)	0.25
Corneal reflexes, N (Nmissing)	18 (6)	23 (3)	0.68
Motor response, N (Nmissing)	14 (6)	22 (3)	0.52
SSEP absence, N (Nmissing)	0 (14)	0 (15)
Unfavorable outcome
N	6	35
ROSC (min), M ± SD	34.16 ± 14.55	28.22 ± 13.78	0.42
Age (y), M ± SEM	59.66 ± 10.27	66.00 ± 13.42	0.25
Gender (male), N (Nmissing)	4 (0)	23 (0)	1
Time to 1st EEG (h), M ± SEM	22.83 ± 5.84	20.07 ± 7.03	0.37
Time from 1st to 2nd EEG (h), M ± SD	23.80 ± 1.32	24.08 ± 4.99	0.81
Temperature (°C), M ± SEM	36.22 ± 0.51	35.82 ± 0.91	0.18
Propofol (mg/kg/h), M ± SD (n)	3.08 ± 1.024 (6)	2.18 ± 1.41 (13)	0.20
Midazolam (mg/kg/h), M ± SD (n)	0.03 ± 0.02 (2)	0.09 ± 0.05 (16)	0.04
Fentanyl (μg/kh/h), M ± SD (n)	0.50 ± 0.24 (3)	1.06 ± 0.73 (17)	0.05
FOUR	2.40 ± 1.62	2.38 ± 2.04	0.98
Eye	0 ± 0	0.08 ± 0.27	0.34
Motor	0 ± 0	0.46 ± 1.08	0.17
Brainstem	2.00 ± 1.41	2.15 ± 1.50	0.87
Respiratory	0.50 ± 0.50	0.62 ± 0.49	0.74
EEG reactivity 1st day, N (Nmissing)	5 (1)	25 (7)	1
EEG discontinuity 1st day, N (Nmissing)	5 (0)	26 (3)	1
Electrographic epileptic activity 1st day, N (Nmissing)	6 (0)	13 (3)	0.02
Pulmonary etiology of CA, N (Nmissing)	2 (0)	5 (3)	0.30
Pupillary reflexes, N (Nmissing)	3 (1)	18 (10)	0.62
Corneal reflexes, N (Nmissing)	3 (1)	9 (10)	0.36
Motor response, N (Nmissing)	1 (1)	2 (10)	0.43
SSEP absence, N (Nmissing)	3 (3)	9 (18)	0.24

In the test dataset, highest accuracy was provided by the time-variance of modularity (0.71, z = 8.65, p < 0.001) with 37 correct predictions out of 52. PPV correctly predicted 30 out of 45 FOs (0.67, z = 6.94, p < 0.001). NPV was at 1.00, reflecting correct identification of UO in 7 out of 7 patients (z = 8.68, p < 0.001). Sensitivity was maximal, correctly identifying all patients with FO (1.00, z = 8.69, p < 0.001), while specificity was non-significant (0.32, n.s.).

Clinical characteristics

We compared the clinical characteristics of all patients falling above and below the threshold for the time-variance of each of the topological features and the correlation between each clinical characteristic and topological features (Supplemental Material 5) in order to exclude that clinical differences could explain the prediction results. We performed the comparisons using t-tests or Fisher’s exact tests. Here, we present the results based on path length time-variance (Table 2 and Supplemental Material 4 for the comparison based on all the other time-variance metrics). Epileptiform activity was present in all patients with UO with high path length (n = 6).

Discussion

We investigated the connectivity properties of resting state EEG in comatose patients after CA as a function of their outcome in a large dataset collected across four hospital sites on the first day of coma. Functional connectivity analyses highlighted differences in the topological organization of brain networks between FO and UO in the alpha band at 10 Hz. Favorable outcome was associated with higher modularity and lower clustering than poor outcome. These features of functional connectivity in patients with FO are consistent with a network characterized by densely connected sub-communities at mesoscopic scale as captured by the modularity (and less densely connected than UO at local scale as indicated by the clustering coefficient), yet integrated at global level by few sparse links (Fig. 2). The group-level analyses of the topological features time-variance yielded significant differences for the path length, that is to say the average communication distance was more variable in FO than UO patients. In terms of long-term outcome prediction, both the topological metrics and their time-variance produced several significant results (Table 1, Supplemental Material 2, 3). Within the best prediction results, the path length variance had a PPV of 0.79 and 0.83 for FO and specificity of 0.84 and 0.86 for UO, in training and test datasets, respectively. By excluding patients with epileptiform activity (n = 6), we would have no false positive predictions (Table 2), thus maximizing PPV similarly to a previous study (Tzovara et al., 2016). These predictions are comparable with previously reported prediction based on the progression of auditory discrimination between the first two days of coma (Tzovara et al., 2016) and consistent with previously observed prediction based on oscillatory power at 10 Hz for long term coma outcome in a partially overlapping cohort (Kustermann et al., 2019).

Previous literature about EEG functional connectivity in comatose patients

Previous attempts to characterize the functional networks in comatose patients were either based on functional magnetic resonance imaging (Achard et al., 2012, Malagurski et al., 2017, Forgacs et al., 2017) or on low-density EEG (Beudel et al., 2014). In the latter study, functional networks were derived from adjacency matrices of shared frequencies across the whole frequency spectrum, a connectivity metric that is influenced by volume conduction and thus prevents its interpretation in terms of communication between distinct neural regions. Despite these differences, our functional connectivity analysis based on non-instantaneous phase coupling across 63-channels EEG around 10 Hz yielded results consistent with this prior EEG work: as patients with FO exhibited lower clustering and low path length was specific of patients with UO (see Supplemental Material 2 for the specificity in UO patients). As suggested in this previous paper, these functional characteristics are compatible with the disruption of long range and metabolically demanding connections as a consequence of severe anoxia in patients with UO. This disruption will likely affect patients to a variable degree depending on the extent of the affected brain region and on the duration of anoxia. Our results suggest that the degree of this disruption can be well captured by the topological properties of functional networks by identifying early in time which patients will be able to reverse these impairments and eventually survive. To the best of our knowledge, our study is the first to investigate the temporal variability of functional connectivity during early coma using neurophysiologically interpretable connectivity metrics and to report significant predictive performance for long-term patients’ outcome.

Functional connectivity as a correlate of consciousness

Our analyses were inspired by a long-standing literature investigating the neural correlates of consciousness during rest (Engemann et al., 2018, Boveroux et al., 2010, Deco et al., 2014). In particular, modularity and characteristic path length are among the possible measures capturing the segregation and integration of functional brain networks (Rubinov and Sporns, 2010), characteristic properties of neural activity underlying conscious experience in humans (Tononi and Edelman, 1998). In our study, the modularity was higher in patients with FO than UO, suggesting that the functional segregation at the mesoscopic level is already present in patients with FO before regaining consciousness. On the other hand, the characteristic path length, related to the global degree of communication among close and distance modules, was not significantly different between FO and UO, yet a low value of path length was specific in patients with UO, thus suggesting that high level of integration was specific within UO patients (Rubinov and Sporns, 2010). We should also note that other EEG studies quantified the segregation through the clustering coefficient and the integration through the participation coefficient (Rizkallah et al., 2019, Lee et al., 2019), which in our study would suggest a higher segregation in patients with UO than FO and no difference in the degree of integration. These differences can be explained by the scale at which the segregation and integration are quantified, with the clustering related to a local characteristic of densely connected close nodes, the modularity capturing the mesoscopic architecture with modules with up to few tens of nodes (Fig. 2).

Our results in deeply unconscious patients, as measured by the FOUR score during the first day of coma, suggest that these connectivity features, capturing segregation and integration, are not related to consciousness levels at least in this context of comatose patients after cardiac arrest. Their interpretation as neural correlates of consciousness should be considered with caution in those clinical contexts where the results of the functional network analysis is compared to the score obtained with clinical scales that can produce a high rate of mislabeling, (as for examples in patients with unresponsive wakefulness syndrome, Schnakers et al., 2009).

Temporal variability of the topological features was generally higher for patients with FO than UO and provided robust predictors for patients’ outcome. This parallels previous evidence that loss of consciousness results in a reduction of the repertoire of functional connectivity configurations during anesthesia (Lee et al., 2017, Solovey et al., 2015) sleep (Jobst et al., 2017) and in disorders of consciousness (Piarulli et al., 2016). Along the same lines, for the differences in terms of characteristic path length as a function of patients’ outcome, our results indicate that the variability of network configurations might constitute an antecedent to consciousness recovery at least in the specific experimental setting of coma after cardiac arrest.

EEG sources underlying functional connectivity

Our connectivity analysis is sensitive to neural synchronization of non-overlapping neural sources while disregarding their instantaneous dependencies. Despite the advantage of overcoming the volume conduction problem, we acknowledge its sensitivity to the presence of electrode bridges and artifacts of non-neural origins that would both produce spurious zero elements in the connectivity matrices. In addition, our analysis does not allow for the identification of the brain regions forming the functional networks observed at electrode level. In a combined EEG and fMRI study, the global EEG synchronization in the alpha frequency was located in brain areas involved in the ‘default mode network’ (Jann et al., 2009). Such networks have been suggested as a plausible candidate for neural correlates of consciousness and might underlie the differences we observed between patients with FO and UO during the first day of coma. This claim is consistent with previous functional neuroimaging studies showing that the level of functional disruption of the ‘default mode network’ (DMN) is related to the severity of the clinical condition during acute coma (Norton et al., 2012). Other studies have also suggested that the DMN is informative of the consciousness level in DOC patients (Vanhaudenhuyse et al., 2010) and that it is consistently altered across different varieties of unconscious states such as deep sleep and anesthesia (for a review see (Guldenmund et al., 2012).

Clinical differences in prediction analyses

Clinical differences based on path length time-variance were found in the discontinuity of the EEG (Table 2 – favorable outcome), in the dosage of Midazolam (Table 2 – unfavorable outcome) and in the presence of epileptiform activity. In particular, in the favorable outcome group, there were more patients with discontinuous EEG and low path length time-variance (i.e. below threshold) than patients with discontinuous EEG and high path length time-variance. Interpreting the relation between phase-based functional connectivity and the discontinuity (based on signal amplitude (Hirsch et al., 2013) is not straightforward. EEG discontinuity refers to a broad spectrum of background EEG categories with voltage activity of less than 10 μV along up to 50% of the recording. One explanation for this difference is that sustained periods of highly suppressed EEG would contribute to lower the temporal variability of the functional connectivity, reflected in a smaller path length time-variance. The higher value of Midazolam in unfavorable outcome patients with low vs high path length variance could be explained by the suppression of the EEG signal by the sedative agent, which reduces the repertoire of possible configurations.

Epileptiform activity occurred more often in patients with unfavorable outcome and low vs high time variance path length. The epileptic activity might result in a high and uniform phase coupling across electrodes leading to a stable pattern of connectivity over time and consistent with low time-variance path length.

Limitations and perspectives

In its present form, the implementation of this network analysis requires a semiautomatic rejection of noisy data following visual inspection from an experienced electrophysiologist. The evaluation of its predictive performance based on a fully automatic implementation goes beyond the scope of the present study, yet automation would represent a key step for its potential use in a clinical setting (Jas et al., 2017). One main limitation of the present study is that the – predominantly muscle – artifacts provoked the exclusion of 28% of the patients. This represents a higher rate in comparison to our previous quantitative EEG studies (Tzovara et al., 2016, Pfeiffer et al., 2018, Cossy et al., 2014), owed to the increased sensitivity of connectivity metrics to noise perturbations (as in similar studies of connectivity in DOC patients (Rizkallah et al., 2019). Considerably longer registration and administration of neuromuscular blocking agents to combat muscle noise could help circumventing such limitations in further studies. Future investigations will focus on the added value of combining several EEG markers, including power spectral analysis, optimal composition of the proposed topological features and clinical markers. One main challenge here is to improve the outcome prediction results without losing results interpretability, which may represent an obstacle for its use in a clinical routine.

Conclusions

We show for the first time that comatose patients following CA exhibit different functional network configurations on the first day of coma depending on their long-term clinical outcome. Both functional topological features and their variability over time provide early quantitative prediction of patients outcome. Patients with FO exhibit a wider set of configurations of functional networks than patients with UO during the first 24 h of coma. Future studies will investigate the evolution over time after CA of these functional network properties, how they parallel the clinical progression and how they are influenced by sedative agents and targeted temperature management.

Data analysis code

Scripts used for the connectivity analysis are available at: https://github.com/DNC-EEG-platform/ConnectivityComa

CRediT authorship contribution statement

Thomas Kustermann: Conceptualization, Methodology, Software, Formal analysis, Data curation, Writing - original draft, Visualization, Investigation. Nathalie Ata Nguepnjo Nguissi: Resources, Writing - review & editing, Data curation, Investigation. Christian Pfeiffer: Investigation, Resources, Writing - review & editing, Project administration. Matthias Haenggi: Resources, Writing - review & editing, Project administration, Data curation. Rebekka Kurmann: Resources, Data curation, Investigation. Frédéric Zubler: Investigation, Resources, Data curation, Writing - review & editing. Mauro Oddo: Resources, Writing - review & editing. Andrea O. Rossetti: Resources, Data curation, Writing - review & editing. Marzia De Lucia: Conceptualization, Methodology, Validation, Data curation, Writing - original draft, Supervision, Project administration, Funding acquisition.