Introduction

Parkinson’s disease (PD) is a common neurodegenerative disorder among middle-aged and elderly populations, primarily characterized by the progressive degeneration of dopaminergic neurons in the substantia nigra of the midbrain1,2,3. Patients with PD mainly exhibit motor symptoms such as resting tremor, freezing gait, and bradykinesia, and may also experience non-motor symptoms such as cognitive impairment and sleep disturbances4,5,6,7. The primary treatment for PD patients involves the use of dopaminergic drugs, but long-term medication can lead to side effects such as dyskinesia and reduced drug efficacy8. Therefore, there is an urgent need to find new adjunctive therapeutic methods to delay the progression of the disease.

Repetitive transcranial magnetic stimulation (rTMS) is a non-invasive neuroregulatory technique that can modulate cortical excitability9. It is characterized by its non-invasiveness, ease of operation, and high safety. Relevant studies have shown that high-frequency (≥5 Hz) rTMS (HF-rTMS) can enhance cortical excitability, while low-frequency (≤1 Hz) rTMS (LF-rTMS) has an inhibitory effect10,11,12. rTMS treatment can serve as a non-pharmacological adjunctive therapy to improve the clinical symptoms of patients with neurological and psychiatric diseases such as Parkinson’s disease, schizophrenia, and depression13,14.

In recent years, a patterned form of rTMS known as theta burst stimulation (TBS) has emerged. It is a variant of rTMS that simulates neural oscillations for pulse delivery15,16. Intermittent TBS (iTBS) is an excitatory stimulus similar to HF-rTMS. For example, delivering 600 pulses, 10Hz-rTMS would take approximately 500 s, while iTBS could complete it in 192 s, significantly improving clinical practicality. Studies have shown that iTBS can produce similar or stronger excitatory effects compared to HF-rTMS17,18,19. However, the results are not consistent. Currently, the therapeutic effects of these two stimulation modes on PD patients are mostly evaluated using clinical scales and behavioral tasks, and further research is needed on the impact of these two stimulation modes on brain function in PD patients. Changes in brain function are crucial for assessing the therapeutic effects of rTMS on PD patients.

The human brain is a complex network system composed of numerous neural cell connections, and the information transmission process between neural cells generates brain electrical signals. EEG can record changes in cortical neuroelectrophysiological activity, characterized by its non-invasiveness and high temporal resolution20. Brain network analysis, as an emerging EEG signal analysis method, can accurately reflect changes in brain functional connectivity. Degeneration of brain function in PD patients can lead to changes in brain functional connectivity21. Traditional dynamic brain network studies often use fixed sliding time windows, and the choice of window size can affect the research results22. Fixed-size windows also cannot accurately reflect various activity stages of the brain. Brain network analysis mostly focuses on studying changes in the overall network of the brain and the functional connectivity attributes between specific brain regions, which can miss local sensitivity and regular information between brain regions. Brain activity always maintains a briefly stable state for about 60–120 milliseconds before rapidly transitioning to another state. This briefly stable state is vividly called a “thought atom”, or functional microstate23,24. Different brain functional networks determine different types of microstates, and currently, four types of microstate categories (A–D) have been repeatedly identified in most research results. These four types of microstate topographies can explain most of the global changes25,26,27. Previous studies have also found that microstate categories are related to known resting-state networks defined by functional magnetic resonance imaging (fMRI), with microstate A related to speech processing networks, B related to visual networks, C related to salience networks, and D related to attention networks28. Guan et al.29 found that dynamic brain networks based on EEG microstates can better study the neural mechanisms of cross-task psychological load. Chu et al.30 discovered that the spatiotemporal variability of dynamic brain functional networks based on EEG microstates can reflect the pathological changes in the brains of PD patients and can classify PD patients and healthy subjects with a high classification rate. Therefore, we have reason to speculate that the spatiotemporal variability of dynamic brain functional networks based on EEG microstates can fully reflect the impact of rTMS on brain functional connectivity in PD patients.

In this paper, we construct brain functional networks based on the time windows of EEG microstates and calculate their spatiotemporal variability, comparing the differences in spatiotemporal variability between PD patients and healthy subjects. Then, we divide PD patients into two groups, treating them with iTBS and 10Hz-rTMS, respectively, and compare the changes in clinical assessment results and spatiotemporal variability before and after treatment. In addition, we explore the correlation between clinical assessment results and spatiotemporal variability in PD patients.

Results

Demographic characteristics and clinical assessment results

