Criticality is considered a dynamic signature of healthy brain activity that can be measured on the short-term timescale with neural avalanches and long-term timescale with long-range temporal correlation (LRTC). It is unclear how the brain dynamics change in adult moyamoya disease (MMD). We used BOLD-fMRI for LRTC analysis from 16 hemorrhagic () and 34 ischemic () patients and 25 healthy controls. Afterwards, they were examined by EEG recordings in the eyes-closed (EC), eyes-open (EO), and working memory (WM) states. The EEG data of 11 and 13 patients and 21 healthy controls were in good quality for analysis. Regarding the 4 metrics of neural avalanches (e.g., size (), duration (), value, and branching parameter ()), both MMD subtypes exhibited subcritical states in the EC state. When switching to the WM state, remained inactive, while surpassed controls and became supercritical (). Regarding LRTC, the amplitude envelope in the EC state was more analogous to random noise in the MMD patients than in controls. During state transitions, LRTC decreased sharply in the controls but remained chaotic in the MMD individuals (). The spatial LRTC reduction distribution based on both EEG and fMRI in the EC state implied that, compared with controls, the two MMD subtypes might exhibit mutually independent but partially overlapping patterns. The regions showing decreased LRTC in both EEG and fMRI were the left supplemental motor area of and right pre-/postcentral gyrus and right inferior temporal gyrus of . This study not only sheds light on the decayed critical dynamics of MMD in both the resting and task states for the first time but also proposes several EEG and fMRI features to identify its two subtypes.

1. Introduction

Moyamoya disease (MMD) is a chronic cerebrovascular disease characterized by progressive stenosis and occlusion of the terminal portion of the bilateral internal carotid arteries and their main branches. The vascular pathology leads to widespread and continuous cerebral hypoperfusion and gradual formation of compensation from collaterals as a response [1]. The blood supply of neuron activities is often disrupted, followed by less regional neural activities and cognitive impairment [24]. The MMD is the only known chronic cerebrovascular disease with both hemorrhagic () and ischemic () clinical outcomes [5, 6]. Besides the general pathophysiology mentioned above in common, the two clinical subtypes are suspected to exhibit different but unknown processes, which are proved to be irrelevant with cerebral infarct or hematoma [3, 79]. However, their morphological and hemodynamic differences are difficult to be discerned and rarely reported [6, 9].

Healthy brain networks in the resting state are generally characterized by well-balanced excitatory and inhibitory synaptic activities [10]. These balanced brain states may sit in a dynamic state close to the criticality. Neural networks in a critical state have been found to be efficient for information transmission and other functions [11, 12]. Such critical dynamics are associated with a sustained neural activity that exhibits scale-free avalanche distribution [13]. Cognitive tasks may shift the brain activity from the critical point to a supercritical state by activating some brain regions and suppressing others [14]. The critical dynamics can be characterized quantitatively on two timescales: neural avalanches in the short term (milliseconds) and long-range temporal correlation (LRTC) in the long term (seconds to hours) [15]. During a neuronal avalanche, spontaneous activation of one neuronal group can trigger consecutive activations of other neuronal groups within just a few milliseconds, propagating cascading waves of activity. This phenomenon has been revealed by multiple neuroimaging modalities, such as electroencephalography (EEG) and fMRI in healthy subjects [1618]. Additionally, brains in disease states such as unconsciousness and schizophrenia exhibit faded critical dynamics [19, 20].

LRTC is another notable measurement to evaluate the critical dynamics of the neuronal system [21, 22]. Although neural oscillations present themselves with large variability in both frequency and amplitude, their fluctuations reflect a tendency toward self-organized criticality [23]. The oscillatory activity at any time is influenced by previous activities, and the LRTC is built up though local interactions. Such oscillations, reflecting short-term and long-term interactions, are observed to exist throughout the entire system. Recent studies have shown that healthy brains in the resting state are associated with a relatively large value of LRTC, while decreased LRTC has been reported in various diseased brain states such as a major depressive disorder and Alzheimer’s disease [24, 25].

