Thirty-five hospitalized DOC patients (13 females and 22 males, 51.7 ± 14.1 years old) were included in this study (Table 1). These patients did not show cortical lesions in the precuneus. Data collection was completed in the Zhengzhou Central Hospital Affiliated to Zhengzhou University. With the CRS-Rs, the score of the patient’s consciousness was assessed (Giacino et al., 2004) and all CRS-R were performed by trained physicians at T0, T1, T2, and T3. Each patient was evaluated by the same physician. In this study, we divided the MCS patients into MCS+ and MCS−.

The patients were divided into UWS (three females and 12 males, 51.0 ± 9.6 years old) and MCS (10 females and 10 males, 52.3 ± 16.9 years old) groups according to their CRS-R scores before treatment. DOC patients who had previously undergone tDCS treatment within the past month were excluded from the present study. The patients with aneurysm clips, pacemakers, other implanted devices, or other drugs that may influence the patients’ EEGs were also excluded from this study. This study was approved by the ethics committee of the Zhengzhou Central Hospital Affiliated to Zhengzhou University. Assessments on whether any of the included patients had been returned to the normal level were completed in June 2019.

A 32-channel EEG recorder was used for data collection (Nicolet EEG V32, Natus, United States; Channel names: Fp1, Fp2, F3, F4, C3, C4, P3, P4, O1, O2, F7, F8, T7, T8, P7, P8, A1, A2, Fz, Cz, Pz, FC5, FC1, FC2, FC6, POz, CP1, CP2, CP6, Fpz, CP5, and Oz). Resting-state EEGs were recorded at four time-points (T0, T1, T2, and T3). The 32 channels EEG signals were continuously recorded covering the whole scalp consistent with the International 10/20 system with the Cz as reference. The signals were notch filtered at 50 Hz in the recorder. The sampling rate of EEG signals was set to 1,000 Hz. During the data collection process, the impedance of all electrodes was kept no more than 5 KΩ and the patients were behaviorally awake. The EEG recording lasted approximately 10 min for each session.

All EEG data from the 32 channels were retained for analysis using EEGLAB (Delorme and Makeig, 2004). EEG signals were bandpass filtered at 0.5 and 45 Hz with the zero-phase shift FIR filter and then were segmented into 10 s epochs (approximately 60 epochs). Relative to the mean voltage of the entire epoch, every epoch was baseline-corrected at each time point. To select the abnormal channels, the normalized variance of each channel was calculated and then manually interpolated with the surrounding electrodes. We used independent component analysis (ICA) to reject other artifacts (EMG, EOG, and ECG; Jung et al., 2000). The artifacts were firstly selected with an automatic EEG artifact detector (Mognon et al., 2011). The EEG signals were re-referenced to the common average reference (CAR; Nunez et al., 2001). The first 40 clean epochs from each patient were preserved for further analysis.

The debiased weighted phase lag Index (dwPLI) was used for estimating spectral connectivity between pairs of channels. The dwPLI measuring of phase relationships is an estimator of scalp-level connectivity that is more robust and partially invariant to volume conduction in comparison to other estimators (Peraza et al., 2012). In the present study, the dwPLI measure was computed to estimate the functional connectivity between electrodes (Vinck et al., 2011) using the FieldTrip toolbox in MATLAB (Oostenveld et al., 2011).

First, we calculated cross-spectral decompositions from cleaned EEG datasets. Then, to achieve the dwPLI measure, the cross-spectrum between the spectral decompositions of every pair of channels was calculated. Furthermore, we restricted the analysis to the delta (1–4 Hz), theta (4–8 Hz), alpha (8–13 Hz), beta (13–30 Hz), and gamma (30–45 Hz) bands. Within each band, the dwPLI values across all epochs were averaged to represent the connectivity between channel pairs. For each patient’s dataset, the dwPLI values across all channel pairs were used to construct symmetric 32 × 32 × 4 dwPLI connectivity matrices for the delta, theta, alpha, beta, and gamma bands. Finally, each dwPLI matrix estimated above was proportionally thresholded. Values below this threshold were set to zero and other values were set to one, and 35% of the connectivity remained after thresholding. The threshold for each DOC patient was determined based on all four condition pools (T0, T1, T2, and T3) and was then applied to each condition separately. The binarized dwPLI connectivity matrices were used for further analysis. Due to the A1 and A2 electrodes being next to the face and vulnerable to noise, these two electrodes were excluded from the subsequent analysis. Finally, we obtained a symmetric 30 × 30 × 4 binarized dwPLI matrix for each band of each patient.