As shown in Table 1, there were no significant differences in gender, age, and education level among the subjects (P > 0.05). There were no statistically significant differences between the iTBS and 10Hz-rTMS groups in terms of illness duration or daily levodopa equivalent dose (LEDD) (P > 0.05).

Table 1 Demographic data
Full size table

Table 2 presents the clinical assessment results of PD patients. Compared to baseline, both the iTBS and 10Hz-rTMS groups showed significant improvements in clinical scale scores (UPDRS-III and MoCA) and behavioral task performance (SSRT and SS-180) after treatment (P < 0.05).

Table 2 Clinical assessment results of PD patients
Full size table

Temporal variability of dynamic brain functional networks based on EEG microstates

We calculated the temporal variability of dynamic brain functional networks for EEG microstates in different groups. Figure 1a–c displays the topological distribution and differences in temporal variability within each microstate category among three types of group comparisons (PD vs HC, Pre-10Hz-rTMS vs Post-10Hz-rTMS, and Pre-iTBS vs Post-iTBS). In the microstate A network, PD patients exhibited significantly lower temporal variability in the prefrontal cortex extending to the right temporal lobe compared to HC. After 10Hz-rTMS treatment, the temporal variability in the prefrontal cortex significantly increased, and after iTBS treatment, the temporal variability in the frontal lobe significantly increased. Regarding the microstate B network, there were smaller differences in temporal variability between PD patients and HC. After 10Hz-rTMS treatment, the temporal variability in the left frontal and right occipital lobes significantly decreased, and after iTBS treatment, the temporal variability in the frontal and occipital lobes significantly decreased. The temporal variability of the microstate C network in PD patients was generally lower than that of HC, with significant differences extending from the left frontal lobe to the right parietal and occipital lobes. After both iTBS and 10Hz-rTMS treatments, the temporal variability in the frontal, parietal, and occipital lobes significantly increased. In contrast to the microstate C network, the temporal variability under the microstate D network in PD patients was generally higher than that of HC, with significant differences throughout the frontal, parietal, and occipital lobes. After both iTBS and 10Hz-rTMS treatments, the temporal variability in the frontal, parietal, and occipital lobes significantly decreased.

Fig. 1: Topological distribution of temporal variability.
figure 1

The distribution of dots represents the ___location of channels, with small dots indicating P < 0.05 and large dots indicating P < 0.01. a Topological distribution of temporal variability between PD patients and HC, where red dots indicate significantly higher temporal variability in PD patients compared to HC, and blue dots indicate significantly lower temporal variability in PD patients compared to HC (P < 0.05 and r > 0.353). b Topological distribution of temporal variability before and after treatment in the 10Hz-rTMS group, where red dots indicate significantly lower temporal variability after treatment compared to before treatment, and blue dots indicate significantly higher temporal variability after treatment compared to before treatment (P < 0.05 and r > 0.3). c Topological distribution of temporal variability before and after treatment in the iTBS group, where red dots indicate significantly lower temporal variability after treatment compared to before treatment, and blue dots indicate significantly higher temporal variability after treatment compared to before treatment (P < 0.05 and r > 0.3).

Full size image

Spatial variability of dynamic brain functional networks based on EEG microstates

Figure 2a–c presents the topological distribution and differences in spatial variability within each microstate category among three types of group comparisons (PD vs HC, Pre-10Hz-rTMS vs Post-10Hz-rTMS, and Pre-iTBS vs Post-iTBS). Concerning the microstate A network, we observed that PD patients exhibited significantly higher spatial variability in the right frontal and right temporal lobes compared to HC. After 10Hz-rTMS treatment, the spatial variability in the right prefrontal region significantly decreased, and after iTBS treatment, the spatial variability in the right frontal lobe significantly decreased. In the microstate B network, PD patients showed significantly higher spatial variability in the left frontal and right occipital lobes compared to HC. After iTBS treatment, the spatial variability in the prefrontal and right parietal lobes significantly decreased. The differences between PD patients and HC under the microstate C network mainly existed in the parietal and bilateral temporal lobes, where PD patients showed significantly lower spatial variability compared to HC. After both treatment modalities, the spatial variability in the parietal and bilateral temporal lobes of PD patients significantly increased. In the microstate D network, PD patients generally exhibited higher spatial variability compared to HC, with significant differences in the prefrontal, parietal, and bilateral temporal lobes. After 10Hz-rTMS treatment, the spatial variability in the parietal and bilateral temporal lobes significantly decreased, and after iTBS treatment, the spatial variability from the left frontal to the right parietal and bilateral temporal lobes significantly decreased.

Fig. 2: Topological distribution of spatial variability.
figure 2

