Neural signature of coma revealed by posteromedial cortex connection density analysis.

Posteromedial cortex (PMC) is a highly segregated and dynamic core, which appears to play a critical role in internally/externally directed cognitive processes, including conscious awareness. Nevertheless, neuroimaging studies on acquired disorders of consciousness, have traditionally explored PMC as a homogenous and indivisible structure. We suggest that a fine-grained description of intrinsic PMC topology during coma, could expand our understanding about how this cortical hub contributes to consciousness generation and maintain, and could permit the identification of specific markers related to brain injury mechanism and useful for neurological prognostication. To explore this, we used a recently developed voxel-based unbiased approach, named functional connectivity density (CD). We compared 27 comatose patients (15 traumatic and 12 anoxic), to 14 age-matched healthy controls. The patients' outcome was assessed 3 months later using Coma Recovery Scale-Revised (CRS-R). A complex pattern of decreased and increased connections was observed, suggesting a network imbalance between internal/external processing systems, within PMC during coma. The number of PMC voxels with hypo-CD positive correlation showed a significant negative association with the CRS-R score, notwithstanding aetiology. Traumatic injury specifically appeared to be associated with a greater prevalence of hyper-connected (negative correlation) voxels, which was inversely associated with patient neurological outcome. A logistic regression model using the number of hypo-CD positive and hyper-CD negative correlations, accurately permitted patient's outcome prediction (AUC = 0.906, 95%IC = 0.795-1). These points might reflect adaptive plasticity mechanism and pave the way for innovative prognosis and therapeutics methods.

Introduction

Over the last few years, the posteromedial cortex (PMC) has received an increasing amount of attention. In fact, this architectonically discrete region, has been recognized as a critical site integrating an important range of multimodal information (Dehaene and Changeux, 2011). Actually, this highly dynamic functional core seems to participate in multiple transitional connectivity networks seemingly playing a critical role in internally/externally directed high-level cognition (Cavanna and Trimble, 2006, Leech and Sharp, 2014). In particular, converging data from physiological, pharmacological (Heine et al., 2012a, Heine et al., 2012b) and pathological models (Hannawi et al., 2015), suggest the implication of PMC and its long-range functional connections in conscious processing. For example, in pathological conditions, patients with disorders of consciousness (DOC) consistently demonstrated a reduced activity (He et al., 2014, Silva et al., 2010, Tsai et al., 2014) or diminished connectivity between this posterior brain structure and other cortical hubs (Hannawi et al., 2015, Vanhaudenhuyse et al., 2010, Qin et al., 2015, Wu et al., 2015), in particular the medial prefrontal cortex (mPFC) (Lant et al., 2016, Silva et al., 2015).

Interestingly, a growing body of literature on animal and human studies, suggests a significant heterogeneity in cytoarchitectonic (Vogt and Laureys, 2005, Vogt et al., 2006), structural (Parvizi et al., 2006, Zhang et al., 2014) and functional connectivity maps characterizing different sub-region of the PMC (Bzdok et al., 2015, Zhang and Li, 2013). Recent studies highlighted that the “metastable” functional connectivity detected in this region follows a complex ventral/dorsal-anterior/posterior gradient, partially overlapped across anatomically defined sub-regions (i.e. Precuneus (PreCu) and Posterior Cingulate Cortex(PCC)) (Bzdok et al., 2015, Cauda et al., 2010, Margulies et al., 2009, Zhang and Li, 2013). Nevertheless, it must be noted that the DOC neuroimaging literature traditionally explored PMC as a homogenous structure and failed to describe such a functional segregation in pathological conditions (Laureys et al., 1999, Norton et al., 2012, Silva et al., 2010, Vanhaudenhuyse et al., 2010). This important issue is probably due to seed-based approaches that are currently used in this setting to evaluate the functional connectivity among non-parcelled brain regions, by using correlation analyses of spontaneous fluctuations of brain activity in resting state conditions (Hannawi et al., 2015).

Therefore, we suggest that a better understanding of intrinsic PMC functional topology, (Silva et al., 2015) could significantly expand our understanding of how this cortical hub contributes to the generation and the maintenance of conscious awareness and might considerably improve DOC patient's clinical management. To explore this, we used a recently developed voxel-based unbiased approach that does not rely on a priori selection of the seed regions, named functional connectivity density (CD) (Tomasi and Volkow, 2010). This voxel-based method, accurately enables the identification of functional connectivity hubs and permit to specifically investigate within brain regions parcellation, in both healthy and pathological conditions. Thus, we aimed to investigated the functional impact of acute brain injures responsible of coma at the level of PMC and intended to study in this setting: (i) the specific interactions of the PMC anatomical (PCC and PreCu) or functional (ventral/dorsal gradient) sub-regions, with a distant cortical hub (mPFC) in resting state conditions (ii) a complete assessment of the whole range of increase/decrease of both positive/negative, i.e. corresponding to positive/negative correlation, connection patterns that could theoretically be detected by this approach, (iii) the impact of injury mechanisms (i.e. traumatic or anoxic), on brain functional connectivity patterns (iv) the prognostic value of functional connection density data for neurological recovery.

Materials and methods

Participants