Clinically, rapid and accurate differential diagnosis between acute ischemia and hemorrhage is crucial for early medical and interventional treatment but the optimal time window of treatment is often missed because of its high reliance on a CT or MR scan. The MMD is expected as the promising disease template to develop a more rapid and accurate differential diagnostic tool. According to one published fMRI study of ours, a dynamic measurement of entropy was proposed as an index of the critical dynamics to describe quantitatively the spatiotemporal changes of neural communication in adult MMD [26]. It found that critical dynamics faded not only in the diseased brain but also with disease progression. Therefore, this study was performed to examine two issues as the first of its kind. One is to explore the faded critical dynamics of MMD in both the resting and task states directly through a combination of EEG (high temporal resolution) and fMRI (high spatial resolution). The other is to investigate whether the two subtypes with a similar extent of cognitive impairment exhibit different neuronal dynamics and generate several features for future rapid and bedside identification of acute cerebral ischemia and hemorrhage.

2. Materials and Methods

2.1. Participants

This study was approved by the Institutional Review Board in our hospital and conducted in accordance with the Helsinki declaration. Informed consent was signed by all the subjects of this study. From March 2017 to August 2018, 50 adult patients with MMD (16 and 34 ) were recruited. The inclusion criteria were as follows: (1) Chinese nationality, right handed, and age of over 18 years; (2) diagnosis through digital subtraction angiography, with a Suzuki grade of III or IV [1]; (3) no evidence of infarct and hematoma larger than 8 mm in the maximum dimension on structural brain images [2, 27, 28]; (4) physical ability to undergo cognitive testing and EEG tasks; (5) no severe systemic or other cerebrovascular diseases; and (6) no medical history of neurosurgery. Twenty-five healthy subjects with no cerebrovascular or mental diseases were enrolled as controls. The Mini-Mental State Examination (MMSE) was adopted for the neuropsychological evaluation of global cognitive states.

2.2. Data Acquisition and Preprocessing

EEG data were acquired at a sampling rate of 1000 Hz in a sound-attenuated room by using a 64-channel actiCHamp Brain Products recording system (Brain Products GmbH Inc., Munich, Germany). The impedance of all channels was below 10 KΩ. The experimental paradigm was presented in E-Prime 2.0, and preprocessing was performed using the MATLAB R2017b software plug-in EEGLAB 14.0.0. Data were filtered to the frequency range of 0.5–100 Hz. Then, the interference at a frequency of 50 Hz was removed using a notch filter. Independent component analysis and the ADJUST toolbox were used to remove the components of eyeblink and cardio ballistic artifacts. EEG recording was started with a 5-minute eyes-closed (EC) resting state, a 5-minute eyes-open (EO) resting state, and then 30 trials of a delayed-response working memory (WM) task for around 20 minutes and ended with a procedure composed of a 5-minute EC and a 5-minute EO resting state to examine the consistency of the data recording quality (Supplementary Figure S2). Ultimately, the EEG data from 24 patients (11 and 13 ) and 21 healthy controls were selected for good quality.