The distribution of dots represents the ___location of channels, with small dots indicating P < 0.05 and large dots indicating P < 0.01. a Topological distribution of spatial variability between PD patients and HC, where red dots indicate significantly higher spatial variability in PD patients compared to HC, and blue dots indicate significantly lower spatial variability in PD patients compared to HC (P < 0.05 and r > 0.353). b Topological distribution of spatial variability before and after treatment in the 10Hz-rTMS group, where red dots indicate significantly lower spatial variability after treatment compared to before treatment, and blue dots indicate significantly higher spatial variability after treatment compared to before treatment (P < 0.05 and r > 0.3). c Topological distribution of spatial variability before and after treatment in the iTBS group, where red dots indicate significantly lower spatial variability after treatment compared to before treatment, and blue dots indicate significantly higher spatial variability after treatment compared to before treatment (P < 0.05 and r > 0.3).

Full size image

Correlation between spatiotemporal variability and clinical assessment results

Figure 3 shows the correlation between the spatiotemporal variability of PD patients before treatment and clinical assessment results. The temporal variability of the occipital lobe microstate B network is significantly negatively correlated with MoCA (P < 0.05), meaning that the better the cognitive function of PD patients, the lower the temporal variability of the occipital lobe microstate B network. Although there is no significant difference, we found that the temporal variability of the occipital lobe microstate B network in PD patients significantly decreased after both treatment modalities. Regarding the microstate C network, the temporal variability of the right parietal and right temporal lobes is significantly negatively correlated with SSRT (P < 0.05), indicating that the better the inhibitory control ability of PD patients, the higher the temporal variability of these two regions. However, we found that the temporal variability of these two regions in PD patients is significantly lower than that of HC. In the microstate D network, the correlation between the spatiotemporal variability of related brain regions and clinical assessment results is more obvious. The temporal variability of the parietal lobe region is significantly positively correlated with UPDRS-III (P < 0.05), meaning that the more normal the motor function of PD patients, the lower the temporal variability of the parietal lobe region. However, we found that the temporal variability of the parietal lobe region in PD patients is significantly higher than that of HC. Additionally, we found that the spatial variability of the left frontal lobe is significantly positively correlated with SSRT (P < 0.05), indicating that the better the inhibitory control ability of PD patients, the lower the spatial variability of the left frontal lobe. However, the results show that the spatial variability of the left frontal lobe in PD patients is significantly higher than that of HC. The spatial variability of the occipital lobe region is significantly positively correlated with the completion time of the SS-180 task (P < 0.05), meaning that the better the motor function of PD patients, the lower the spatial variability of the occipital lobe region. However, the spatial variability of the occipital lobe region in PD patients is significantly higher than that of HC. We also found that the spatial variability of the right temporal lobe is significantly negatively correlated with MoCA (P < 0.05), indicating that the better the cognitive function of PD patients, the lower the spatial variability of the right temporal lobe. However, the spatial variability of the right temporal lobe region in PD patients is significantly higher than that of HC. There are no significant correlations between brain regions in other states and clinical scales and behavioral tasks (P > 0.05).

Fig. 3
figure 3

Correlation between spatiotemporal variability and clinical assessment results in PD patients before treatment, the distribution of purple dots represents the ___location of channels.

Full size image

Discussion

rTMS is a non-invasive and effective neuroregulatory technique that plays a supportive role in the treatment of neuropsychiatric diseases. However, the differences in intervention effects between the two different stimulation modes (iTBS and 10Hz-rTMS) on PD remain unknown. Compared to traditional methods of constructing brain functional networks, constructing dynamic brain functional networks based on microstate categories can more accurately represent changes in various stages of the brain. Analyzing the spatiotemporal variability of brain regions also avoids losing information on local sensitivity and inter-regional regularity. This study used dynamic brain functional networks constructed based on EEG microstates to investigate changes in spatiotemporal variability in PD patients compared to age-matched HC, before and after iTBS treatment, and before and after 10Hz-rTMS treatment. This was done to assess abnormalities in brain functional connectivity in PD patients and the impact of different stimulation modes on brain functional connectivity in PD patients. Additionally, this paper explored the relationship between spatiotemporal variability and clinical assessment results in PD patients.