Patients were included from three intensive critical care units affiliated with the University Teaching Hospital (Toulouse, France) between January 2013 and February 2014. We compared 27 patients, 15 with traumatic and 12 with anoxic brain injury, who met the clinical definition of coma (Glasgow Coma Scale score (Teasdale and Jennett, 1974) at the admission to hospital < 8, with motor responses < 6; age range: 19–70 years) to 14 approximately age-matched healthy controls (age range: 22–37 years). Patients underwent rs-fMRI scanning at least 2 days (4 ± 2days) after complete withdrawal of sedation and under normothermic condition. Standardized clinical examination was performed on the day of the scanning using the Glasgow Coma Scale and the Full Outline of Unresponsiveness (Wijdicks et al., 2005) and 3 months later using Coma Recovery Scale-Revised (Schnakers et al., 2008).

Image acquisition

In all participants, we acquired 11 min resting state fMRI using a 3 T magnetic resonance scanner (Intera Achieva; Philips, Best, the Netherlands). Two hundred and fifty multislice T2*- weighted images were retrieved with a gradient echo-planar sequence using axial slice orientation (37 slices; voxel size: 2 × 2 × 3.5 mm; TR = 2600 ms; TE = 30 ms; flip angle = 90°; FOV = 240 mm). A 3D T1-weighted sequence (in-plane resolution 1 × 1 × 1 mm, 170 contiguous slices) was also acquired in the same session, which was later used for visual assessment of the structural integrity of the PMC.

Pre-processing

The rs-fMRI data was preprocessed using SPM 8 (http://www.fil.ion.ucl.ac.uk/spm/) and CONN toolbox ver. 13f (http://www.nitrc.org/projects/conn) (Whitfield-Gabrieli and Nieto-Castanon, 2012). In order to reduce the motion effects on our data, we only included subjects characterized by motion parameters smaller than 3 mm translation and 3° rotation. First, the echo-planar images were realigned (motion corrected), slice-time corrected and normalized to the Montreal Neurological Institute echo-planar imaging template. Second, non-neuronal sources of noise were estimated and removed using the anatomical CompCor method (aCompCor) integrated in the CONN toolbox. Principal components of the signals from the white matter and the CSF voxels, alongside the motion parameters estimated during realignment were removed with regression. Finally, a temporal band-pass filter was applied to the residual blood oxygen level–dependent (BOLD) time course in order to obtain a low-frequency range of interest (0.008 Hz < f < 0.09 Hz).

Region of interest selection

Using a home-made MATLAB (MATLAB and Statistics Toolbox Release 2011a, The MathWorks Inc., Natick, Massachusetts, United States) script, the BOLD time series was extracted from voxels in two main regions of interest (ROI), the Posterior Medial Cortex and the Medial Prefrontal Cortex, defined by the Automated Anatomical Labeling atlas (Tzourio-Mazoyer et al., 2002) (voxel size 2 × 2 × 2 mm). The PMC (size = 12,862 voxels) consisted of the Precuneus L/R (size = 11,222 voxels) and the Posterior Cingulate Cortex L/R (size = 1640 voxels), and the mPFC (size = 13,389 voxels) comprised the Frontal Superior Medial L/R (size = 8373 voxels) and the Anterior Cingulate Cortex L/R (5016 voxels). As preliminary step, T2* mean images were used the extract the mean value of voxels in the PMC region, and a two-sample t-test was performed to compare values between the control and the patient group.

Additionally, to investigate if the potential changes in connection density are specific to the PMC-mPFC interactions, we have also included the bilateral Calcarine L/R (size = 7134 voxels; part of the primary visual cortex) as a control region, currently not considered to be relevant for conscious awareness.

Voxel-based connection density

The data analysis pipeline is presented in Fig. 1.

Fig. 1

Overview of the data analysis pipeline. 1) First, the rs-fMRI data was pre-processed using SPM8 and CONN13f, respectively. 2) Using a home-made MATLAB script, the BOLD time series was extracted from voxels in two main regions of interest then used in 3) calculation of Pearson correlation coefficients between the BOLD time course of all the possible pairs of voxels from PMC and mPFC. Then, a subject-specific threshold of p ≤ 0.05 was applied to include only the significant connections in further analysis. 4) The obtained significant connections were split on positive and negative (based on the sign of the normalized z coefficients), binarized and used to 5) obtain the density of connection between PMC and mPFC voxels. 6) A Z-score was calculated for each single PMC voxel as explained in the figure. 7) Voxels with a sig. Z score were summed and characterized as hypo/hyper-CDP and hypo/hyper-CDN as presented in the figure.

Pearson correlation coefficients were computed between the BOLD time course of all the possible pairs of voxels from PMC and mPFC. Correlation coefficients were then normalized using Fisher's r-to-z transformation. A subject-specific threshold of p ≤ 0.05 was applied to each correlation coefficient in order in order to retain only a subset of connections with higher in further analysis. We labelled the obtained connections as positive and negative depending on the sign of the obtained z coefficients and we treated them separately. Importantly, to avoid the artificial induction of negative correlations (anticorrelations) by the global signal regression (Murphy et al., 2009, Weissenbacher et al., 2009) we adopted the CompCor method in the preprocessing step, which has been reported as a reliable approach for the exploration of both positive and negative correlations (Chai et al., 2012, Chai et al., 2014).

The resulting z coefficients were binarized and then summed to obtain the density of connections between every voxel belonging to PMC and every voxel within the mPFC. These would represent the number of significant positive and negative connections from PMC to mPFC. The number of significant connections obtained for a single PMC voxel was then used to compute a Z-score for every voxel in the patient group using the mean (M) and the standard deviation (SD) of the healthy control group. This will represent the significant deviation of each voxel from the corresponding control.