Two network measures, clustering coefficient (Watts and Strogatz, 1998) and global efficiency (Latora and Marchiori, 2001) were used to explore the properties that refer to the information processing in the human brain network. The brain-network measurement algorithms were implemented by the Brain Connectivity Toolbox (Rubinov and Sporns, 2010).

The clustering coefficient was computed as follows:

in which i, j refers to the electrode number, v_ij is 1 if there is suprathreshold connectivity between electrodes i and j and 0 otherwise, and n is the number of electrodes with suprathreshold connections to electrode i. The clustering coefficient reflects the integrity and interconnectedness of a smaller network. Namely, it reflects local information processing. To get the average clustering coefficient of a network, the clustering coefficient values over all 30 channels were averaged for each patient.

Global network efficiency is closely related to the characteristic path length. The formula to compute global efficiency is as follows:

where l_ij indicates the shortest path length from the node j to the node i. Global network efficiency could reflect the efficiency of global information exchange in brain networks.

To test whether the clustering coefficient and global efficiency were increased at T1, T2, and T3 compared to T0 in both VS and MCS groups, the single tail paired t-test was used. Ten repeated statistical tests (2 network properties × 5 EEG bands) were performed at each time point for the two groups. To reduce the false-positive rate for multiple comparisons, statistical tests were corrected with the BHFDR method (Benjamini and Hochberg, 1995). The similar statistical analysis method and corrected methods were then used to test whether the nodal dwPLI values were increased at T1, T2, and T3 compared to T0, only the repeated times were different, which was 30 in this test. When testing the significant changes of CRS-R scores, the beta/gamma-band average clustering coefficient, and global network efficiency between the two groups, a two-sample t-test, and BHFDR correction methods were used and the repeated times is five.

For the UWS patients, as shown in Figure 2A, the average clustering coefficient values in the delta band decreased with increasing sessions of HD-tDCS; however, this trend was not significant (p > 0.05, FDR corrected). The average clustering coefficient values in theta, alpha, beta, and gamma bands did not show any significant change at the T3 phase compared with those of the T0 phase. The results in Figure 3A were the same as the above results.

For the MCS patients, as shown in Figure 2B, the average clustering coefficient values in beta (p < 0.05, FDR corrected) and gamma (p < 0.05, FDR corrected) bands increased with sessions of HD-tDCS. The average clustering coefficient values in the delta band decreased at the T3 phase; however, this trend was not significant (p > 0.05, FDR corrected). The average clustering coefficient values in alpha and beta bands did not show a significant change at the T3 phase compared with those at the T0 phase. As shown in Figure 3B, the global efficiency values in the delta band decreased (p > 0.05, FDR corrected). In the beta and gamma bands, the global efficiency values increased with increasing sessions of HD-tDCS (p < 0.05, FDR corrected). In the delta and alpha bands, the global efficiency values did not show a significant change.

We analyzed the changes in CRS-R scores, average clustering coefficient, and global network efficiency at T1, T2, and T3 for UWS and MCS groups. As showed in Figure 4A, although the two groups both showed a score increase, the MCS patients had a higher score increase. The increased CRS-R scores in the MCS group significantly higher than in the UWS group at the T3 phase (p < 0.05, FDR corrected). In Figures 4B,C, the beta band changes of the average clustering coefficient and global network efficiency increased in T1 and T2 phases but decreased in T3 for the UWS group. For the MCS group, although it didn’t show a significant difference at the T3 phase between the two groups, the changes increased steadily with the HD-tDCS treatment. In Figures 4D,E, the gamma band changes of average clustering coefficient and global network efficiency were similar to above beta band results, but it showed significant difference (p < 0.05, FDR corrected) at T3 between the two groups.

The average clustering coefficient and global efficiency values could not estimate which brain regions function increased or decreased after 14 days of HD-tDCS therapy. Therefore, the average nodal ΔdwPLI topographs of UWS and MCS patients in the five frequency bands were plotted. To obtain the average nodal ΔdwPLI of each node, the symmetric 30 × 30 × 4 binarized dwPLI matrix was averaged by rows; then, we obtained the average nodal connection-strength matrix (30 × 4). Finally, the nodal ΔdwPLI matrix (30 × 1) was obtained by calculating the difference between the four phases (T1-T0, T2-T0, and T3-T0).