In the microstate A network, PD patients exhibited significantly lower temporal variability in the prefrontal and right temporal regions compared to HC, and significantly higher spatial variability in the right frontal and right temporal regions. Microstate A is associated with the brain’s speech processing network31, which is responsible for the storage and control of language information. The dorsolateral prefrontal cortex is involved in short-term storage and processing of information32, and the temporal lobe is related to memory and emotion33. Some PD patients experience non-motor symptoms such as memory decline and depression, and we speculate that the abnormal changes in spatiotemporal variability in the frontal and temporal regions reflect the decline in short-term language memory and processing functions in PD patients. After treatment, the temporal variability in the frontal regions of PD patients significantly increased, and the spatial variability in the right frontal region significantly decreased. Previous studies have shown that 10Hz-rTMS treatment in the supplementary motor area of PD patients can change the functional connectivity of the prefrontal region34. Therefore, we speculate that the changes in spatiotemporal variability in related brain regions after treatment reflect improvements in the patients’ short-term memory and storage processing functions. In the future, the sample size can be increased to study the relationship between the prefrontal region under the microstate A network and language memory and processing functions in PD patients by subtype.

Before treatment, there was no significant difference in temporal variability in the microstate B network between PD patients and healthy subjects, and differences in spatial variability were mainly distributed in the left frontal and right occipital regions. Microstate B is closely related to the visual network35, and the acquisition of visual information is an important part of cognitive processes. PD patients are mainly characterized by motor symptoms, and some patients’ non-motor symptoms are not obvious or relatively mild in the early stages of the disease, which may lead to no significant difference in temporal variability under the microstate B network. The temporal variability of the occipital region is significantly negatively correlated with the MoCA score, indicating that the higher the temporal variability of the occipital region, the worse the patient’s cognitive function. After treatment, the temporal variability in the left frontal and right occipital regions of patients in the 10Hz-rTMS group significantly decreased, and the iTBS group involved changes in a wider range of brain regions, including most of the frontal lobe and the entire occipital region. The occipital lobe is the main area of the brain for processing visual information36, and the higher spatial variability in the right occipital region of PD patients may reflect abnormalities in the related visual network. Related studies have shown that reduced dopamine secretion in the brain of PD patients causes cognitive decline37, and rTMS stimulation can increase dopamine secretion in the brain38. Therefore, we speculate that after treatment, the increase in dopamine levels in the brains of patients in both groups relieved the fatigue of the visual network and improved the stability of the visual network39,40. As a result, the overall cognitive level of patients within the group was improved, reflected in the significant reduction in temporal variability in the frontal and occipital regions of patients. Therefore, after treatment, the brain network connection pattern in the frontal and occipital regions under the microstate B network of PD patients became more stable, and the temporal variability in the occipital region under the microstate B network can reflect changes in cognitive function in PD patients.

The results of this study show that the temporal variability in the left frontal, right parietal, and right occipital regions under the microstate C network in PD patients is significantly lower than that in HC, and the spatial variability in the bilateral temporal and parietal regions is significantly lower than that in HC. Low spatiotemporal variability indicates that the functional connectivity patterns or synchronicity of the above-related brain regions in PD patients are stronger. We also found that the temporal variability in the right parietal and right occipital regions of PD patients is significantly negatively correlated with SSRT; that is, the lower the temporal variability in the right parietal and right occipital regions, the worse the SST performance. Previous studies have shown that microstate C is significantly related to the cognitive level of PD patients41, and the SST task includes motor inhibition and cognitive processes of recognition. The parietal region is involved in functions related to motor control42, and the occipital region is related to the visual network36. The acquisition of visual information in the SST task is a key step in completing the task, and patients need to stare at the screen and make judgments on changing information. Therefore, we speculate that the reduction in spatiotemporal variability in related brain regions under the microstate C network in PD patients is related to the decline in motor inhibition and complex cognitive functions. After treatment, the spatiotemporal variability in abnormal brain regions significantly increased. Combined studies of fMRI and EEG have shown that microstate C is closely related to brain regions associated with cognitive function and executive control43. iTBS stimulation of the M1 area can change the blood oxygenation level-dependent (BOLD) of brain regions related to movement44. Therefore, we speculate that after treatment, changes in BOLD in the brains of patients in both groups improved the patients’ motor control and cognitive functions, thereby increasing the spatiotemporal variability in related brain regions. The spatiotemporal variability of the microstate C network can reflect changes in motor and cognitive functions in PD patients.