Thus, the number of PMC voxels with a Z-score higher or equal to 2SD or less or equal to − 2SD (for positive and negative connections) were counted for each patient individually and further used to investigate the changes (significant deviations) in functional connectivity between the PMC to the mPFC. The same analysis was repeated for the control pathway (i.e. representing the connections between the bilateral Calcarine and the mPFC).

Furthermore, in order to facilitate the discussion of our results (Fig. 1), a significant decrease (Patients < Controls, Z-score ≤ − 2) in connection densities is designated as “hypo” connection density or hypo-CDP and hypo-CDN, for positive and negative correlation based connection densities, respectively. The opposite - increase in connection densities (Patients > Controls, Z-score ≥ 2) is denoted by “hyper” connection density e.g. hyper-CDP and hyper-CDN for positive and negative correlation based connections, respectively.

Lastly, to allow an easier description of our findings and to test our hypothesis on the differences in functional connectivity changes between the PreCu and PCC, we present the results from these two ROIs separately.

Spatial homogeneity

In addition to the individual results, we were interested in the spatial similarity (overlap) of these changes between patients within groups (intra-group) and between different coma aetiology groups (inter-group). This allowed us to investigate if these changes were spatially scattered thus highly heterogeneous between subjects or if they were organized in functional clusters shared among multiple participants within/between the two aetiology groups.

The intra-group spatial congruity of single voxel Z-score results was explored at two thresholds: a criteria of 33% and 67% was used to define the total number of disconnected voxels spatially shared between 1/3 and 2/3 of the patient group, respectively. This was explored in the group with all the patients and with two coma aetiology groups separately.

To explore the inter-group similarity, we have calculated the Jaccard similarity coefficient (index) to further test the spatial overlap in connectivity changes between these two groups. The Jaccard index (JI) is defined as the intersection divided by the union of the number of voxels representing significant changes in PMC-mPFC connection density in given groups. The values of JI range from 0 (0%), indicating no overlap, to 1 (100%) suggesting full spatial overlap.

Statistical analysis

We have conducted the Mann-Whitney U test to investigate the differences in PMC-to-mPFC connection density changes between the PreCu and PCC. The same test was used to compare the changes in functional connectivity between our two coma aetiology groups – traumatic and anoxic brain injury patients. This analysis was also repeated for the control pathways. The resulting p values were corrected for multiple comparison using a false discovery rate (FDR p = 0.05).

We have also calculated the effect size for the Mann-Whitney U test by dividing the z value (test statistic) by the square root of N (sample size). This is analogous to Cohen's d for parametric group testing.

Spearman's correlation analysis was performed to explore the link between the number of PMC voxels with significant changes in connection density and the CRS-R score. The same correlation analysis was done to test the association between the connection density changes in the control pathway - Calcarine-mPFC - and the CRS-R score. The resulting p values were corrected for multiple comparison using a false discovery rate (FDR p = 0.05).

Additionally, in order to explore the non-linear association between the outcome and the significant changes in PMC-mPFC connection density, and thus the predictive value of these changes in relation to recovery, we have conducted the binary/logistic regression analysis. Patients were divided into two outcome groups (based on the CRS-R score): good outcome - comprising patients who had recovered or progressed to minimally conscious state (MCS) (N = 12); and the bad outcome – incorporating patients with unresponsive wakefulness syndrome/vegetative state (UWS/VS) (N = 15). This analysis was not done separately for the two aetiologies due to small sample size.

All of the above mentioned analyses were performed using the IBM SPSS (IBM SPSS Statistics for Windows, Version 20.0. Armonk, NY: IBM Corp.) statistical package.

Results

Induced functional topological changes in the PMC

Anatomical sub-regions — precuneus and posterior cingulate cortex

In Fig. 2 we report the induced changes in functional connectivity of the entire PMC and the differences between its anatomical sub-regions, the PreCu and PCC, in the patient group compared to controls. The changes in connection density are presented as the percentage of voxels more or less connected within a given region/sub-region in order to take into account the differences in size (total number of voxels) of different ROIs.

Fig. 2

PMC (PreCu/PCC) connection density Z-score results for individual patients. Patients exhibit a significant number of hypo-CDP (panel A), hyper-CDP (panel B) and hyper-CDN (panel D) in comparison with the control group. There seem to be some differences depending on the coma aetiology. The x axis represent the number of patients (each patient is coded with the same number in all four panels: A, B, C, D). The y axis reflects the percentage of voxels with a sig. Z-score (− 2 ≤ Z ≥ + 2) of the total number of voxels in the PMC (12,862 voxels) and the PreCu (11,222 voxels) and PCC (1640 voxels) separately. The changes in PreCu and PCC are presented separately for two coma aetiologies using different colours - dark blue/light blue for TBI and dark green/light green for anoxic BI.