All MRI data were collected using a 3.0 Tesla scanner (GE Healthcare, GE Asian Hub, Shanghai, China) with a 32-channel intraoperative head coil. The fMRI data were acquired using gradient echo-planar imaging with the following parameters: 3.2 mm slice thickness, 2000 ms repetition time, 30 ms echo time, 90° flip angle, field of view, matrix size, and voxel size. Each scan lasted for 8 min and collected 240 volumes. Preprocessing was performed using the SPM12 (https://www.fil.ion.ucl.ac.uk/spm) and DPARSF [29]. Volumes were corrected for slice timing and head motion, and subjects were excluded if their head movement exceeded 3 mm or 3 degrees. Here, we used the Friston 24-parameter model [30] to regress out the head motion effects from the realigned data (i.e., 6 head motion parameters, 6 head motion parameters on the time point before, and the 12 corresponding squared items) that were based on recent reports that higher-order models demonstrated benefits in removing head motion effects [31]. We also obtained the mean framewise displacement (FD) derived from Jenkinson’s relative root mean square algorithm [32] for each participant to assess the effects of head motion on group-level statistics. After spatial normalization (Montreal Neurological Institute space), resampling (3 mm isotropic voxels), and spatial smoothing (4 mm, full-width, half-maximum Gaussian kernel), volumes were preprocessed using the linear trend subtraction and temporal filtering (0.01–0.08 Hz). In addition, the signals from the cerebrospinal fluid and white matter were regressed out to reduce the effect of nonneuronal BOLD fluctuation [33].

2.3. Short-Term Timescales of Critical Dynamics: Neuronal Avalanches
2.3.1. Transfer of EEG Data from Amplitude to Mean Frequency

The time sequence of each channel was translated from amplitude to mean frequency by the short-term Fourier transform (Figures 1(a) and 1(b)). The mean frequency of each 1-sec epoch was calculated as with a 2 ms time shift in each channel, where and represent the frequency and power of the power spectrum of the signal, respectively. The average baseline of each channel was replaced with the group ensemble average of the channel to avoid noise interference caused by individual difference perturbations. Given that the EEG amplitude is mainly dominated by the low-frequency component of synchronized slow oscillations and some high-frequency activities are submerged in the slow-wave activities with high amplitudes, the EEG frequency may be a better sequence for reflecting the intensity of the neural activity rate.

2.3.2. Overthreshold Event Detection and Neuronal Avalanche Extraction

For each channel, the threshold for events was defined as the mean SD (standard deviation) of the EC state averaged overall controls. It is noted that for frequency events, 1 SD fluctuation above the mean frequency rate is strong enough to be a burst event (two and three SD were also tested, but the events were very sparse). The function “findpeaks” in MATLAB 2015b was used to find the peaks higher than the threshold (with a duration over 20 ms), which were defined as the events. Neuronal avalanches were defined as continuous time bins with at least one event in each channel, beginning and ending with event-free time bins. The time bin used was 16 ms (although the conclusion held true with other time bins from 4 to 20 ms). Binarized data were used for event detection.

2.3.3. Calculation of the Avalanche Size () and Duration ()

As shown in Figures 1(c) and 1(d), we first obtained the probability density function (PDF) of the avalanche size (number of events) and duration (number of time bins) and then fitted it to the power law distribution using maximum likelihood estimation to obtain the slope, denoted by the size () and duration () [34, 35]. The exponent was estimated by calculating , which means that the best-fitting exponent was calculated by maximizing the log-likelihood function. The higher absolute value of the power law exponent means the steeper fitted power law curve and indicates the higher probability of small size and duration avalanche [18]. Since the large size and duration avalanche means a rather large range and longer time of neuronal dynamics, the comparison of these parameters between the MMD and controls can be utilized to reflect the aberrant neuronal dynamics of MMD.

2.3.4. Calculation of the Value

The value was used to estimate how far the system was from criticality. This nonparametric measure quantifies the difference in the subjects’ cumulative density function (CDF) from the standard critical reference CDF with a power law exponent of 1.5 (Figure 1(e)) [11, 36]:

The value of this parameter at approximately 1 indicates a critical state of the system, while values above or below 1 indicate the super- or subcritical states, respectively.

2.3.5. Calculation of the Branching Parameter ()

The branching parameter is defined as the average number of subsequent events triggered by a single preceding event in an avalanche [34]. Thus, it was first calculated by averaging the ratio of the numbers of events in the second bin of an avalanche to the number of events in the first bin and then averaging all avalanches of a subject’s time series: where is the total number of cascades in the time series and is the number of the events in a certain time bin. For the single-bin cascades, the branching parameter is equal to 0. Therefore, the branching parameter of each cascade can vary over a wide range and the number of the cascades must be at least 700 in this study.

2.4. Long-Term Timescales of Critical Dynamics: LRTC
2.4.1. LRTC Based on EEG Data

Detrended fluctuation analysis was used to calculate the Hurst exponent of LRTC [37, 38]. The EEG data were bandpass filtered (finite impulse response filter) to the delta (1–4 Hz), theta (4–8 Hz), alpha (8–12 Hz), beta (12–25 Hz), and gamma (25–100 Hz) bands using a filter order of 2/minimum frequency of the band. Given a consecutive time series , the first step was to subtract the mean signal and obtain the cumulative sum of the signal: where denotes the time average. Then, the signal was divided into a set of time windows of the same length with 50% overlap. For each window in , the linear trend was removed using the least-squares method. Then, we computed the average root-mean-square fluctuation, , of all the time windows in with an identical length . The time window length was logarithmically spaced with a lower bound of 4 samples and an upper bound of 10% of the signal length [38]. According to our timescale in this work, we used 29 different window lengths ranging from 100 samples to 10000 samples (0.1 ~ 10 s). The Hurst exponent was defined as the coefficient of the linear regression of the sequence plotted on a log-log scale. For a scale-free signal , in the case of , the time series was correlated. When ranged from 0.5 to 1, the signal was considered to exhibit positive autocorrelation, while indicated an uncorrelated signal.

2.4.2. LRTC Based on fMRI Data

The voxel-wise Hurst exponent was adopted to evaluate the LRTC of fMRI data through the classical rescaled range (R/S) analysis [39]. The BOLD signal time series was divided into multiple shorter time series (length, ). The number of subtime series was , and was the length of the full-time series. For the subtime series of length , , where , we calculated the rescaled range :

The is the difference between the extremes of subtime series . is the cumulative deviation from to of the subtime series: where . is the difference between the maximum and minimum of and minimum:

The is the standard deviation of subtime series :

Calculating all the subtime series, we obtained the sequence corresponding to the sequence , where was all the possible values of the subtime series length.

2.5. Statistical Analysis

For EEG data analysis, one-way ANOVA was used and all the pairwise comparisons were assessed using Tukey-Kramer’s multiple comparison method. For fMRI data analysis, one-sample -tests were first performed on individual Hurst exponent maps for each group in a voxel-wise manner and to explore the within-group LRTC patterns using Rest version 1.8 ( for voxel-wise and cluster-wise , corrected through Gaussian random field theory) [40]. Afterwards, the averaged whole-brain Hurst exponent values were generated and compared among the three groups using the one-way ANOVA. Then, two-sample -tests were performed between each pair of subgroups on individual Hurst exponent maps ( for voxel-wise and cluster-wise ).

3. Results

3.1. Clinical Information

Table 1 shows the clinical information of the involved subjects. The demographic differences among the three groups were not significant (). However, the mean MMSE score of the healthy controls was significantly higher than that of the patient groups (). Among the subjects who were selected for the EEG analysis, differences of the demographics and cognitive and task performances among the three groups were not significant ().

3.2. Faded Critical Dynamics of Neuronal Avalanches
3.2.1. Cascade Size () and Duration ()

Figure 2 indicates that the three groups exhibited significant differences in both (EC, , ; EO, , ; and WM, , ) and (EC, , ; EO, , ; and WM, , ). In the subgroup analysis, and exhibited significant differences in in the EO ( ) and WM ( , ) states and in in the WM state ( , ). The control and groups showed significant differences in both ( , ) and ( , ) in the EC state. Additionally, the control and groups exhibited significant differences in in both the EC (, , ) and EO (, ) states and in in the EO state ( , ).

3.2.2. Deviations from Criticality with κ and the Branching Parameter of

The three groups exhibited significant differences in both (EC, , ; EO, , ; and WM, , ) and (EC, , ; EO, , ; and WM, , ) in all three states. and showed significant differences in in the EO ( ) and WM ( , ) states and in in the WM state ( , ). The control and groups showed significant differences in both ( , ) and ( , ) in the EC state. Additionally, the control and showed significant differences in in both the EC ( , ) and EO ( ) states and in in both the EC ( , ) and EO ( , ) states. The results kept the same tendency when utilizing a different parameter to calculate.

3.2.3. State Transition

The three groups exhibited a similar transition tendency between different states in all parameters (Figure 2(d)). For example, considering , both subtypes of MMD presented with larger values than controls in the EC state, implying a subcritical state. Similarly, MMD exhibited less activity than controls when switching to the EO state. Interestingly, when changing into the WM state, the group remained inactive but the group exhibited supercritical dynamics that were even more active than those of the controls.

3.3. Faded Critical Dynamics of LRTC
3.3.1. LRTC Based on EEG Data

The alpha band signal was taken for Hurst exponent analysis. Figure 3(a) indicates that the amplitude envelopes of the MMD groups were more analogous to the random noise than the control amplitude envelope was (the EC state, with a central-parietal channel as an example). The averaged Hurst exponents of the and groups were approximately 0.63 and 0.70, respectively, while that of the controls was 0.87 (Figure 3(b)).

The Hurst exponents of all the recording channels were calculated and mapped in Figure 3(c). Then, channels with significant differences in Hurst exponents among the three groups were mapped (Figure 3(d), Supplementary Table S1, ANOVA, ). Afterwards, the mean Hurst exponent value of the generated channels was calculated and compared among the three groups (Figure 3(e)). In the EC state, the control group exhibited significantly higher values in these channels than either (, ) or (, ). In the EO state, more channels presented a significant group difference and the showed a significantly higher value than the controls in these channels (). Furthermore, even more channels showed significant group differences in Hurst exponents in the WM state and the controls showed significantly lower values than the group (, ) or the group (, ).

When switching from the EC to EO state, 18 channels had significant Hurst exponent changes in the compared to the control group (Figures 4(a), ). Interestingly, the Hurst exponents of these 18 channels all decreased in the controls, while the group exhibited a chaotic pattern of change in these channels. In comparison, only 4 channels showed significant changes in the group compared to controls (Figure 4(b), ). Similarly, the Hurst exponents of these 4 channels all decreased in the controls, while the also showed a chaotic changing pattern in these channels. When switching from the EC to WM state, 27 channels in the and 28 channels in the showed significant changes in their Hurst exponents compared to those in the controls (Figures 4(c) and 4(d), ). Similarly, all these channels exhibited a large decline in Hurst exponents in the controls, while presenting with no regular pattern of change in either or .

3.3.2. LRTC Based on fMRI Data

The Hurst exponent patterns of the three groups are presented in Figures 5(a)5(c). Visual inspection indicated that in all three groups, the bilateral orbital frontal gyrus (OFG) and left precuneus (PCu) exhibited high values, while the bilateral fusiform gyrus (FFG), left inferior temporal gyrus(ITG), and left insular gyri (IG) showed low values. In addition, the bilateral supplemental motor area (SMA), precentral gyrus (PreCG), and postcentral (PoCG) gyrus of the controls; the left medial superior frontal gyrus (SFGmed) and right PCu of patients; and the bilateral SFGmed, dorsolateral prefrontal gyrus (DLPFC), SMA, and right PCu of patients all exhibited high values. In addition, the bilateral caudate nucleus (CN), hippocampus (HIP), parahippocampal gyrus (PHG), and thalamus (THA) of both controls and showed low values. Similar to the EEG results in the EC state (Figure 3(e)), the controls exhibited the highest mean Hurst exponent value among the three groups (Figure 5(d), ANOVA, , ), showing as controls .

Compared with controls, significant decreases in the Hurst exponent in the group were found in the left SMA, left DLPFC, left PCu, left superior parietal gyrus (SPG), and left middle occipital gyrus (MOG). Afterwards, significant decreases in the group compared to controls were found in the bilateral DLPFC, left SMA, right PreCG, right PoCG, and right ITG. No regions in MMD exhibited significantly higher Hurst exponent values than the controls. In addition, the left PreCG exhibited a significantly lower value in than in but no region showed a significant difference in the opposite direction (Figure 6, Table 2).

Since head micromovements could introduce artefactual interindividual differences in resting-state fMRI metrics [41, 42], we also measured the difference of head motion among the three groups. Although the exhibited the highest FD (Welch’s ANOVA, , Supplementary Figure S5), the mean Hurst exponent was not correlated significantly with mean FD across participants (, Supplementary Figure S6). The impact of head motion on the group-level differences was then tested by adding mean FD as a covariate, and results confirmed that head motion was not responsible for the differences among the three groups (Supplementary Figure S7 and Table S2).

3.3.3. LRTC Colocalization Patterns Based on EEG and fMRI Data in the EC State

The EEG channel placement was projected onto the cortical surface and converted to the Talairach Stereotactic System based on a published Brodmann’s area (BA) atlas [43]. Afterwards, we assessed the colocalization of regions with LRTC abnormities based on EEG data and fMRI data in the EC state. Compared to controls, regions with a significant Hurst exponent decrease in the group based on overlapping data were found in the left SMA (BA6, Figure 7(a)). In addition, regions with a significant decrease in the group based on overlapping data compared to controls were found in the right PreCG and PoCG (BA4 and BA6) and right ITG (BA19, Figure 7(b)).

4. Discussion

The criticality theory provides a novel insight into the neuronal dynamics underlying brain disorders. This study was the first to apply multiscale critical dynamics analysis to examine multimodal dynamical features in two moyamoya subtypes as compared to healthy controls. The neuronal avalanches on both fast and slow timescales were analyzed during rest and task performance, and several critical EEG features were derived. Both hemorrhagic and ischemic MMD exhibited particularly low EEG frequency activity and distinct subcritical dynamics, which could be distinguished easily from those of healthy controls. In addition, the decreased long-term correlations revealed in both high temporal (EEG) and spatial (fMRI) resolution were observed to reflect distinct neurophysiological processes associated with abnormal vascular network patterns in hemorrhagic and ischemic brains. Besides, this study provided clues for further rapid differential diagnosis between acute stroke and hemorrhage at the very early phase by use of EEG instead of CT and MR, which had greater advantages of rapidness, convenience, low cost, and radiation safe. Undoubtedly, time is the brain in treatment of acute ischemic stroke [44].

Previous investigations have suggested that the healthy brain in the resting state is usually characterized by well-balanced excitatory and inhibitory synaptic activities. These balanced levels of excitation and inhibition drive irregular spontaneous firing activities that exhibit scale-free avalanche distributions in the brain. Such a scale-free state can be effectively described by criticality [13]. Any input stimulus could effectively drive the brain into a supercritical state with additional excitatory activity, while relaxed low-signal states such as sleep can slow down the activity of the brain and shift it into a subcritical state. We noted that both subtypes of MMD exhibited subcritical states in the EC state and these suppressed dynamics prevented adaptive switching of brain function from introspective to extrospective states. This phenomenon might result from serious metabolic decrease and low neural activity rates caused by the chronic steno-occlusive angiopathy of MMD.

When switched to the EO state, both healthy controls and MMD presented with more neural activity and the MMD group remained less active than the controls. However, the ischemic moyamoya brains demonstrate stronger neural activities than hemorrhagic ones in the EO state. When switched to the WM state, all three groups exhibited more neural activities than that in the EO state, as was expected. Interestingly, the ischemic subtype surpassed the controls, while the hemorrhagic subtype still remained the least active. All parameters of neuronal avalanches exhibited similar results and were mutually verified. For healthy subjects, this phenomenon is reasonable because working memory is a behavioral state and requires effective neuronal activity to accomplish tasks [45]. However, we note that patients with ischemic MMD need to maintain a supercritical status to achieve similar scores to patients with hemorrhagic MMD and healthy subjects in certain WM tasks. In addition, these EEG features of neuronal avalanches during tasks may provide valuable clues for understanding the different neurophysiological processes of and .

For healthy brains in the EC state, spatially distributed neuronal activity may oscillate in phase with each other and result in high LRTC values [4648]. Compared with controls, the MMD group exhibited a reduced LRTC value close to 0.5 in the EC state, implying less correlated and more random brain activity. When subjects switched to EO and WM states, however, the LRTC value decreased sharply in the controls but both the and groups remained chaotic. The controls exhibited the lowest LRTC value in both the EO and WM states. In addition, spatial patterns of LRTC differences among the three groups were mapped based on EEG channels in both resting and WM states and under state transition. The healthy brain breaks the overall non-task-related dynamic balance into a set of subnetworks, such as the default-mode network, to deal with input signals effectively during state transition [49, 50]. However, the abnormal LRTC values and distribution of MMD in this study imply a completely different breakdown process for long-range neuronal fluctuations.

To further locate the regions with significant LRTC abnormities in MMD, we also examined BOLD fluctuations on fMRI due to the high spatial resolution of this modality. The results indicate that in the EC state, the patterns of the LRTC decreases in the two moyamoya subtypes are mutually independent but overlap in the left DLPFC of the executive control network and the left SMA of the salience network. Nevertheless, all regions of these patterns are key nodes involved in planning or direct control of movement, language, and visual information [51]. Referring to the pathophysiological nature of MMD, chronic stenosis/occlusion of the anterior circulation (bilateral internal carotid arteries and their main branches) is often followed by the collaterals from bilateral external carotid arteries and posterior circulation (bilateral vertebrobasilar arteries). Thus, the mismatch of the anterior circulation degradation and collateral development often results in a seemingly random and individualized cerebral hypoperfusion [7, 27, 52]. However, previous fMRI studies of MMD all revealed that patterns of functional deterioration are not random and key nodes of brain network such as the DLPFC, left SMA, are always involved [3, 4, 26, 28, 53]. Thus, this paper provides a crucial evidence that to output a similar extent of cognitive impairment, the neurophysiological processes of the two moyamoya subtypes may be mutually independent but overlap in some key nodes of the brain networks. Furthermore, we wondered whether there were potential links between spatial delay and temporal decay of neuronal oscillations in MMD and we attempted to trace their anatomical basis through both EEG and fMRI in the EC state. The identified regions are believed to play key roles in the neurophysiological processes of cognitive impairment in both ischemic and hemorrhagic MMD.

This study has several limitations that must be addressed. First, the EEG and fMRI data were not acquired at the same time. In order to generate a more stable and reliable result, the simultaneous EEG-fMRI technology should be used in future studies. Second, the study is based on a small sample size because completing EEG tasks is difficult for moyamoya patients with executive dysfunction. More patients are in need not only to increase the statistical power but to involve patients with Suzuki grading I–II and V–VI so as to obtain more knowledge of disease progression. Nevertheless, this study is the first of its kind to characterize the variability of brain dynamics in MMD on both short-term and long-term timescales and to show different neurophysiological features of its hemorrhagic and ischemic subtypes.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare no competing financial interests.

Authors’ Contributions

Y Mao, YX Gu, and YG Yu provided the idea of the article. Y Lei, YZ Li, and LC Yu wrote the article. LZ Xu, GX Zheng, and XY Qi analyzed the data. X Zhang, L Chen, and W Zhang collected the data. Yu Lei, Yuzhu Li, and Lianchun Yu contributed equally to this work.


This work was supported by the National Natural Science Foundation of China (grant numbers 81801155, 81771237, and 81761128011), the Shanghai Science and Technology Committee (grant numbers 18511102800 and 16JC1420100), the Shanghai Health and Family Planning Commission (grant number 2017BR022), the Shanghai Municipal Science and Technology Commission 795 Major Project (grant number 2018SHZDZX01), the Program for Professor of Special Appointment (Eastern Scholar) (grant number SHH1140004), and the Fundamental Research Funds for the Central Universities (grant number lzujbky-2015-119).

Supplementary Materials

Figure S1: the schematic instruction of the process of EEG data. EEG data was preprocessed and then transferred to mean frequency data to extract the cascade through event detection. The cascade size and duration were qualified to provide critical information used for clinical diagnosis. The artifact-free data was also bandpass filtered to measure the long-range temporal correlation altered in the patients. Figure S2: the EEG recording procedure. The EEG recording was started with a 5-minute eyes-closed (EC), a 5-minute eyes-open (EO) resting state, and then 30 trials of a delayed-response working memory (WM) task for around 20 minutes and ended with a procedure composed of a 5-minute EC and a 5-minute EO resting state to examine the consistency of the data recording quality. Figure S3: abnormal functional network connectivity mode of the three groups (alpha band). WM: working memory. Figure S4: Difference of dynamic (left) and static (right) mean frequencies with nonoverlapping 10 seconds’ epoch of EEG data among the three groups of hemorrhagic MMD, ischemic MMD, and healthy control in the resting and working memory states. Figure S5: differences of head motion among the three groups. The mean framewise displacement (FD) derived from Jenkinson’s relative root mean square (RMS) algorithm [32] was used to measure the extent of head motion. The groups exhibited the highest FD (Welch’s ANOVA, ): . Figure S6: the correlation between head motion and Hurst exponent. Red lines indicate the linear fit. Figure S7: between-group difference of LRTC maps with the consideration of head motion effects. A, the -statistical difference maps of LRTC between the healthy and hemorrhagic MMD groups; the positive value is termed as the stronger effect in controls. The statistical threshold is set as for a voxel-wise threshold of 0.001; the minimum cluster size for a voxel threshold of 0.001 and cluster threshold of 0.05 is 24 voxels. B, the -statistical difference maps of LRTC between the healthy and ischemic MMD groups; the positive value means a stronger effect in controls. The statistical threshold is set at for a voxel-wise threshold of 0.001; the minimum cluster size for a voxel threshold of 0.001 and cluster threshold of 0.05 is 21 voxels. C, the -statistical difference maps of LRTC between the ischemic and hemorrhagic MMD groups; the positive value indicates a stronger effect in ischemic MMD. The statistical threshold is set at for a voxel-wise threshold of 0.003; the minimum cluster size for a voxel threshold of 0.001 and a cluster threshold of 0.05 is 29 voxels. R: right; L: left. Table S1: channels with significant difference of Hurst exponent among the three groups in all three states. Table S2: comparison of brain regions with significant difference in LRTC among groups with and without framewise displacement covariate. (Supplementary Materials)