In this study, before treatment, the temporal variability in the frontal, parietal, and occipital regions under the microstate D network in PD patients was significantly higher than that in HC, and the spatial variability in the prefrontal, parietal, and bilateral temporal regions was significantly higher than that in HC. This indicates that the stability of functional connectivity patterns in the frontal, parietal, and occipital regions under the microstate D network in PD patients is lower than that in HC, and the synchronicity of dynamic functional connections between the prefrontal, parietal, and bilateral temporal regions and other brain regions is also lower than that in HC. Regarding specific brain region correlations, we found that the spatiotemporal variability of some brain regions in PD patients is significantly correlated with clinical assessment results. The spatial variability in the left frontal region is significantly positively correlated with SSRT; that is, the greater the temporal variability in the frontal region of PD patients, the longer the SSRT, and the worse the inhibitory control ability. The spatial variability in the right temporal region is significantly negatively correlated with MoCA; that is, the greater the spatial variability in the temporal region of PD patients, the worse the cognitive function. Previous studies have found that cognitive impairment in PD and schizophrenia patients is related to reduced dopamine secretion in the brain, causing abnormalities in microstate D45,46. The frontal lobe is related to complex higher cognitive processes32, and the temporal lobe is responsible for memory-related cognitive processes33, so the increase in spatial variability in the left frontal and right temporal regions under the microstate D network may reflect the decline in complex higher cognitive functions in PD patients. Reduced dopamine can cause symptoms such as tremors and bradykinesia in PD patients, affecting the patients’ motor functions47. The temporal variability in the parietal region is significantly positively correlated with UPDRS-III; that is, the greater the temporal variability in the parietal region of PD patients, the higher the UPDRS-III score, and the worse the motor function. The spatial variability in the parietal region is significantly positively correlated with the SS-180 task; that is, the greater the spatial variability in the parietal region of PD patients, the longer the completion time of the SS-180 task, and the worse the motor function. The parietal region is mainly responsible for motor functions and maps the motor function network42, and the increase in temporal variability in the parietal region may be related to the decline in motor functions in PD patients. After treatment, the temporal variability in the frontal, parietal, and occipital regions of patients significantly decreased, the spatial variability in the bilateral temporal and parietal regions of patients in the 10Hz-rTMS group significantly decreased, and the spatial variability in the left frontal, bilateral temporal, and parietal regions of patients in the iTBS group significantly decreased. PD patients can significantly change the characteristic parameters of microstate D after receiving dopamine drug treatment or rTMS treatment48,49. Therefore, we speculate that the abnormal motor and cognitive functions of PD patients cause the spatiotemporal variability of some brain regions under the microstate D network to significantly increase. After treatment, the dopamine level in the brain of PD patients increased, reducing the spatiotemporal variability of related brain regions. Comparing the results after treatment in both groups, except for the bilateral temporal and parietal regions, the spatial variability in the left frontal region of patients in the iTBS group significantly decreased. Combined with previous correlation analysis, we speculate that iTBS can better improve the complex higher cognitive functions of PD patients. The spatiotemporal variability of the microstate D network can become a powerful indicator of changes in motor and cognitive functions in PD patients.

In summary, both stimulation modes can improve the clinical symptoms of PD patients. Compared to the 10Hz-rTMS mode, the iTBS mode can significantly reduce the temporal variability in the right frontal and parietal regions under the microstate B network. Under the microstate D network, it can significantly reduce the spatial variability in the left frontal region, and the spatial variability in the left frontal region is significantly correlated with SSRT in PD patients. The iTBS mode can better improve the brain functional connectivity of abnormal brain regions in PD patients.

This study has several limitations that warrant acknowledgment. First, the relatively small sample size (51 PD patients and 21 healthy controls) may limit the generalizability of our findings. Future research should recruit larger cohorts to validate the relationship between spatiotemporal variability in identified brain regions and clinical symptomatology in PD patients. Second, the integration of multimodal imaging modalities (e.g., fMRI) could provide more comprehensive spatiotemporal insights, potentially offering stronger evidence for the neural mechanisms underlying rTMS effects in PD. Finally, our analysis was confined to the short-term effects of iTBS and 10Hz-rTMS. Longitudinal follow-up studies are needed to assess the enduring therapeutic benefits of these interventions in PD patients over extended periods.

The results of this study indicate that, compared to the 10Hz-rTMS mode, the iTBS mode can achieve the same therapeutic effect, and there were no differences in clinical assessment results between the two groups of patients after treatment. In the microstate B and D networks, the iTBS mode improved the brain functional connectivity of a larger range of brain regions in PD patients. Therefore, the iTBS mode may achieve the same or better therapeutic effect in a shorter time for PD patients. The spatiotemporal variability of dynamic brain functional networks based on EEG microstates can reveal abnormalities in brain functional activity in PD patients and can reflect the effects of different modes of rTMS on neural regulation in PD patients. We speculate that the spatiotemporal variability of dynamic brain functional networks based on EEG microstates can also be used in the study of other neuropsychiatric diseases.

Methods

Subjects