The Mann-Whitney U test with percentage-wise results showed that the PCC exhibited significantly more voxels with both hypo-CDP (U = 202, p = 0.004; Fig. 3.A.1.) in comparison to PreCu, with a medium effect size (r = 0.38), which is in accordance to previous research. The latter showed a tendency toward a higher number of voxels with hyper-CDP (Fig. 3.A.2) in comparison to PCC, but these results were not statistically significant (U = 261, p = 0.074; small effect size (r = 0.24)). PCC did not significantly differ from PreCu in hyper-CDN (U = 248, p = 0.044; small effect size (r = 0.27) after the multiple comparison correction (p > 0.004; FDR corr.).

Fig. 3

Differences in changes in connection density between PreCu and PCC and between traumatic and anoxic brain injury. Panel A) PCC showed significantly more hypo-CDP (sub-panel 1) in comparison with the PreCu. Panel B) Traumatic BI patients had more hyper-CDP (sub-panel 2) within the PCC and hyper-CDN voxels in both PreCu and PCC (sub-panel 4) in comparison to anoxic BI patients; Boxplots represent medians with interquartile range and whiskers signify minimum and maximum values (excluding the outliers) (*p < 0.05, (**p < 0.005, ns: nonsignificant).

No significant differences were observed in the number of voxels with hypo-CDN between the PCC and the PreCu (U = 256.5, p = 0.055; small effect size r = 0.26)), which is not unexpected given that these changes were barely present in both sub-regions (Fig. 3.A.3).

Functional sub-regions — ventro/dorsal segregation of the PMC

In order to investigate the spatial architecture of the observed changes in the density of connections, in Fig. 4 we report the spatial topography of the obtained connections. The figure represents the summed changes in functional connectivity in the entire patient group. The individual spatial topography maps of these changes (with aetiology and prognostic information) are presented in the Supp. Fig. 1 and Supp. Fig. 2.

Fig. 4

Spatial maps of changes in PMC-to-mPFC connection density in the patient group. Patients seem to show significant hypo-CDP changes in the ventral portions of both PreCu and PCC (panel A). Hyper-CDP voxels are more specific to the dorsal sub-regions of the PreCu and PCC, and slightly border on the ventral PCC (panel B). Hypo-CDN changes are barely noticeable (panel C), and hyper-CDN voxels seem wide-spread and highly heterogeneous between patients (panel D). Individual results of patients are summed and presented in a single image; Purple and blue outlines reflect the borders of the anatomical PreCu and PCC respectively. Gradient bars (%) reflect the percentage of patients sharing the same voxel with a sig. Z-score (hypo/hyper CDP/CDN) at given anatomical location (spatial homogeneity).

It can be noted that: (i) the ventral PCC and the ventral PreCu seemed to form a functional cluster of hypo-CDP voxels (Fig. 4.A), (ii) the hyper-CDP voxels were primarily located in the dorsal PreCu and a part of the dorsal PCC (Fig. 4.B.), alongside a portion of the ventral PCC (Fig. 4.B.), (ii) the hyper-CDN voxels seemed to be widespread, suggesting high inter-individual variability, and covered large portions of both PMC sub-regions (Fig. 4.D), including almost the entire PCC (97% for the entire patient group, see Table 1), (iv) finally, the hypo-CDN PMC voxels were barely present and mostly found in the posterior ventral PreCu (Fig. 4.C).

Table 1

Spatial homogeneity of changes in PMC connection density in the patient group. Threshold – at least 33% or 67% of subjects spatially share the same voxel with sig. Z-score. The percentages present the proportion of voxels with significant Z-score out of the total number of voxels in PMC spatially shared between patients.

Group		Hypo-CDP	Hyper-CDP	Hyper-CDN
PatientsN = 27	No threshold	8% (1054)	51% (6599)	68% (8803)
	≥ 33%	3% (347)	3% (396)	9% (1138)
	≥ 67%	0% (7)	0% (9)	0% (5)
TBIN = 15	No threshold	8% (1011)	46% (5866)	67% (8668)
	≥ 33%	2% (217)	7% (942)	23% (2959)
	≥ 67%	0% (1)	0.4% (51)	0.8% (102)
ANOXICN = 12	No threshold	6% (730)	29% (3740)	33% (4305)
	≥ 33%	4% (473)	2% (244)	3% (339)
	≥ 67%	1% (143)	0% (3)	0% (6)

The spatial topography of our results indicated a possible significant spatial overlap between voxels with hypo-CDP and hyper-CDN which lead us to calculate the exact number of these intersecting voxels. Thus, our calculation demonstrated a complete overlap of these two types of disconnections, which was equivalent to almost the total number of hypo-CDP voxels (1048 voxels or 11% of summed hypo-CDP and hyper-CDN voxels; see Table 1.) in the entire patient group. However, these changes were not homogeneous and were characterized by important inter-individual differences (SD = 155.3), as demonstrated in the Supp. Table 1.

Control pathways — Calcarine-mPFC

The changes in connection density are fully presented in the Supp. Fig. 3. The results showed only a small percentage of Calcarine hypo-CDP voxels representing reduced connectivity between this region and the mPFC (median = 0.04%, min = 0, max = 0.7%) in patients. The same result was found for the hypo-CDN voxels (median = 0, min = 0, max = 0.2%). Nevertheless, hyper-CDN and hyper-CDP voxels were present in this control pathways (median = 22%, min = 0.08%, max = 63% and median = 4%, min = 0.04%, max = 22%, respectively).

Impact of brain injury mechanisms

PMC-mPFC connection density changes

The impact of the patients' aetiology on the observed changes reported in Fig. 3.B. In accordance with the single-subject level (Fig. 2), the group-level statistical analysis showed significant differences between the two aetiologies in the percentage of PMC voxels with hyper-CDN (U = 32, p = 0.004) with large effect size (r = 0.55). There were no significant difference in hypo-CDP (U = 63, p = 0.200; small effect size (r = 0.25)), hyper-CDP (U = 49, p = 0.045, FDR correction p < 0.005; medium effect size (r = 0.39)) and hypo-CDN voxels (U = 66, p = 0.249, small effect size (r = 0.23)). These differences were reflected in a significantly higher percentage of hyper-CDN voxels in both PreCu (U = 32, p = 0.004) and PCC (U = 33.5, p = 0.005) in TBI in comparison to the anoxic patient group (Fig. 3.B.4). These differences were associated with large effect size in both cases (PreCu, r = 0.55; PCC, r = 0.53). Additionally, the results showed a higher percentage of hyper-CDP voxels in the PCC sub-region in the TBI group (U = 27, p = 0.001; Fig. 3.B.2.) with large effect size (r = 0.59), which was not reflected in the comparison done using the total percentage of PMC hyper-PCD voxels (p > 0.005, FDR corr.).