For the UWS patients, as shown in Figure 5A, there were no significantly increased nodal dwPLI values in the delta, theta, alpha, beta, or gamma bands between the T0 and T3 phase. Although the nodal dwPLI values increased in some regions in terms of theta, alpha, beta, and gamma bands for the T2 phase. But the results didn’t show statistical significance.

For the MCS patients, as shown in Figure 5B, the mean nodal dwPLI values decreased in the delta band. In the theta and alpha band, there were no significant changes in the mean nodal dwPLI values. In the beta band, the mean nodal dwPLI values increased in FC2, CP5, CP6, T8, P3, P7, P8, and POz (p < 0.05, FDR corrected). In the gamma band, the mean nodal dwPLI values increased in all electrodes (p < 0.05, FDR corrected), except Fp1, F3, F4, O1, O2, and Oz.

The changes in the nodal connection strength for MCS patients had a consistent increasing/decreasing trend during the HD-tDCS intervention. In contrast, for UWS patients, the changes in the node connection strength were not consistent. Additionally, the connection strength may have increased at the T2 phase, and then decreased at the T3 phase.

In this study, we found that there was an increasing trend in local and global information processing of resting-state functional brain networks during the 14 days of HD-tDCS interventions. Although UWS patients did not show any significant change from the T0 phase to the T3 phase, we still observed that information processing was increased. The electrophysiological findings were consistent with clinical statistical results (Table 1). For the MCS group, there were more patients exhibited significant improvements at T3 compared with the UWS group, and most of the patients with improved consciousness states had behaviorally improved. This finding was consistent with conclusions from previous studies. A previous study showed that the brain networks of DOC patients had reduced local and global efficiencies compared with those of healthy controls (Chennu et al., 2014). The increased local and global efficiencies after 14 days of HD-tDCS interventions suggested that HD-tDCS may be useful for treating DOC patients. Our recent study also demonstrated that long-lasting HD-tDCS over the precuneus is promising for the treatment of DOC patients (Guo et al., 2019). Rizkallah et al. (2019) reported that resting-state brain networks of DOC patients showed significantly increased local information processing when compared to that of controls. Similar results were reported by Chennu et al. (2017), in which they found that there was a linear trend with an increasing level of consciousness in local information processing of alpha-band brain networks.

The reasons why HD-tDCS interventions had better effects on the MCS group compared to the UWS group need to be further investigated. One possible reason is that the MCS patients have more residual functional connectivity of the default mode network (DMN; Naro et al., 2017). Such functional connectivity is an indicator of covert consciousness, and the HD-tDCS may indirectly modulate other brain networks through the residual DMN.

Surprisingly, our findings mainly revealed significant differences via HD-tDCS interventions in terms of beta and gamma bands. The beta frequency range is related to selective attention, large-scale neuronal integration, states of the motor system, and environmental stimuli (Gilbertson et al., 2005; Engel and Fries, 2010; Donner and Siegel, 2011; Bonfiglio et al., 2014); gamma oscillatory activity is correlated with the level of awareness (Naro et al., 2016; Koch et al., 2016; Siclari et al., 2017). In our present study, the increased network of local and global information processing in beta and gamma bands suggests that there was the amelioration of symptoms in the included MCS patients.

Another important finding was that the recovered functional connections in the beta and gamma bands were not only in the local brain area surrounding the stimulation position but were distributed across a more global brain area. This is not entirely inconsistent with our expected results; we originally hypothesized that the brain areas close to the stimulus position would be more likely to recover. Such results may be caused by the following two reasons. First, coma patients have been found to have a lower number of significant functional connections (Chennu et al., 2014; Malagurski et al., 2019), whereas MCS patients have been demonstrated to retain whole-brain functional connectivity (van den Heuvel et al., 2017; Malagurski et al., 2019). As such, stimulus information may be transmitted to other brain regions through these reserved connections. Additionally, although some specific brain regions or networks are specialized for performing specific tasks (He et al., 2018), integration of information cannot be done without whole-brain functional connectivity for complex cognition and behavior. Hence, the recovery of consciousness relies on strengthening of whole-brain functional connectivity.

Our study has two limitations. On the one hand, due to the patient’s family’s strong expectations for the recovery of consciousness, the sham stimulation group that lasted for 2 weeks did not recruit enough volunteers. So this study is not a random double-blind sham-controlled trial, thus the evidence obtained in the current study is not strong enough. For future HD-tDCS studies, adding a control group is necessary. On the other hand, although we have already recruited a lot of volunteers, the number is not big enough, more volunteers should be recruited for future study for making the results more convincing.