From 2021 to 2024, 51 PD patients from the Neurosurgery Outpatient Department of Tianjin Huanhu Hospital and 21 healthy subjects of the same age group were enrolled. The inclusion criteria for PD patients were as follows: One. PD subjects had been diagnosed with Parkinson’s disease. Tow. They could complete the sit-to-stand walking action without the assistance of tools. Three. Cognitive function was basically normal, with no communication barriers and good compliance. Four. There was no history of other neuropsychiatric diseases. Five. There were no metal implants in the head. Six. There were no changes in the type and dosage of medication during the treatment period.

All subjects were right-handed, and patients and their families jointly signed the “Treatment Informed Consent Form.” This study was approved by the Biomedical Ethics Committee of Hebei University of Technology, and all procedures followed the “Helsinki Declaration.”

Research design

The experimental procedure is shown in Fig. 4a. During the baseline assessment phase, all PD patients underwent EEG recording during eyes-closed resting state, followed by clinical scale evaluations and behavioral task assessments conducted by professional clinicians. The clinical assessments included the Unified Parkinson’s Disease Rating Scale Part III (UPDRS-III) and Montreal Cognitive Assessment (MoCA), while behavioral tasks consisted of the Standing Start 180° Turn Test (SS-180) and Stop Signal Task (SST). PD patients were then randomly assigned to either iTBS or 10Hz-rTMS treatment groups, receiving daily sessions (10 sessions total over two weeks). Post-treatment reevaluation included repeat EEG recordings during eyes-closed resting state and repetition of clinical/behavioral assessments. Additionally, we collected eyes-closed resting-state EEG data from 21 healthy controls.

Fig. 4: Experimental workflow diagram.
figure 4

a Design and flow of participants through the study. b Schematic of SS-180. c Schematic of SST. d Schematic of 10Hz-rTMS Pulse Delivery Pattern. e Schematic of iTBS Pulse Delivery Pattern.

Full size image

Behavioral tasks

We further tested the motor function of PD patients using the SS-180, as shown in Fig. 4b. Patients sat on a chair, with an obstacle placed 4.2 m away from the chair. Upon hearing the “start” command, patients stood up, walked to the obstacle, bypassed it, and then returned to the original position to sit down. To ensure the reliability of the test results, the entire process was recorded with dual cameras, and the completion time was obtained by replaying the video. The average time of two test completions by the subjects was taken as the final time.

Response inhibition refers to the suppression of an already formed response impulse. The SST is one of the commonly used paradigms for studying response inhibition and is widely applied in research on response inhibition impairments in patients with PD, attention deficit hyperactivity disorder, and other conditions. The response inhibition task is illustrated in Fig. 4c. Patients sit in front of a computer screen, where a fixation point “+” first appears in the center of the screen. After a fixed time, a response signal appears, which is a white arrow pointing either left or right, instructing patients to respond quickly upon seeing the signal. If the arrow turns red after the response signal appears, it serves as a stop signal, requiring patients to suppress the urge to press the button and make no response. We use the Stop-signal reaction time (SSRT) as the performance measure for the SST, calculated as the difference between the Go task reaction time (Go RT) and the Stop-signal delay time (SSDT). To make the measurement results more accurate, we use a “tracking method” to set the SSDT: after a successful inhibition of the button press impulse in a particular SST, the SSDT will increase by 50 ms in the next SST, thereby increasing the difficulty of inhibition. Conversely, it will decrease by 50 ms. A complete set of SST tasks consists of 100 trials, including 26 stop tasks. After completing three sets of complete SSTs, the average value is taken to obtain the SSRT.

rTMS protocol

This study utilized the Rapid2 stimulator produced by Magstim Company in the UK, equipped with an “8”-shaped air-cooled coil. The stimulation target was the bilateral primary motor cortex (M1) of the patients, with the center of the coil closely adhering to the target and positioned at a 45° angle to the sagittal suture. The stimulation target was selected based on the ___location where at least five out of ten consecutive single pulses resulted in an electromyographic amplitude of 50 mV or greater. The stimulation intensity for both modes was set at 90% of the resting motor threshold. The 10Hz-rTMS stimulation protocol is shown in Fig. 4d: a pulse frequency of 10 Hz, lasting 4 s, with a 26-s interval, continuously stimulating for 30 cycles. Unilateral stimulation lasted 15 min, totaling 1200 pulses, with bilateral stimulation lasting a total of 30 min. The iTBS stimulation protocol is illustrated in Fig. 4e: an inter-burst frequency of 5 Hz, lasting 2 s, with an 8-s interval, and an intra-burst frequency of 50 Hz, meaning each burst contained three pulses. A single iTBS session lasted 192 s, totaling 600 pulses, with unilateral stimulation administered twice, each time followed by a 10-min rest interval, and bilateral stimulation totaling four times, ensuring that the total number of pulses delivered by both modes was the same. To minimize external environmental interference, the treatment room was kept quiet during sessions, patients wore earplugs throughout, and treatment was administered once daily for a total of 10 sessions over 2 weeks.