There were no significant differences between these two groups in the hyper-PCD in the PreCu sub-region (U = 50, p = 0.053; medium effect size (r = 0.38)), hypo-CDP in PreCu (U = 64, p = 0.213; small effect size (r = 0.24)) and PCC (U = 60, p = 0.148; small effect size (r = 0.28)), hypo-CDN in PreCu (U = 73, p = 0.417), small effect size (r = 0.16) and PCC (U = 62, p = 0.118; medium effect size (r = 0.30)).

Connection density changes in control pathways

Additionally, we wanted to test if the hyper-CDP and hyper-CDP voxels, found in the control pathways, were more specific to one of the aetiologies. The group-level statistical analysis showed a significantly higher percentage of Calcarine-mPFC hyper-CDP ((U = 37, p = 0.009; large effect size (r = 0.5)), and hyper-CDN voxels ((U = 33, p = 0.004; large effect size (r = 0.54)) in traumatic in comparison to anoxic brain injury patients (Supp. Fig. 3, Supp. Fig. 4).

The intra-group spatial similarity in voxels with changes in connection density are presented in Table 1. The hypo-CDN voxels were not included in the analysis, because there were no significant differences between the control and patient group. The results are reported without threshold and with the threshold of 33% and 67% of patients having the same voxel changes.

First, there seemed to be a lower spatial congruity in the entire group of patients, as opposed to the groups of same aetiology (Table 1). It can be observed, that even though there is spatial similarity between 1/3 of patients in the disconnected voxels, this number seems very low when considering 2/3 of the patient group.

Second, in comparison with traumatic patients, anoxic brain injury patients seem to have a slightly higher spatial homogeneity in PMC voxels with hypo-CDP. Third, TBI patients showed a higher intra-group similarity in hyper-CDP and hyper-CDN in comparison to anoxic patients. The latter changes were present even at the threshold of 67%, meaning that > 2/3 of patients were spatially sharing voxels with hyper negative connection density. The spatial homogeneity for separate aetiologies is presented in Fig. 5 with a threshold of 33% of patients needing to have the same hypo-CDP and hyper-CDN changes. An overlap between these two types of disconnection was additionally calculated for each subgroup showing an overlap of 1005 voxels in the TBI group (SD = 194; 10% of hypo-CDP/hyper-CDN voxels), and 665 voxels (SD = 67.6; 13% of hypo-CDP/hyper-CDN voxels) in the anoxic group.

Fig. 5

Intra-group spatial homogeneity differences between the traumatic and anoxic brain injury patients (threshold 33%). Anoxic BI patients show a higher intra-group similarity in hypo-CDP in comparison to TBI patients, while TBI patients show a higher spatial congruity in hyper-CDN voxels. Purple and blue outlines reflect the borders of the anatomical PreCu and PCC respectively. Gradient bars (%) reflect the percentage of patients sharing the same voxel with a sig. Z-score (hypo/hyper CDP/CDN) at given anatomical location. The minimum spatial homogeneity is set to at least 33% of patients in a given group.

The Jaccard index showed a similarity of 65% for hypo-CDP, 46% for hyper-CDP and 47% for hyper-CDN spatial maps between traumatic and anoxic BI patients. The similarity between the spatial maps of anoxic BI and the spatial maps of all patients (including TBI) was 69% for hypo-CDP, 57% for hyper-CDP and 49% for hyper-CDN voxel results. Finally, the comparison between the TBI and the group of all patients, suggested a similarity of 96% for hypo-CDP, 89% for hyper-CDP and 98% for hyper-CDN spatial maps. These results, along the spatial homogeneity results within different sub-groups, indicates that the results are more spatially heterogeneous and more widely distributed in the TBI group in comparison to the anoxic groups, further confirming differences in the pathological processes.

Prognostic value

Spearman correlation analysis with PMC-mPFC voxels

The number of PMC voxels with hypo-CDP showed a significant negative association with the CRS-R score assessed 3 months after the coma onset (rs = − 0.72, p = 0.00002), indicating that if a patient had a higher number PMC voxels with hypo-CDP he was less likely to recover (Fig. 6.A). We did not find a significant correlation between hyper-CDP (rs = 0.29, p = 0.14) and hyper-CDN (rs = − 0.35, p = 0.07) voxels and the CRS-R score for the entire patient group. The number of voxels with hypo-CDN was not included in the analysis, because there were no significant differences between the control and patient group.

Fig. 6

The prognostic value of changes in PMC-to-mPFC connection density. Panel A) There was a sig. Negative link between the number of PMC voxels with hypo-CDP and the patient outcome (rs = − 0.72; p = 0.00002); Panel B) TBI patients showed a significant negative association between the number of voxels with hypo-CDP (rs = − 0.80; p = 0.0004; panel A), hyper-CDN and the CRS-R score (rs = − 0.86; p = 0.00005); Panel C) In TBI patients, there was highly significant negative association between the number of spatially overlapped hypo-CDP and hyper-CDN voxels and the 3-month outcome (rs = − 0.73, p = 0.002; Fig. 6.C). The x axis represents the 3-month CRS-R score, the y axis represent the number of voxels with changes in connection density.

