The Temporal Dynamics of EEG Microstate Reveals the Neuromodulation Effect of Acupuncture With Deqi

The electroencephalography (EEG) microstate has recently emerged as a new whole-brain mapping tool for studying the temporal dynamics of the human brain. Meanwhile, the neuromodulation effect of external stimulation on the human brain is of increasing interest to neuroscientists. Acupuncture, which originated in ancient China, is recognized as an external neuromodulation method with therapeutic effects. Effective acupuncture could elicit the deqi effect, which is a combination of multiple sensations. However, whether the EEG microstate could be used to reveal the neuromodulation effect of acupuncture with deqi remains largely unclear. In this study, multichannel EEG data were recorded from 16 healthy subjects during acupuncture manipulation, as well as during pre- and post-manipulation tactile controls and pre- and post-acupuncture rest controls. As the basic acupuncture unit for regulating the central nervous system, the Hegu acupoint was used in this study, and each subject’s acupuncture deqi behavior scores were collected. To reveal the neuroimaging evidence of acupuncture with deqi, EEG microstate analysis was conducted to obtain the microstate maps and microstate parameters for different conditions. Furthermore, Pearson’s correlation was analyzed to investigate the correlation relationship between microstate parameters and deqi behavioral scores. Results showed that: (1) compared with tactile controls, acupuncture manipulation caused significantly increased deqi behavioral scores. (2) Acupuncture manipulation significantly increased the duration, occurrence, and contribution parameters of microstate C, whereas it decreased those parameters of microstate D. (3) Microstate C’s duration parameter showed a significantly positive correlation with acupuncture deqi behavior scores. (4) Acupuncture manipulation significantly increased the transition probabilities with microstate C as node, whereas it reduced the transition probabilities with microstate D as node. (5) Microstate B→C’s transition probability also showed a significantly positive correlation with acupuncture deqi behavior scores. Taken together, the temporal dynamic feature of EEG microstate could be used as objective neuroimaging evidence to reveal the neuromodulation effect of acupuncture with deqi.

The electroencephalography (EEG) microstate has recently emerged as a new whole-brain mapping tool for studying the temporal dynamics of the human brain. Meanwhile, the neuromodulation effect of external stimulation on the human brain is of increasing interest to neuroscientists. Acupuncture, which originated in ancient China, is recognized as an external neuromodulation method with therapeutic effects. Effective acupuncture could elicit the deqi effect, which is a combination of multiple sensations. However, whether the EEG microstate could be used to reveal the neuromodulation effect of acupuncture with deqi remains largely unclear. In this study, multichannel EEG data were recorded from 16 healthy subjects during acupuncture manipulation, as well as during pre-and post-manipulation tactile controls and pre-and postacupuncture rest controls. As the basic acupuncture unit for regulating the central nervous system, the Hegu acupoint was used in this study, and each subject's acupuncture deqi behavior scores were collected. To reveal the neuroimaging evidence of acupuncture with deqi, EEG microstate analysis was conducted to obtain the microstate maps and microstate parameters for different conditions. Furthermore, Pearson's correlation was analyzed to investigate the correlation relationship between microstate parameters and deqi behavioral scores. Results showed that: (1) compared with tactile controls, acupuncture manipulation caused significantly increased deqi behavioral scores. (2) Acupuncture manipulation significantly increased the duration, occurrence, and contribution parameters of microstate C, whereas it decreased those parameters of microstate D. (3) Microstate C's duration parameter showed a significantly positive correlation with acupuncture deqi behavior scores. (4) Acupuncture manipulation significantly increased the transition probabilities with microstate C as node, whereas it reduced the transition probabilities with microstate D as node. (5) Microstate B→C's transition probability also showed a significantly positive correlation with acupuncture deqi behavior scores. Taken together, the temporal dynamic feature of EEG microstate could be used as objective neuroimaging evidence to reveal the neuromodulation effect of acupuncture with deqi.