EEG data collection

We used the NeuSen W-series 32-channel wireless EEG acquisition system produced by Neuracle Company for EEG data collection. The electrode placement positions were referenced to the international standard 10–20 system, with bilateral mastoids selected as the reference electrode placement points. Before collection, the impedance of each electrode was reduced to below 10 KΩ, and the sampling frequency was set to 1000 Hz. Subjects were placed in a quiet and comfortable room, wore earplugs throughout the process to avoid external environmental interference. Before collection, subjects were reminded to maintain a closed-eye resting state, relax, and minimize unnecessary movements. Each collection session lasted 6 min.

EEG data preprocessing

EEG signals are inherently weak and contain various background noises, resulting in a low signal-to-noise ratio in raw EEG data. We preprocessed the EEG signals using the EEGLAB toolbox. The preprocessing pipeline included the following steps: (1) removal of bilateral mastoid reference electrodes and downsampling to 500 Hz; (2) manual rejection of segments with significant artifacts, with spherical interpolation applied to heavily artifact-contaminated electrodes; (3) average referencing followed by 0.5–80 Hz bandpass filtering; (4) 50 Hz notch filtering to eliminate powerline interference; (5) independent component analysis (ICA) for artifact and noise removal50; (6) 2–20 Hz bandpass filtering followed by extraction of a 4-minute EEG segment for microstate analysis.

Construction of dynamic brain functional networks based on EEG microstate windows

Figure 5a illustrates the construction process of dynamic brain functional networks based on EEG microstate windows. The first step involves calculating the global field power (GFP) of the preprocessed raw EEG signals51:

$$\mathrm{GFP}=\sqrt{\frac{\mathop{\sum }\limits_{i=1}^{N}{({u}_{i}-\overline{u})}^{2}}{N}}$$
(1)

Where N is the number of electrodes, ui is the voltage value at electrode i, and \(\bar{u}\) is the average voltage across all electrodes. EEG topographic maps at the peaks of the GFP curve are selected for clustering, as these maps have a higher signal-to-noise ratio52. The second step involves using a modified K-means algorithm for clustering53, with the number of clusters set from 2 to 8. Each cluster number is run 100 times to minimize inter-run error. The optimal number of clusters is determined using cross-validation criteria54, resulting in several template topographic maps for microstate classes within each group. The third step is to calculate the global map dissimilarity (GMD) between each original EEG signal topographic map within each group and each template topographic map52:

$$\mathrm{GMD}=\sqrt{\frac{\mathop{\sum }\limits_{i=1}^{N}{(\frac{{u}_{i}-\overline{u}}{{\mathrm{GFP}}_{u}}-\frac{{v}_{i}-\overline{v}}{{\mathrm{GFP}}_{v}})}^{2}}{N}}$$
(2)

Where N is the number of electrodes, \({u}_{i}\) is the voltage value of electrodei in topographic map u, \({v}_{i}\) is the voltage value of electrode i in topographic map v, \(\bar{u}\) is the average voltage of all electrodes in topographic map u, and \(\bar{v}\) is the average voltage of all electrodes in topographic map v. The original EEG signal topographic maps are fitted to the templates with the smallest GMD. At this point, the original EEG signals correspond to a sequence of alternating microstate categories. The fourth step involves dividing the EEG signals into non-overlapping windows for each subject based on the microstate sequence, with each window corresponding to a fixed microstate category. The functional connectivity network adjacency matrix is then calculated for each window using the phase locking value (PLV)55:

$$\mathrm{PLV}=\left|\frac{1}{N}\mathop{\sum }\limits_{t=1}^{N}{e}^{i({\theta }_{{xt}}-{\theta }_{{yt}})}\right|$$
(3)

Where N represents the total number of sample points in the EEG signal, \({\theta }_{{xt}}\) denotes the instantaneous phase of signal x(t) at time t, and \({\theta }_{{yt}}\) denotes the instantaneous phase of signal y(t) at time t. The PLV ranges from 0 to 1, with a higher PLV indicating stronger phase synchronization between the two channel signals. The fifth step involves grouping adjacency matrices belonging to the same microstate category into the same set.