At the aetiology level, it is worth noting that in the traumatic brain injury group, we found a highly significant link between number of hypo-CDP voxels (rs = − 0.80, p = 0.0004) and the CRS-R score (Fig. 6.A). Interestingly, a highly statistically relevant association was found between the number of hyper-CDN voxels and the CRS-R score (rs = − 0.86, p = 0.00005). Therefore, TBI patients with a high number of hypo-CDP and hyper-CDN PMC voxels had less chance for recovery 3 months after the initial fMRI scan (Fig. 6.B). The analysis in the anoxic group did not show a significant association between the CRS-R score and the hyper-CDP (rs = 0.43, p = 0.16), hyper-CDN (rs = 0.20, p = 0.54) nor the hypo-CDP voxels (rs = − 0.38, p = 0.22). However, this result should be interpreted in caution due to small sample size of anoxic BI patients.

Additionally, given the significant spatial overlap between voxels with hypo-CDP and hyper-CDP voxels, and the significance of these voxels in relation to the neurological recovery in TBI patients, we performed the same analysis on overlapping voxels. The results showed the same highly significant negative association between the number of spatially overlapping hypo-CDP and hyper-CDN voxels and the 3-month outcome (rs = − 0.73, p = 0.002; Fig. 6.C).

Spearman correlation analysis with the control pathway voxels

The Spearman correlation analysis did not show any significant association between the CRS-R score assessed 3 months after the coma onset and any of the number of voxels with significant connection density changes in the Calcarine-mPFC control pathways (FDR p > 0.05). This was true for the entire group of patients and for the anoxic and traumatic BI groups separately.

Logistic regression analysis with PMC-mPFC voxels

A logistic regression analysis was conducted to ascertain the effects of the number of hypo-CDP, hyper-CDN and hyper-CDP PMC voxels on the likelihood that the patients had a good outcome (recovery and MCS) or bad outcome (UWS/VS). Due to significant collinearity between the hyper-CDP and hyper-NCD (r = − 0.91), two models were tested, the first including hypo-CDP and hyper-CDP and the second hypo-CDP and hyper-CDN as independent/predictor variables.

The first model was significant (chi square = 10.204, p = 0.006 with df = 2), however the number of hyper-CDP voxels did not significantly contribute to the prediction (Wald, p = 0.296). Thus, we decided to keep the second model, where both predictors had a significant effect on the outcome.