INTRODUCTION
The dynamics and functions of the nervous system could be efficiently modulated by external stimuli, which has received increasing attention from neuroscientists (Marder, 2012). Acupuncture, as an external stimulation to the human nerve system, has been applied to manage various clinical diseases, such as neuropathic pains (Miranda et al., 2015), epilepsy (Nicolaou and Georgiou, 2012), Alzheimer's disease (Liang et al., 2014), and Parkinson's disease (Aroxa et al., 2017), which has received increasing attention by the public. Moreover, previous studies showed that acupuncture could modulate the human brain by balancing the parasympathetic and sympathetic systems (Wu et al., 1999;Ernst, 2006;Takahashi, 2011). Although the World Health Organization announced that more than 43 diseases could be treated with acupuncture (Hsu et al., 2011), there is a lack of objective neuroimaging evidence for revealing the neuromodulation effect of acupuncture.
According to traditional Chinese medicine (TCM), a fine needle is inserted into the subcutaneous muscle at specific acupoints during acupuncture. Acupuncture could stimulate the afferent nerve and regulate the central nervous system, in which subjects experience the sensation of acupuncture with deqi (Takahashi, 2011). Moreover, acupuncture to achieve the deqi effect is called effective acupuncture, in which the deqi effect is created by a combination of different unique senses (Hui et al., 2007). And, the deqi effect is positively correlated with the therapeutic effect (Takahashi, 2011). Studies have shown that the deqi effect is an essential factor in acupuncture's regulation and treatment of various diseases (Fang et al., 2009;Asghar et al., 2010). However, the current method for evaluating the acupuncture effect mainly relies on the subjective deqi behavioral questionnaire (Hui et al., 2007;Yu et al., 2012). Therefore, there lacks a reliable objective neuroimaging evidence to reveal the neuromodulation effects of acupuncture with deqi.
To reveal the neural mechanism of acupuncture, functional magnetic resonance imaging (fMRI) and electroencephalograph (EEG) were employed in previous studies (Paraskeva et al., 2004;Onton et al., 2005;Tyvaert et al., 2009). The high-spatialresolution fMRI studies helped reveal the specific brain areas involved in acupuncture. It was found that acupuncture could activate the insula (Wu et al., 1999;Fang et al., 2009;Chae et al., 2013) and anterior cingulate cortex (Wu et al., 1999;Fang et al., 2009;Chae et al., 2013), whereas it could deactivate the medial prefrontal cortex, caudate, amygdala, and posterior cingulate cortex (Wu et al., 1999;Fang et al., 2009;Chae et al., 2013). These brain response patterns indicated that acupuncture could modulate cognitive functions such as pain perception, memory, and emotional dimensions (Fang et al., 2009;Chae et al., 2013). However, due to fMRI's poor temporal resolution, fMRI studies failed to reveal the temporal dynamic changes during acupuncture. In addition, studies have found that acupuncture can strengthen brain connectives (Yu et al., 2017(Yu et al., , 2018. However, these studies only answer the connectivity between static network nodes, and the study of dynamic characteristics is still unclear. Previous EEG studies mainly focused on the specific frequency band's neural oscillation changes, in which EEG studies found increased power changes of delta and alpha frequency bands during acupuncture by using frequency analysis methods (Chen et al., 2006;Bian et al., 2011;Yu et al., 2018). Although the high temporal resolution EEG could capture the brain dynamics, previous EEG studies neglected the temporal characteristics during acupuncture. In summary, previous fMRI and EEG neuroimaging studies showed that acupuncture could modulate the neural activities of specific brain areas and regulate the EEG frequency bands. However, previous studies have not revealed the temporal dynamics of the human brain during acupuncture. Therefore, how the human brain is dynamically modulated by acupuncture manipulation remains largely unknown.
The human brain is essentially a dynamic network system (Bassett and Sporns, 2017). The communication within and between large-scale brain networks will fluctuate spontaneously at all times in a well-organized way (Bassett and Sporns, 2017;Kucyi et al., 2017;Hsu et al., 2018). Acupuncture is a way to regulate the neural activity of the brain. Currently, more studies focus on the response of different network nodes and brain regions using fMRI (Wu et al., 1999;Fang et al., 2009;Chae et al., 2013). However, the temporal dynamics of the brain during acupuncture with deqi remain unclear. Brain dynamics can be explored by using different neuroimaging methods with different timescales. While the high-spatial-resolution fMRI captures the slow neural dynamics on a second timescale, the faster millisecond timescale neural dynamic changes could be explored by using the EEG microstate analysis (Michel and Koenig, 2018;Schumacher et al., 2019). Previous studies showed that the dynamic patterns of EEG activities are empirically divided into a series of alternating types of quasi-stable EEG topographic maps, which are termed "EEG microstates" (Lehmann et al., 1987;Khanna et al., 2014;Michel and Koenig, 2018). It was proved that there existed four discrete EEG microstates, i.e., microstates A, B, C, and D, which were found to be reproducible across participants (Koenig et al., 2002). The four microstates (A, B, C, and D) showed maximum activity from right frontal to left posterior, left frontal to right posterior, midline frontal occipital, and central frontal (Lehmann et al., 2009). And they are related to four brain networks [auditory network, visual network, salience network, and dorsal attention network (DAN)] . Each EEG microstate remains stable for a short period of 80-120 ms before rapidly transitioning into another microstate (Lehmann et al., 1987;Khanna et al., 2014;Michel and Koenig, 2018), which could detect the discontinuous and nonlinear dynamic changes of brain activities (Lehmann et al., 1987). Particularly, the EEG microstate has been considered as the basic units of human cognition and was suggested as the "atom of thought, " which gained increasing popularity in recent years (Koenig et al., 1999(Koenig et al., , 2002Michel and Koenig, 2018).
Importantly, the parameters of the alternating EEG microstate time series, such as duration, occurrence, and transition probability, were proved to be effective markers both for distinguishing healthy subjects' different cognition states (Milz et al., 2016) and for detecting the specific brain changes associated with neuropsychiatric diseases (Koenig et al., 1999;Lehmann et al., 2005;Nishida et al., 2013). For example, different mental states, such as different meditation states  or sleep conditions (Brodbeck et al., 2012;Xu et al., 2020), were found to be correlated with different EEG microstate parameters for healthy subjects. Furthermore, several neuropsychiatric diseases, such as Alzheimer's disease (Nishida et al., 2013), Parkinson's disease (Chu et al., 2020), dementia (Schumacher et al., 2019), and schizophrenia (Tomescu et al., 2014;Murphy et al., 2020), could also be characterized by using EEG microstate parameters. Microstates as a way to describe brain dynamics may be affected by acupuncture, thus providing an objective measure of its effect on brain dynamics (Hsu et al., 2018;Abreu et al., 2021). However, to our knowledge, there is no previous research employing the EEG microstate theory to investigate the brain dynamics during acupuncture with deqi. Therefore, whether the EEG microstate parameters could be used as markers to reveal the dynamic neuromodulation effects of acupuncture with deqi remains largely unclear.
To reveal the temporal dynamics of brain activity during acupuncture with deqi, multichannel EEG data were recorded from 16 healthy subjects during acupuncture manipulation, as well as during pre-and post-manipulation tactile controls and pre-and post-acupuncture rest controls. In addition, the Hegu point is one of the most commonly used acupuncture points, located at the midpoint on the radial side of the second metacarpal, which is also known as large intestine 4 (LI4) (Kong et al., 2015). Considering that the Hegu acupoint is the basic neural acupuncture unit that modulates the central nervous system (Zhang et al., 2012), the Hegu point was explored in this study. The subject's acupuncture deqi behavior scores were collected. EEG microstate analysis was conducted to obtain the microstate maps and microstate parameters for different conditions. Furthermore, Pearson's correlation was analyzed to investigate the correlation relationship between microstate parameters and deqi behavior scores. We speculate that the EEG microstate theory could provide a new perspective to help understand the temporal dynamics of the human brain during acupuncture with deqi.

Subjects
Sixteen healthy right-handed participants (24.0 ± 1.9 years old [mean ± SD], age range from 22 to 28 years, 10 males and 6 females) were recruited for this study. All the participants had normal intelligence, and none of the following conditions existed: (a) history of mental illness, (b) brain damage caused by longterm use of drugs, and (c) experience of acupuncture treatment. All participants signed the informed consent approved by the Institutional Review Board and Ethics Committee of Tianjin University. All participants agreed to provide data anonymously. Supplementary Table S1 shows the details of the subjects.

Acupuncture Experiment Paradigm
Each subject was required to sit in a comfortable chair during the entire process of the acupuncture experiment in a bright and quiet room. Throughout the acupuncture experiment, subjects were asked to keep awake and relaxed during the preparation stage. All acupunctures were performed by an experienced acupuncturist. Once the subject is found to be uncomfortable or not awake, the subject will be timely reminded and asked about the situation.
During the preparation stage, each subject was asked to be familiar with the visual analogue scale (VAS) behavior questionnaire for evaluating the acupuncture deqi sensations, which includes the following six sensations: soreness, numbness, distention, heaviness, spread, and dull pain according to the TCM (Hui et al., 2007;Yu et al., 2012). The intensity of each sensation is rated by a point, ranging from 0 to 10 (0-3 mild, 4-6 moderate, 7-8 strong, 9 severe, and 10 intolerable), on a continuous horizontal line on the VAS questionnaire (Lukacz et al., 2004;Hui et al., 2007). Considering that the unique feeling of deqi appears through the combination of different detailed sensations (soreness, numbness, distention, heaviness, spread, and dull pain), the comprehensive index could comprehensively characterize the deqi effect (Hui et al., 2007).
Considering that the quasi-experimental design could compare the condition difference before and after time (Baernighausen et al., 2017), the quasi-experimental design was used in this experiment (Figure 1), which was in line with previous acupuncture researches (Yu et al., 2017(Yu et al., , 2018. The whole acupuncture experiment included the following five stages: (I) pre-acupuncture rest control: each subject was asked to maintain awake and rest for 5 min without movement before the needle was inserted by the acupuncturist. (II) Premanipulation tactile control: a single-use sterile acupuncture needle was inserted by the acupuncturist at the Hegu acupoint to a certain depth without twist randomly on the left or right hand for each subject. During the tactile control stage, each subject was asked to maintain awake and rest for 5 min before the acupuncture twist starts by the acupuncturist. (III) Acupuncture manipulation: the lifting thrusting manipulation was performed for 2 min at the Hegu acupoint with the acupuncture needle twist frequency of 2 Hz. (IV) Post-manipulation tactile control: after the acupuncture manipulation period, each subject was asked to maintain awake and rest for 5 min. During this post-twist tactile control stage, the needle was kept at the Hegu acupoint without the twist. (V) Post-acupuncture rest control: after the acupuncturist removed the needle, each subject was asked to maintain awake and rest for 5 min without movement.
Immediately after the acupuncture experiment, the VAS deqi sensation behavior questionnaire was given to each subject to evaluate his/her deqi sensation for the following three stages: pre-manipulation tactile control, acupuncture manipulation, and post-manipulation tactile control (Supplementary Table S2).

The Deqi Behavioral Index Calculation and Statistical Analysis
To comprehensively characterize the needing sensations, the deqi index was used for quantifying the deqi effect. The deqi behavioral index was defined as the summation of all the six individual FIGURE 1 | Acupuncture experiment paradigm. During the preparation stage, each subject was asked to be familiar with the behavior questionnaire. There were five conditions in the experiment: 5-min pre-acupuncture rest control, 5-min pre-manipulation tactile control, 2-min acupuncture manipulation, 5-min post-manipulation tactile control, and 5-min post-acupuncture rest control. After the experiment, subjects were asked to evaluate the needling sensations for the pre-manipulation tactile control, acupuncture manipulation, and post-manipulation tactile control. The twisting frequency in the acupuncture manipulation was 2 Hz.
needing sensation value according to the following equation (Supplementary Table S2): where S i indicates the rating points for the ith sensation and N represents the total number of different deqi sensations.
To evaluate if there exists significant deqi sensation difference between different conditions, the two-sample t-test was used on all subjects' deqi indexes in MATLAB version R2016a for the following three stages: pre-manipulation tactile control, acupuncture manipulation, and post-manipulation tactile control. The p-Values were corrected with Bonferroni correction for multiple comparisons.

Electroencephalography Recordings
The EEG signals were recorded by a Neuroscan SynAmps2 amplifier with 64 channels in the international 10-20 system. The reference electrode (REF) was placed on the right earlobe, and the ground electrode (GND) was placed on the forehead. The sampling frequency of the EEG signal was 1,000 Hz. For removing power line noise during the EEG data acquisition, the 50-Hz hardware notch filter was used.

Electroencephalography Data Preprocessing
All the EEG data were preprocessed by using the EEGLAB toolbox (biosig extension) in MATLAB R2016a environment (The MathWorks, Natick, MA, United States). The data were first 0.5-60 Hz bandpass filtered with a 6,600th-order finite impulse response (FIR) filter and then were filtered by a 50-Hz notch filter to remove the power line noise (Kikuchi et al., 2011;Pedroni et al., 2017). Then, the blind source separation (BSS) algorithm was used to remove the interference of eye movements, electromyography (EMG), and other artifacts from the EEG data affected by artifacts (Gao et al., 2018). The BSS is based on second-order blind identification (SOBI)/canonical correlation analysis (CCA) to automatically remove electrooculogram (EOG)/EMG as described by previous studies (De Clercq et al., 2005;Gomez-Herrero et al., 2006). The EEG data were then rereferenced to the average reference and were bandpass filtered between 2 and 20 Hz with a 1,650th-order FIR filter (Koenig et al., 2002). After filtering, the data were segmented into 2-s epochs, in which epochs with an amplitude of more than ±80 µV at any electrode were rejected (Koenig et al., 2002;Gao et al., 2018).

Electroencephalography Microstate Analysis
The EEG microstate analysis was processed with the EEGLAB's microstate toolbox 1 in the MATLAB R2016a environment. Microstate analysis followed the procedures described in Koenig et al. (1999); Brodbeck et al. (2012).

Global Field Power Calculation
As the first step, the global field power (GFP) at each time point was calculated on the preprocessed EEG data, which describes the spatial standard deviation of the EEG signal across all electrodes. The time series of GFP reflects the total energy change in the EEG topographic map over time. The GFP peaks have the highest energy with the optimal topographic signalto-noise ratio (Lehmann and Skrandies, 1980;Koenig et al., 2002), which can characterize momentary quasi-stable voltage topography (Skrandies, 1989;Koenig and Brandeis, 2016).
(N = number of electrodes, V i (t) = measured voltage of electrode i at time t, and i = electrode i).

The k-Means Clustering of Electroencephalography Topographies
To get individual cluster maps, the k-means clustering analysis was first calculated on each individual subject for each condition separately. The GFP peaks were used as the initial topographic maps for clustering to maximize the topographic signal-to-noise ratio (Figure 2A). For each subject, the topographic map at GFP peaks was sent to the unsupervised k-means clustering algorithm (Koenig et al., 2002;Murray et al., 2008). The number of microstate classes was set to 4 in line with previous classical EEG microstate studies (Koenig et al., 2002;Brodbeck et al., 2012;Mishra et al., 2020), which could explain the most variance of map topographies during EEG microstate clustering (Koenig et al., 2002).
To identify the specific group maps for each condition, the individual cluster maps of all subjects were then averaged FIGURE 2 | Schematic illustration of electroencephalography (EEG) microstate analysis. (A) The global field power (GFP) of spontaneous EEG was calculated, and all topographic maps at the maximum of local GFP were clustered by k-means algorithm. (B) Group template microstate maps according to the k-means clustering algorithm. (C) By computing the spatial correlation between the topographical at each GFP peak and each of the group microstate maps, the GFP peak was assigned to the microstate class with the highest correlation value. Then the microstate class sequence was obtained, which could be used to calculate the microstate parameters of each subject for each condition.
using a permutation algorithm. The algorithm determines the solution with the largest average correlation across all four classes of subjects, which is consistent with the previous study (Koenig et al., 1999;Milz et al., 2016). Finally, a specific set of group template microstate maps for each condition is obtained separately ( Figure 2B).

Spatial Correlation of Electroencephalography Map Topographies and Sequence
To obtain the class label of each topography map at the GFP peak of original data and the corresponding class label sequence, the group template microstate maps were then fit back to each subject's topographic maps at the GFP peaks by computing the spatial correlation between topographies ( Figure 2C). To calculate the spatial correlation between the EEG microstate maps in our study with the corresponding classical microstate map templates (Koenig et al., 2002), the spatial correlation analysis method was applied in this study, which was consistent with the previous study (Brodbeck et al., 2012). The spatial correlation used the spatial Pearson's product moment correlation coefficient (Brandeis et al., 1992), which was computed as follows: (C = spatial correlation, N = number of electrodes, u = measured voltage map u, v = measured voltage map v, and i = electrode i).
The GFP peak was assigned to the microstate class with the highest spatial correlation value, and the middle time point between the two GFP peaks was set as the start and end points of each microstate (Koenig et al., 1999;Murray et al., 2008). In addition, the microstates at the beginning and end of each epoch were likely to be incomplete and were excluded from further analysis (Koenig et al., 2002).
Finally, the microstate time series of all subjects in each condition could be obtained through the microstate class label of each GFP peak in the original data ( Figure 2C).

Electroencephalography Microstate Parameters
According to the calculated microstate time series, the following four microstate parameters can be calculated for each condition separately.
-Duration: the average time length during which each microstate remains stable. The microstate duration parameter is considered to reflect the average time for a set of neural generators to maintain synchronized activities (Lehmann et al., 1987;Strelets et al., 2003). -Occurrence: during a given period, the average number of occurrences per second for each microstate. The microstate occurrence parameter represents the tendency of a set of neural generators to coordinate their activities over time (Lehmann et al., 1987;Khanna et al., 2014;Michel and Koenig, 2018). -Contribution (time coverage): the proportion of each microstate's occurrence time over the total time. The microstate contribution parameter is interpreted as reflecting the time coverage of its underlying neural generator relative to other neural generators (Lehmann et al., 1987). In addition, for the three parameters duration, occurrence, and contribution, knowing any two of them could calculate the remaining parameter (Lehmann et al., 2005). -Transition probability ("microstate syntax"): percentage of transitions from one microstate class to another over all transitions occurring during a given period. The microstate transition probability is usually interpreted as different activation sequences of neural components that generate microstates (Lehmann et al., 2005;Khanna et al., 2015). If the transition from the previous microstate class to the next microstate class occurred randomly, then the relative occurrence of the microstate class would be proportional to the transition value (Lehmann et al., 2005).

Statistical Analysis
Statistical Analysis of Microstate Parameters (Duration, Occurrence, and Contribution) For each of the three microstate parameters (duration, occurrence, and contribution), separate 5 × 4 two-way repeatedmeasures analysis of variance (rmANOVA) was conducted using condition (pre-acupuncture rest control, pre-manipulation tactile control, acupuncture manipulation, post-manipulation tactile control, or post-acupuncture rest control) and microstate class (A, B, C, or D) as repeated measures (Table 1) in SPSS (Statistical Package for the Social Sciences v26, IBM). Mauchly's test of sphericity was conducted during rmANOVA, and the Greenhouse-Geisser correction was applied when the sphericity assumption was not valid. Then, the post hoc paired t-tests were employed when the main effects or interactions in the rmANOVA were significant. To minimize the risk of type I errors, the p-Values were corrected with Bonferroni correction for multiple comparisons (Tomescu et al., 2014;Faber et al., 2017;Fu et al., 2018).

Statistical Analysis of Microstate Parameter (Transition Probability)
For each of the 12 transition probability parameters (A→B, A→C, A→D, B→A, B→C, B→D, C→A, C→B, C→D, D→A, D→B, and D→C), separate one-way five-level rmANOVA was conducted using condition (pre-acupuncture rest control, premanipulation tactile control, acupuncture manipulation, postmanipulation tactile control, or post-acupuncture rest control) as repeated measures ( Table 2). Mauchly's test of sphericity was conducted during rmANOVA, and the Greenhouse-Geisser correction was applied when the sphericity assumption was not valid. The above statistical analysis was computed using SPSS (Statistical Package for the Social Sciences v26, IBM).
To further compare the transition probability parameters for different conditions, the two-sample t-tests were performed on the transition probabilities between acupuncture manipulation and all controls, and those between tactile controls and rest controls in MATLAB version R2016a, separately. To minimize the risk of type I errors, the p-Values were corrected with Bonferroni correction for multiple comparisons.

Correlation Analysis Between Neural Responses and Acupuncture Behavior Scores
To further investigate if the neural responses are significantly correlated with acupuncture's behavior scores, Pearson's correlation and linear regression analysis were conducted between deqi behavior data and EEG microstate parameters (duration, occurrence, contribution, and transition probability) for acupuncture manipulation, and pre-and post-manipulation tactile controls, separately. To keep the sample size consistent during Pearson's correlation calculating under different conditions, the x-axis (behavior) and y-axis (microstate parameters) values of tactile controls' each data sample were obtained by averaging the corresponding x-axis and y-axis data of the pre-and post-manipulation tactile controls, respectively. During the linear regression analysis, the 95% confidence intervals for the mean of the polynomial evaluation were also computed.

Behavior Deqi Results for Different Conditions
To validate the behavior effectiveness of acupuncture manipulation, we compared the deqi behavior scores of  the acupuncture manipulation with those of the preand post-manipulation tactile controls (Figure 3 and Supplementary Table S2). The subjects' deqi indexes of the acupuncture manipulation were significantly larger than those of pre-manipulation tactile control (p = 0.036, paired t-test, Bonferroni correction) and were significantly larger than those of post-manipulation tactile control (p = 0.037, paired t-test, Bonferroni correction). There was no significant difference between pre-manipulation tactile control and post-manipulation tactile control. These results indicated that acupuncture manipulation condition could significantly cause behavioral effects compared with pre-and post-manipulation tactile control conditions.

The Electroencephalography Microstate Results Were Consistent With Previous Classical Findings
To investigate whether this study's EEG microstate maps are consistent with those of previous results, microstate map comparisons and the spatial correlations were analyzed. The microstate analysis was performed on the entire data of all five conditions, and the corresponding microstate maps were obtained (Figure 4), in which the four microstate maps were marked as A, B, C, and D according to the maximum spatial correlation value with the previous classical microstate maps (Koenig et al., 2002). Microstate map A showed a left occipital to right frontal orientation, and microstate B showed a right occipital to left frontal orientation. Microstate C showed an approximately symmetric occipital to frontal orientation, and microstate D was characterized by a frontal-central maximum. The characteristics of this study's four microstate maps were in line with previous findings (Koenig et al., 2002;Britz et al., 2010). Furthermore, the spatial correlations on rest condition's microstate maps between this study and previous publication ( Figure 4A) were analyzed to investigate whether our microstate results were consistent with previous classical results (Koenig et al., 2002). It was showed that our four microstate maps of the pre-acupuncture rest control had relatively high spatial correlation values (mean ± SD, 93% ± 3.42%; A = 96%, B = 88%, C = 96%, and D = 91%) with maps described in Koenig et al. (2002) (Figure 4A).

Correlation Results Between Microstate Maps of Different Conditions
To investigate the microstate map changes during different conditions, spatial correlations of map topographies between the pre-acupuncture rest control and the next four conditions were analyzed for each microstate class separately. It was shown that compared with the pre-acupuncture rest control's microstate maps, the other four conditions' microstate maps overall showed relatively high spatial correlation (the premanipulation tactile control, mean ± SD, 98% ± 1.12%; the post-manipulation tactile control, mean ± SD, 97% ± 2.29%; acupuncture manipulation, mean ± SD, 95% ± 1.22%; and postacupuncture rest control, mean ± SD, 99% ± 0.43%) ( Figure 4B). Interestingly, the acupuncture manipulation's four class maps showed the lowest concordance with the pre-acupuncture rest control maps (mean ± SD, 95% ± 1.22%), in which the spatial correlation percent for microstates C and D was the smallest among all results (C = 94% and D = 94%) ( Figure 4B). On the other hand, the post-acupuncture rest control maps showed the overall highest concordance with the pre-acupuncture rest control's microstate maps (mean ± SD, 99% ± 0.43%). These phenomena showed that microstate maps were first modulated by acupuncture manipulation and then recovered back to the baseline resting state. All these results indicated that the EEG microstate maps could be a promising objective marker to reveal the acupuncture effects.

Acupuncture Manipulation Increased the Duration, Occurrence, and Contribution Parameters of Microstate C, Whereas It Decreased Those of Microstate D
To further investigate the spatial correlation changes for microstates C and D, the time series of microstate C's and D's parameters (duration, occurrence, and contribution) changes over time in different conditions were obtained ( Figure 5). Results showed different parameters' change patterns for microstates C and D under acupuncture manipulation. Microstate C's duration, occurrence, and contribution parameters had tremendous increases during acupuncture manipulation than the other four controls (Figures 5A,C,E).
On the contrary, microstate D's duration, occurrence, and contribution parameters had obvious decreases during acupuncture manipulation as compared with other controls (Figures 5B,D,F). These results indicated that the microstate dynamics parameters could reveal the acupuncture's neural modulation effects on the human brain.

Significant Condition × Microstate Class Interaction Effect for Duration, Occurrence, and Contribution Parameters
To further quantify the microstate parameter changes between different conditions (pre-acupuncture rest control, premanipulation tactile control, acupuncture manipulation, post-manipulation tactile control, and post-acupuncture FIGURE 4 | The spatial correlation percent and electroencephalography (EEG) microstate maps during different conditions. (A) The percent spatial correlation of this study's pre-acupuncture rest control microstate topography maps with the rest of the maps of previous publications (Koenig et al., 2002). The percent spatial correlation is given below each microstate class map. (B) The EEG microstate maps during different conditions. The percent spatial correlation between pre-acupuncture rest control maps with the other four conditions maps is given below each map. The average correlation for each condition is given at the top of each column with bold type. Colors represent the relative potential distribution.

The Duration, Occurrence, and Contribution Parameters of Microstate C Were Significantly Increased by Acupuncture Manipulation
The post hoc rmANOVA results showed that microstate C's parameters (duration, occurrence, and contribution) of acupuncture manipulation were significantly larger than those of other conditions (Figure 6 and Table 1). First, microstate C's duration of acupuncture manipulation was significantly larger than that of the other conditions (pre-acupuncture rest control, p = 0.008 * ; pre-manipulation tactile control, p = 0.010 * ; post-manipulation tactile control, p = 0.011 * ; post-acupuncture rest control, p = 0.004 * * ; alpha = 0.05; Bonferroni correction) ( Figure 6A and Table 1). Second, microstate C's occurrence of acupuncture manipulation was also significantly higher than that of the other conditions (pre-acupuncture rest control, p = 0.032 * ; pre-manipulation tactile control, p = 0.006 * ; postmanipulation tactile control, p = 0.004 * * ; post-acupuncture rest control, p = 0.008 * ; with Bonferroni correction) ( Figure 6B and Table 1). Finally, microstate C's contribution of acupuncture manipulation was significantly higher than that of the other conditions (pre-acupuncture rest control, p = 0.003 * * ; premanipulation tactile control, p = 0.002 * * ; post-manipulation tactile control, p = 0.003 * * ; post-acupuncture rest control, p = 0.001 * * ; with Bonferroni correction) ( Figure 6C and Table 1). These changes were consistent with the previous time series results of microstate C's parameters (duration, occurrence, and contribution) increasing changes during acupuncture manipulation (Figure 5).

The Duration, Occurrence, and Contribution Parameters of Microstate D Were Significantly Decreased by Acupuncture Manipulation
On the contrary, the post hoc rmANOVA results showed that microstate D's parameters (duration, occurrence, and contribution) of acupuncture manipulation were significantly smaller than those of other conditions (Figure 6 and Table 1). First, microstate D's duration of acupuncture manipulation was smaller than that of other conditions and was significantly lower than that of the post-acupuncture rest control (p = 0.016 * ; with Bonferroni correction) ( Figure 6A and Table 1). Second, microstate D's occurrence of acupuncture manipulation was lower than that of other conditions and was significantly lower than that of the other conditions (pre-acupuncture rest control, p = 0.030 * ; post-manipulation tactile control, p = 0.049 * ; post-acupuncture rest control, p = 0.017 * ; with Bonferroni correction) ( Figure 6B and Table 1). Finally, microstate D's contribution of acupuncture manipulation was significantly lower than that of the other conditions (pre-acupuncture rest control, p = 0.012 * ; pre-manipulation tactile control, p = 0.030 * ; post-manipulation tactile control, p = 0.026 * ; postacupuncture rest control, p = 0.005 * ; with Bonferroni correction) ( Figure 6C and Table 1). These changes were consistent with the previous time series results of microstate D's parameters (duration, occurrence, and contribution) decreasing changes during acupuncture manipulation (Figure 5).

The Duration Parameter of Microstate C Was Significantly Positively Correlated With the Behavior Score During Acupuncture Manipulation
To further investigate if there existed relationships between neural responses (parameters of EEG microstates C and D) and acupuncture behavior scores, Pearson's correlations were performed between the microstate parameters (duration, occurrence, and contribution) and each subject's behavior deqi index. Interestingly, the duration parameter of microstate C was found to have a significant positive correlation with acupuncture's behavior deqi index during the acupuncture manipulation (r = 0.521, p = 0.038), whereas there was no significant correlation for the tactile controls (r = 0.156, n.s.) (Figure 7). Besides, the contribution parameter of microstate C had a positive correlation tendency (r = 0.433, n.s.) with behavior deqi index during acupuncture manipulation compared with tactile controls (r = −0.136, n.s.) (Supplementary Figure S1). Besides, no significant correlation was found for any other conditions. The neural behavior correlation results indicated that the duration parameter of microstate C could be an effective marker to reflect subjects' acupuncture deqi behavior scores during acupuncture manipulation. FIGURE 7 | Correlation analysis between neural responses and acupuncture's behavior scores across all subjects. Pearson's correlations between the duration parameter of microstate C and the deqi behavior index for the acupuncture manipulation and tactile controls. Each red dot is indicated for each subject's data. Dashed lines represent the 95% confidence intervals for the mean of polynomial evaluation. Solid lines represent linear regression (p < 0.05, n = 16).
Interestingly, the acupuncture manipulation showed a significant increase for the transition probabilities with microstate C as the node [A→C, t(78) = 4.49, p = 3.0 × 10 −4 ; B→C, t(78) = 3.89, p = 0.003; C→A, t(78) = 4.40, p = 4.2 × 10 −4 ; and C→B, t(78) = 4.26, p = 6.7 × 10 −4 ; two-sample t-tests with Bonferroni correction], compared with all controls (combining four controls) ( Figure 8B and Table 2). On the other hand, A B FIGURE 8 | Transition probability trend and statistical comparisons for different conditions. (A) Trends of the transition probability for the acupuncture manipulation and the remaining four controls, in which the lines represent the switching probability from a microstate class to another. For each of the 12 transition probability parameters, separate one-way five-level repeated-measures analysis of variance (rmANOVA) was conducted using condition as repeated measures (mean ± SEM, *p < 0.05, n = 16). (B) Schematic arrangement and statistical comparisons of the transition probability. Left: comparison between acupuncture manipulation (n = 16) and all controls (combining four controls, n = 64). Right: comparison between tactile controls (combining two tactile controls, n = 32) and rest controls (combining two rest controls, n = 32). Two-sample t-tests with Bonferroni correction for multiple comparisons, *p < 0.05; the number of multiple comparisons = 12.
the acupuncture manipulation showed a significant decrease for the transition probabilities with microstate D as the node [A→D, t(78) = −3.86, p = 0.003; B→D, t(78) = −3.85, p = 0.003; D→A, t(78) = −3.96, p = 0.002; and D→B, t(78) = −4.26, p = 6.8 × 10 −4 ; two-sample t-tests with Bonferroni correction], compared with all controls (combining four controls) ( Figure 8B and Table 2). For the control group, there was no significant difference for the transition probability between tactile controls and rest controls ( Figure 8B). In addition, compared with rest controls (combining two rest controls) or with tactile controls (combining two tactile controls), the acupuncture manipulation also showed significant increase for the transition probabilities with microstate C as the node, but significant decrease for the transition probabilities with microstate D as the node (Supplementary Figure S2).

The Parameter of Transition Probability B→C Was Significantly Positively Correlated With the Behavior Score During Acupuncture Manipulation
To further investigate if there existed relationships between neural responses (parameters of transition probability with microstates C and D as nodes) and acupuncture behavior scores, Pearson's correlations were performed between the microstate's transition probability parameters and each subject's behavior deqi index. Interestingly, the transition probability from microstate B to C was found to have a significant positive correlation with acupuncture's behavior deqi index during the acupuncture manipulation (r = 0.536, p = 0.032), whereas there was no significant correlation for the tactile controls (r = 0.028, n.s.) (Figure 9). Besides, no significant correlation was found for any other conditions. These results further indicated that B→C transition probability parameters could be an effective marker to reflect subjects' acupuncture deqi behavior scores during acupuncture manipulation.

Electroencephalography Microstates Reveal the Neuromodulation Effects of Acupuncture With Deqi
Previous studies have stated that microstates are atoms of thought, which can represent basic cognitive processes in different states and tasks (Koenig et al., 2002;Brodbeck et al., 2012;Milz et al., 2016). To the best of our knowledge, this study was the first to propose the idea of EEG microstates to study the neuromodulation mechanism of acupuncture. And our results showed that EEG microstates could be used as novel markers to reveal the neuromodulation effects of acupuncture with deqi. First, microstate C was significantly increased during acupuncture stimulation (Figures 5, 6). On the other hand, microstate D was significantly decreased during acupuncture stimulation (Figures 5, 6). Second, the transition probabilities with microstate C as nodes were significantly increased (Figure 8 and Supplementary Figure S2), FIGURE 9 | Correlation analysis between neural responses and acupuncture's behavior scores across all subjects. Pearson's correlations between the transition probability parameter B→C and the deqi behavior index for the acupuncture manipulation and tactile controls, respectively. Each red dot is indicated for each subject's data. Dashed lines indicate the 95% confidence intervals for the mean of polynomial evaluation. Solid lines represent linear regression. *p < 0.05, n = 16.
whereas the transition probabilities with microstate D as nodes were significantly decreased during acupuncture stimulation (Figure 8 and Supplementary Figure S2). Moreover, the microstate parameters (microstate C's duration and microstate B→C's transition probability) showed a significantly positive correlation with acupuncture deqi behavior scores (Figures 7, 9). Taken together, the EEG microstates C and D could provide novel neuroimaging evidence for revealing the neuromodulation effects of acupuncture with deqi.

Acupuncture Manipulation Could Cause Significantly Increased Deqi Behavior Effects
The acupuncture deqi sensation that patients experienced could be a valuable behavioral index to reflect the acupuncture effect according to the TCM (Hui et al., 2007;Yu et al., 2012). Behavior results showed that the subjects' deqi indexes of the acupuncture manipulation were significantly larger than those of pre-/post-manipulation tactile control conditions. And there was no significant difference between pre-manipulation and post-manipulation tactile controls (Figure 3). These results indicate that acupuncture manipulation could effectively cause deqi behavioral effects as compared with control conditions. Therefore, this experimental design is effective and usable.

Electroencephalography Microstate Maps and the Corresponding Functional Networks
By comparing our results with previous classical EEG microstate maps, the observed four EEG microstate topography maps (A, B, C, and D) during acupuncture were in line with previous findings (Koenig et al., 2002; Figure 4A). Moreover, previous EEG-fMRI simultaneous recording studies showed that the four EEG microstate topography maps were related to different brain functional networks Musso et al., 2010;Van de Ville et al., 2010). Microstates A and B were related to the auditory and visual networks, respectively ; and microstates C and D were related to the salience network Michel and Koenig, 2018) and DAN, respectively (Smith et al., 2009;Britz et al., 2010;Michel and Koenig, 2018). In this study, the four EEG microstate topography maps during acupuncture manipulation showed the most dramatic changes with the lowest spatial correlation with pre-acupuncture rest control ( Figure 4B). We speculated that the corresponding functional networks could be modulated by acupuncture manipulation. Furthermore, microstates C and D had the smallest correlation values (Figure 4B), which may mean that the corresponding salience network and DAN were modulated by acupuncture manipulation. In addition, the postacupuncture rest control maps recovered back to the baseline condition and showed the overall highest concordance with pre-acupuncture rest control (Figure 4B), which may indicate that the corresponding functional networks were modulated by acupuncture manipulation and then recovered back to the baseline resting state.

Microstate C Could Be an Objective Marker for Acupuncture Manipulation
The acupuncture deqi sensation that patients experienced was a valuable behavioral index to reflect the acupuncture effect according to the TCM (Hui et al., 2007;Yu et al., 2012). However, the current behavioral evaluation of deqi is subjective, which needs an objective marker. Here, first, the EEG microstate C's parameters (duration, occurrence, and contribution) showed the most obvious and significant increase during acupuncture manipulation, compared with other controls (Figures 5, 6). This indicates that microstate C could be a potential marker for evaluate neuromodulation effects of acupuncture with deqi. Second, compared with tactile controls, the duration parameter of microstate C was significantly positively correlated with the acupuncture deqi behavior index during acupuncture manipulation (Figure 7). It meant that microstate C's duration parameter could reflect subject's deqi behavior index. Moreover, the previous study showed that the EEG microstate C had significant changes for adolescents with 22q11 deletion syndrome compared with normal control, indicating that microstate C could be a marker for the early diagnosis of schizophrenia (Tomescu et al., 2014). Therefore, we speculated that our findings also prove that the EEG microstate C could be an objective marker to evaluate neuromodulation effects of acupuncture with deqi.
Acupuncture Modulated the Salience Network Corresponding to Microstate C Previous EEG-fMRI simultaneous recording studies showed that the EEG microstate C was associated with the neural activities of brain regions including the right anterior insula, the anterior cingulate, and bilateral inferior brand gyri . Therefore, the observed microstate C's increase indicated that these underlying brain regions were engaged during the acupuncture manipulation, which was consistent with the previous fMRI finding of the insula and anterior cingulate cortex's activations during acupuncture (Fang et al., 2009;Chae et al., 2013). Furthermore, the anterior insula and anterior cingulate regions corresponding to microstate C belong to the salience network , which indicates that the salience network could be modulated during acupuncture stimulation.
In addition, it has been shown that the salience network could coordinate the neural resources to important self-related cognition and external stimuli of importance (Menon, 2011;Palaniyappan and Liddle, 2012;Uddin, 2015), by mediating the switch between the default mode network (DMN) and the central executive network (CEN) (Uddin, 2015). Therefore, it is reasonable to speculate that by activating the salience network, the external acupuncture stimulation helps to balance the activity of DMN and CEN and finally contributes to coordinate the self-related (internally directed) and goal-oriented (externally directed) cognition (Uddin, 2015). This helps to explain the previous finding that acupuncture at the Hegu acupoint could improve the cognitive ability of Alzheimer's disease patients by modulating the connectivity of large-scale network (Liang et al., 2014).

Acupuncture Modulated the Transition Probabilities With Microstate C as Node
Compared with those of controls, transition probabilities with microstate C as node were significantly increased during acupuncture manipulation (Figure 8). Especially, the transition probability of microstate B→C was significantly positively correlated with the deqi index during acupuncture manipulation (Figure 9). These results add more evidence for microstate C as an objective marker for acupuncture manipulation with deqi.
In addition, previous EEG-fMRI simultaneous study showed that EEG microstates could reveal the dynamic patterns of EEG activities (Lehmann et al., 1987;Khanna et al., 2014;Michel and Koenig, 2018), and especially transition probabilities could speculate the dynamic switching patterns of the underlying functional networks . And microstates A, B, and C were related with auditory network, visual network, and salience network, respectively, by previous EEG-fMRI study (Smith et al., 2009;Britz et al., 2010;Michel and Koenig, 2018). Therefore, we speculate that acupuncture manipulation could enhance the dynamic information exchanges between salience network and auditory/visual sensory networks, thus improving the ability for cognitive information processing (Brodbeck et al., 2012), which could be further evidenced by microstate C's significantly increased transition probabilities with both microstates A and B (Figures 8A,B). This could further explain previous findings that acupuncture stimulation could strengthen brain connectives (Yu et al., 2017(Yu et al., , 2018. Furthermore, the salience network was found to mediate neural resources to self-related cognition and external stimuli of importance (Menon, 2011;Palaniyappan and Liddle, 2012;Uddin, 2015). Therefore, it is reasonable to speculate that by enhancing the connections between salience network and sensory network, acupuncture could help to coordinate the internally or externally directed cognition. Taken together, the significantly increased transition probabilities with microstate C as node not only add more evidence for microstate C as an objective marker but also help to reveal the exchanging information between different functional networks during acupuncture manipulation.

Acupuncture Manipulation Modulated Microstate D and Its Corresponding Dorsal Attention Network
Compared with those of controls, the EEG microstate D's duration, occurrence, and contribution parameters showed significant decrease during acupuncture manipulation (Figures 5, 6). Here, we argue that microstate D could also be a marker to evaluate neuromodulation effects of acupuncture with deqi. Moreover, previous EEG-fMRI simultaneous study showed that the EEG microstate D was related with rightlateralized frontal and parietal cortex . And the frontoparietal cortex corresponding to microstate D belongs to the DAN (Seitzman et al., 2017;Michel and Koenig, 2018), which is responsible for the cognitive selection of sensory stimuli and responses (Corbetta and Shulman, 2002;Fox et al., 2006;Kim, 2014). Therefore, we speculate that the DAN is inhibited during the acupuncture manipulation, which is evidenced by the decrease of microstate D. In addition, the Hegu acupoint is considered to have analgesic effects (White et al., 2003;Kong et al., 2005) and is commonly used in the treatment of painful syndrome (Andersson and Lundeberg, 1995). We further postulate that by inhibiting the DAN, acupuncture manipulation could decrease the top-down attention to perceive pain sensory information, which could explain the analgesic effects of acupuncture at the Hegu acupoint. Taken together, we argue that the EEG microstate D could also be an objective marker to evaluate neuromodulation effects of acupuncture with deqi.

Acupuncture Modulated the Transition Probabilities With Microstate D as Node
Compared with those of controls, the transition probabilities with microstate D as node were significantly decreased during acupuncture manipulation (Figure 8), which adds more evidence for microstate D as an objective marker for acupuncture manipulation. In addition, previous EEG-fMRI simultaneous studies showed that microstate D was associated with DAN (Smith et al., 2009;Britz et al., 2010;Michel and Koenig, 2018). We further speculate that acupuncture could inhibit the information processing of DAN with the auditory and visual sensory network, which could be evidenced by microstate D's decreased transition probabilities with both microstates A and B (Figures 8A,B). By inhibiting the connections of DAN with other sensory networks, acupuncture manipulation could help reduce the subject's top-down attention to perceive sensory stimuli, which could further explain the analgesic effects of acupuncture (Hui et al., 2005;Chae et al., 2013). Taken together, the significantly decreased transition probabilities with microstate D as node not only add more evidence for microstate D as an objective marker but also help to reveal the exchanging information between different functional networks during acupuncture manipulation.

Limitations and Future Works
In addition, there were several potential limitations in this study. First, more participations should be recruited to provide more evidence for further exploring the correlation relationship between neural responses with behavior scores in the following study. For example, although there existed statistical significance, the statistical power and significance of the correlation results between behavior scores with microstate C's duration (Figure 7) and B→C's transition probability (Figure 9) parameters could be further enhanced and validated by recruiting more participations. Second, the temporal dynamics modulation of functional networks during acupuncture revealed by current EEG microstate analysis should be further explored and validated using the fMRI-EEG simultaneous recording method. The fMRI-EEG simultaneous recording could add more evidence for explaining the underlying functional networks changes revealed by different EEG microstates during acupuncture. Third, this study only explored the effects of deqi with EEG microstate mapping on healthy subjects. Future clinical studies could be conducted on patients to further validate the marker value of EEG microstates and to provide more neuroimaging evidences for further explaining the therapeutic effects of acupuncture with deqi. In addition, although the deqi index could characterize the significant difference between different conditions, the behavior questionnaire could be further optimized. Fourth, only the basic acupuncture unit, the Hegu acupoint, was used as an example to explore the neuromodulation effects of acupuncture with deqi. Future studies should further explore the neuromodulation mechanism by EEG microstates with other acupoints. Finally, magnetoencephalography (MEG) experiments could further add more evidence to explain the dynamic processing characteristics of the brain network during acupuncture manipulation.

CONCLUSION
Despite acupuncture's public acceptance and long history, how acupuncture manipulation modulates the human brain remains largely unknown. By using EEG microstates theory and multichannel EEG recordings under acupuncture manipulation and tactile/rest controls, our results showed that EEG microstates could be used as objective markers to reveal the neuromodulation effects of acupuncture with deqi. Compared with that of controls, acupuncture manipulation could significantly increase the duration, occurrence, and contribution parameters of microstate C, whereas it could decrease those parameters of microstate D. Besides, acupuncture manipulation could also significantly improve the transition probabilities with microstate C as node, whereas it could reduce the transition probabilities with microstate D as node. Moreover, the microstate parameters (microstate C's duration and microstate B→C's transition probability) showed a significantly positive correlation with acupuncture deqi behavior scores. Taken together, EEG microstates could be effective markers to provide novel evidence for explaining acupuncture's neuromodulation effects on human central nervous system, which also contributes to the objectification of TCM.

DATA AVAILABILITY STATEMENT
The data and codes that support the findings of this study are available upon request.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Institutional Review Board and Ethics Committee of Tianjin University. The patients/participants provided their written informed consent to participate in this study.