Fig. 5: Calculation process of spatiotemporal variability in dynamic brain functional networks based on EEG microstate windows.
figure 5

a Step 1: EEG microstate analysis and division of brain functional network windows. Step 2: Calculate the PLV adjacency matrix for each window and classify it. b Calculate the temporal and spatial variability of dynamic brain functional networks.

Full size image

Temporal and spatial variability of dynamic brain functional networks

Figure 5b illustrates the calculation process for temporal and spatial variability. For each subject, after selecting a particular microstate, its adjacency matrix set is determined as F = {F1, …, Fr…, Fv}. Taking channel i as an example, its temporal variability is defined as:

$${T}_{i}=1-\frac{1}{v\left(v-1\right)}\mathop{\sum }\limits_{p,q=1,p\ne q}^{v}\mathrm{corr}({F}_{p}\left(i,:\right),{F}_{q}\left(i,:\right))$$
(4)

Where \({F}_{p}\left(i,:\right)\) represents the functional architecture of channel i in the pth adjacency matrix, for example, selecting the 15th channel: \({F}_{p}\left(15,:\right)={[{F}_{p}\left(\mathrm{15,1}\right),{F}_{p}\left(\mathrm{15,2}\right),{F}_{p}\left(\mathrm{15,4}\right),\,\ldots ,\,{F}_{p}\left(\mathrm{15,30}\right)]}^{T}\), \(\mathrm{corr}({F}_{p}\left(i,:\right),{F}_{q}\left(i,:\right))\) denotes the correlation coefficient between \({F}_{p}\left(i,:\right)\) and \({F}_{q}\left(i,:\right)\).

Before calculating the spatial variability, it is first necessary to define a spatial functional connectivity sequence:

$${F}_{s}(i,j)={[{F}_{1}\left(i,j\right),{F}_{2}\left(i,j\right),\,\ldots ,\,{F}_{v}\left(i,j\right)]}^{T}$$
(5)

Taking channel i as an example, its spatial variability is defined as:

$${S}_{i}=1-\frac{1}{29\left(29-1\right)}\mathop{\sum }\limits_{j,h=1,j\ne h\ne i}^{v}\mathrm{corr}({F}_{s}\left(i,j\right),{F}_{s}\left(i,h\right))$$
(6)

Where \(\mathrm{corr}({F}_{s}\left(i,j\right),{F}_{s}\left(i,h\right))\) represents the correlation coefficient between \({F}_{s}\left(i,j\right)\) and \({F}_{s}\left(i,h\right)\).

Temporal variability reflects the temporal characteristics of functional structure changes at specific brain regions, while spatial variability reflects the spatial characteristics of changes in spatial functional connectivity sequences related to specific brain regions.

Statistical analysis

Statistical analyses were conducted using SPSS 26.056. The Shapiro-Wilk test was employed for normality testing, with parametric tests (independent samples t-test/paired samples t-test) applied to normally distributed data and non-parametric tests (Mann–Whitney U test/Wilcoxon signed-rank test) used for non-normally distributed data. Spearman’s correlation analysis was performed for correlation assessment. All p-values underwent Benjamini–Hochberg false discovery rate (FDR) correction, with the significance level set at p < 0.05 after correction. For between-group comparisons, independent sample t-tests and Mann–Whitney U tests were supplemented by Cohen’s d and Cliff’s delta effect size calculations, respectively. For within-group comparisons, Hedges ‘g and Rank-biserial effect sizes were reported with the paired-samples t-test and the Wilcoxon signed-rank test. Using G*Power software (version 3.1.9.7) for sensitivity analysis57, we calculated the effect size threshold with parameters set at α = 0.05, statistical power (1 − β) = 0.8, and sample sizes of N1 = 51 (PD) and N2 = 21 (HC), yielding an effect size threshold of r = 0.353. Therefore, in analyzing the differences in spatiotemporal variability between PD patients before treatment and healthy controls (PD vs HC between-group comparison), we only discussed results showing statistical significance (p < 0.05) with effect sizes r > 0.353. Additionally, we performed a priori sample size calculation for the 10Hz-rTMS and iTBS groups, with parameters set at α = 0.05, statistical power (1 − β) = 0.8, and effect size r = 0.3, which indicated a minimum requirement of 23 patients per group. Both the 10Hz-rTMS group (n = 25) and iTBS group (n = 26) met this requirement. Consequently, when comparing pre- and post-treatment differences in spatiotemporal variability within the 10Hz-rTMS and iTBS groups (within-group comparison), we only analyzed and discussed results demonstrating statistical significance (p < 0.05) with effect sizes r > 0.3.