The full model against the constant only model was statistically significant, suggesting that the predictors as a set reliably distinguished between patients with good and bad outcome (chi square = 16.701, p < 0.001 with df = 2). The model explained 62% (Nagelkerke's R2) of the variance in outcome, and correctly classified 78% of cases (sensitivity – 80%; specificity – 75%). The Wald criterion demonstrated that both variables, hypo-CDP (p = 0.02) and hyper-CDN (p = 0.044) had a weak but significant contribution to prediction. Increase in both predictors was associated with decreased likelihood of good outcome. Accordingly, Exp (B) indicates that the for every extra unit (number of voxels) of hypo-CDP and hyper-NCD, the odds of recovery decrease by a factor of 0.988 and 0.999, respectively, all other factor being equal.

Discussion

Our findings demonstrate that a comparative analysis of voxel-wise connection density disruptions between the PMC and the mPFC, might constitute a reliable and fine-grained study approach, permitting to accurately address the brain functional network changes, which are critically related to consciousness abolition during coma. Crucially, built upon statistical analysis at a single subject level - a necessary condition to make reliable inferences in individual patients - this voxel-wise connection density approach, enabled the description of a novel PMC/mPFC dysconnectional taxonomy, permitting the identification of specific markers related to brain injury mechanism and useful for neuroprognostication. The current study, confirms and expands previous finding through seed-based analysis method (Silva et al., 2015), and permits: (i) a precise topological description of the PMC sub-regions (ventral PreCu and the PCC) who seems implicated in the consciousness impairment observed in this setting, (ii) identification of useful and accurate neuroimaging biomarker for prognostication in this setting, as the number of PMC voxels with hypo-CDP to mPFC, had significant predictive value in relation to patients neurological outcome assessed 3 months after the coma onset in the current study. It is worth noting that such adaptive changes, characterized by a reduction of CDP, were rarely observed in the control pathway (Calcarine/mPFC). Crucially, the number of voxels with any significant connection density changes within this control pathway, was not correlated with patient's neurological outcome. These points strongly support the critical role of PMC/mPFC functional connections in consciousness emergence, and consolidate the specific use of PMC/mPFC functional connectivity to predict coma patient's neurological outcome.

Our data strongly support the concept of a significant structural and functional PMC heterogeneity, which could be implicated in the “tuning” of the whole-brain metastable status (Hellyer et al., 2015). Furthermore, our results are in line with the hypothesis of a dysfunctional mesocircuit (Schiff, 2010) in DOC patients, build upon anterior and posterior critical disruption of neural communication. For example, tract-tracing studies conducted in non-human primates (Vogt et al., 2005) and diffusion tensor tractography in humans (Hagmann et al., 2008) have clearly identified structural connections between dorsal PCC to the mPFC along the cingulum bundle (Greicius et al., 2009). Interestingly, structural injuries within this tract have been described as related to neurological outcome of severely brain injured patients (Fernández-Espejo et al., 2012).

Furthermore, our findings suggest a higher prevalence of hypo-CDP voxels in comatose patients, specifically located in the ventral subdivision of both the PCC and the PreCu. This result seems compatible with the findings of the rs-fMRI study of Silva et al. (2015), suggesting stronger functional connectivity between PCC-centred seed and the mPFC in patient who eventually recovered. However, our voxel-wise connection density method further elucidated the importance of the total number of PMC voxels in prediction of patient's outcome and allowed us to explore the intra-regional topological changes in connectivity which is not possible in seed-based analysis. Additionally, previous resting-state fMRI (Bzdok et al., 2015, Cauda et al., 2010, Margulies et al., 2009, Vogt et al., 2006, Zhang and Li, 2013), and DTI studies (Zhang et al., 2014) indicate that both these PMC subregions are densely connected to other brain regions that participate to DMN network and appears to be implicated in internally directed cognition (i.e. self-referential information processing, social cognition, mind-wandering, episodic memory retrieval) (Andrews-Hanna et al., 2010, Cavanna and Trimble, 2006, Fox et al., 2015, Leech and Sharp, 2014). In contrast with this hypo-CDP restricted to ventral parts of PreCu/PCC, we demonstrated that dorsal PreCu and anterior dorsal (and a part of ventral) PCC, encompasses specific voxels with hyper-CDP in coma patients. These dorsal PMC structures are known to support and connect to other regions involved in externally directed cognitive processes (Margulies et al., 2009, Zhang and Li, 2013) such as spatially-guided behaviour, visual/motor mental imagery (Cavanna and Trimble, 2006), control of attentional focus (Leech and Sharp, 2014) and high-level sensorimotor functions (Balestrini et al., 2015, Herbet et al., 2015). Overall, the observed rostro-caudal functional segregation between increase and decrease of CDP, seems in line with the hypothesis of the central role played by an imbalance between internal and external awareness systems in the genesis of consciousness disorders (Demertzi et al., 2013, Di Perri et al., 2014, He et al., 2014). Actually, as described in chronic DOC patients (He et al., 2014), this imbalance could reflect overcompensation in the external network due to the loss of input from internal self-networks.

Additionally, we sought for functional signatures of the brain injury mechanisms that were responsible of coma. Regarding, the hypo-CDP voxels, our findings indicate that within the PMC, hypo-CDP voxels were slightly more pronounced in anoxic patients (although not statistically significant), probably representing a neural underpinning of the overall worse prognosis of this group (Horsting et al., 2015, Koenig et al., 2014). Furthermore, in line with recent studies, suggesting that hyperconnectivity patterns are a common network response to traumatic brain injury (Bharath et al., 2015, Hillary et al., 2014, Hillary et al., 2015, Stevens et al., 2012) we observed that hyper-CDP (in the PCC sub-region) and hyper-CDN voxels were more frequently detected in the traumatic group. This was true for the PMC-mPFC connectivity patterns, but also for the control pathway, representing connections between the Calcarine and the mPFC.

We hypothesize that the higher prevalence of hyperconnectivity patterns in this group could be the consequence of compensatory brain plastic processes (Di Perri et al., 2014) possibly reflecting resilient connections massively disrupted, but not totally interrupted by traumatic brain injury (i.e. diffuse axonal injury). These resilient connections might have engaged residual critical neural resources, otherwise normally distributed through efficient brain network connections, here disrupted by brain injury, resulting in impaired conscious processing (Hillary et al., 2015). Furthermore, these brain network signature of acute brain injury mechanisms was confirmed by spatial homogeneity analysis (thresholding and Jaccard index), which indicated higher spatial similarity of identified changes in voxel connection density between subjects with same aetiologies in comparison to the whole patient group. However, despite higher intra-group similarity, these changes seemed more spatially scatted in the TBI patients group in comparison to anoxic patients, which is in accordance with the heterogeneity of brain lesion found in this group of patients.

Finally, it is worth noting, that we specifically observed in the traumatic brain injury group, that the loss long-range connection density of positive correlations (hypo-CDP) between the PMC and the mPFC, was balanced with a spatially overlapping increase of negative coactivation (hyper-CDN). The total amount of the later was associated to patient recovery and could depict maladaptive brain plasticity processes. We hypothesize that this potentially compensatory increase of negative connection density may indicate some aberrant inhibitory processes or competitive mechanisms of information processing (Fox et al., 2009, Fox et al., 2012, Gee et al., 2011, Gopinath et al., 2015, Liu et al., 2015), resulting in the breakdown of more efficient network organization (i.e. PMC hypo-CDP). This is further supported by the unfavourable outcome of our patients with a higher number of overlapping PMC voxels with decreased positive connections but yet increased negative connections with mPFC.

Interestingly, the same hyper-CDP and hyper-CDN patterns in the control pathways were not associated with the outcome in TBI patients. This result indicates that despite of the more global pattern of these changes, the increase of negative connection density seems maladaptive only when present in regions critical to conscious processing.

There are several limitations in this study. First the number of patients is relatively small and our results need to be validated in larger cohorts. Second, we performed a study of density of connections between each PMC voxel and all of the MPFC voxels, aiming to explore in detail, the topology and the functional characteristics of the connections established from PMC which was assumed homogeneous. Therefore, future studies should also address the functional segregation and heterogeneity within mPFC in this setting. In future studies, brain damage in the form of lesions in traumatic brain injury should be assessed in the pre-processing step, as though not properly treated lesions can potentially give rise to artificially induced connectivity.

In summary, a complex pattern of decreased and increased connections was observed, and the topography of these changes, seemed to be in agreement with the hypothesis of network imbalance between internal/external processing systems, within PMC during acquired disorders of consciousness. We report a significant link between the PMC (ventral parts of both the PCC and PreCu) and the mPFC functional connectivity and patient recovery. Actually, the number of PMC voxels with hypo-CDP showed a significant negative association with the CRS-R score and a “negative” predictive value in relation to good outcome assessed 3 months after the coma onset, notwithstanding aetiology. Additionally, traumatic brain injury specifically appeared to be associated with a greater prevalence of hyperconnected positively and negatively correlated voxels, and the total amount of latter, was inversely associated with patient neurological outcome. This point might reflect a maladaptive plasticity mechanism through a resilient functional network and an inefficient redistribution of remaining resources and pave the way for innovative prognosis and therapeutics methods in this challenging setting.

The following are the supplementary data related to this article.

Supp. Fig. 1

Individual traumatic brain injury patients' connection density spatial maps. The spatial maps are presented separately for each patient (numbered as in figures embedded in the main text – Fig. 2, Fig. 6.). The changes in connection density are presented separately for hypo-CDP, hyper-CDP and hyper-CDN and the two anatomical PMC sub-regions – PreCu and PCC. The patients' outcome is marked beneath the patients' code – REC – Recovery, MCS – Minimally conscious state; VS/UWS – Vegetative state/Unresponsive wakefulness syndrome.

Supp. Fig. 2

Individual anoxic brain injury patients' connection density spatial maps. The spatial maps are presented separately for each patient (numbered as in figures embedded in the main text – Fig. 2, Fig. 6.). The changes in connection density are presented separately for hypo-CDP, hyper-CDP and hyper-CDN and the two anatomical PMC sub-regions – PreCu and PCC. The patients' outcome is marked beneath the patients' code – REC – Recovery, MCS – Minimally conscious state; VS/UWS – Vegetative state/Unresponsive wakefulness syndrome.

Supp. Fig. 3

Differences in Calcarine-mPFC connection density changes between traumatic and anoxic brain injury. TBI patients showed significantly more Calcarine-mPFC hyper-CDP (sub-panel 2) and hyper-CDN voxels (sub-panel 4) in comparison with the anoxic BI patients. Boxplots represent medians with interquartile range and whiskers signify minimum and maximum values (excluding the outliers) (*p < 0.05, **p < 0.005, ns: nonsignificant).

Supp. Fig. 4

Connection density changes in the Calcarine-mPFC pathway (threshold 33%). TBI patients seem to show more widespread changes and a higher spatial congruity in hyper-CDN and hyper-CDP voxels in comparison to anoxic BI. Gradient bars (%) reflects the percentage of patients sharing the same voxel with a sig. Z-score (hyper CDP/CDN) at given anatomical location within the Calcarine. The minimum spatial homogeneity is set to at least 33% of patients in a given group.

Supp. Fig. 5

Logistic regression model. The graph presents true positive (sensitivity) on the vertical axis and true negative (1-specificity) on the horizontal axis. The area under the ROC curve (AUC) ranges from 0.5 to 1.00 with larger values indicative of a better fit. The area under the curve is 0.906 with 95% confidence interval (0.795–1.00). Also, the AUC is significantly different from 0.5 (50/50 chance) with p-value = 0.000 suggesting that the logistic regression classifies the group significantly better than by chance. PPV – positive predictive value; NPV – negative predictive value; AUC – area under the curve; CI – confidence interval; Std. error – standard error.

Supp. Table 1

Individual patient results – the total number of hypo/hyper-CDP and hypo/hyper-CDN voxels for PMC (Panel A), PreCu (Panel B) and PCC (Panel B). The overlapping hypo-CDP and hyper-CDN voxels are presented only for PMC (Panel A). Number of voxels per region: PMC - 12,862 voxels, PreCu - 11,222 voxels, PCC - 1640 voxels. Abbreviations: TBI – traumatic brain injury, ANOXIC – anoxic brain injury; N – number of voxels with significant Z-score. REC – Recovery, MCS – Minimally conscious state; VS/UWS – Vegetative state/Unresponsive wakefulness syndrome. The patients are numbered as in figures embedded in the main text.

3D Neuroimaging data

The nifti images present summed binarized images with significant changes in PMC-to-mPFC connection density in the patient group (allPAT), and the anoxic BI (ANOX) and traumatic BI (TBI) groups separately. The values range from 0, meaning 0 patients have a significant change in a given voxel, to 27 (patient group), 12 (anoxic BI group) and 15 (TBI group), representing the maximum number of patients in each of the groups and thus reflecting complete spatial homogeneity in connection density changes in a given voxel. Data with “hypo-CDP” and “hyper-CDP” presents decrease and increase in PMC-to-mPFC connection densities based on positive correlation. Data with “hypo-CDN” and “hyper-CDN” reflects decrease and increase, respectively, in connection densities based on negative correlation.

Funding

This study was supported by “Association des Traumatisés du Crâne et de la Face” (2014-A3); “Institut des Sciences du Cerveau de Toulouse” (2016-21); “Institut National de la Santé et de la Recherche Médicale” (2016-DEVIN). The funding sources had no role in the study design, data collection, data analysis, data interpretation, or writing of this report.