Participants and experiment protocol

This study was performed at the University of Michigan. All study participants gave written informed consent, and the study protocol and informed consent documents were approved by the University of Michigan Institutional Review Board (Ann Arbor, Michigan). We confirmed that all methods were performed in accordance with the relevant guidelines and regulations. 17 FM patients (all females, age mean ±SD: 45.7 ± 11.4) were recruited as part of a phase 2 clinical trial (NCT01914679) evaluating a novel noninvasive brain stimulation device (Reduced Impedence Noninvasive Cortical Electrostimulation, RINCE) as a treatment for chronic pain. For this study, we analyzed only baseline EEG data, before any treatments occurred. Ten out of seventeen patients’ data were used for the analysis. Seven subjects were excluded due to missing EEG data or high electromyogram or ocular movement artifact. Inclusion criteria were: FM patients who met the 1990 ACR diagnostic criteria for FM, female, and between the age of 18–65. Exclusion criteria were: current psychiatric disorder, history of suicide attempt in preceding 5 years, any history of bipolar disorder, schizophrenia or other psychotic disorder, a Hospital and Anxiety Depression score greater than 11 for either anxiety or depression, other chronic infection or condition that may cause pain, history of seizure disorder, pregnant or breastfeeding, history of alcohol or drug abuse, BMI greater than 40 kg/m2, or the use of stimulant medication, centrally active analgesics, or anesthetic or narcotic patches.

EEG recording

64-channel sensor net from Electrical Geodesics, Inc. (Eugene, OR) was used to acquire the EEG data with a sampling frequency of 500 Hz and the electrode impedance was kept below 50 KΩ.

The EEG protocol consisted of 10 minutes of resting state (5 minutes eyes open, 5 minutes eyes closed). Clinical pain was assessed immediately before the resting state period with a Visual Analog Scale (VAS) where “0” represented no pain and “100” represented the worst pain imaginable.

Noise and artifacts of the signals were automatically removed by the EEGLAB toolbox based on the power spectrum with a 10 dB threshold. In addition to the automatic rejection, we visually inspected the data and excluded artifacts. Independent component analysis (ICA) was also applied to the signals (function runica.m, EEGLAB, MATLAB toolbox, https://sccn.ucsd.edu/eeglab/index.php) and removed components of eye movements, cardiac signals, and focal noise. Noisy channels were also removed with a high-pass filter from 0.5 Hz, and 4-minute eyes closed artifact-free resting-state EEG signals were re-referenced to an average reference. We focused our analysis on the alpha frequency range (8–13 Hz) because it is the dominant peak frequency in the eyes-closed resting state.

EEG network analysis

To reconstruct the functional brain network, we used the weighted phase lag index (WPLI)26, which is a measure that captures phase locking between two signals with high robustness to volume conduction.

$$WPL{I}_{ij}=\,\frac{|E\{\Im ({C}_{ij})\}|}{E\{|\Im ({C}_{ij})|\}}=\,\frac{|E\{|\Im ({C}_{ij})|sgn(\Im ({C}_{ij}))\}|}{E\{|\Im ({C}_{ij})|\}}$$ (1)

where ℑ(C ij ) is an imaginary part of cross-spectrum C ij between two signals i and j. The cross-spectrum C ij is defined as \({Z}_{i}{Z}_{j}^{\ast }\), where a complex value Z i denotes the Fourier spectra of the signal i for a particular frequency and \({Z}_{j}^{\ast }\) is the complex conjugate of Z j . C ij can be represented as Re iθ, where R is magnitude and θ is the relative phase between signal i and j. If the phases of one signal i always lead or lag those of the other signal j, then WPLI ij equals 1. On the other hand, if the phase lead/lag relationship of two signals is perfectly random, the WPLI ij value is 0.

We segmented the EEG signals into 10-s epochs with 2-s overlap. WPLI matrix was obtained for each 10-s window. We then took the average of WPLI matrices over all windows. A binary adjacency matrix A ij was constructed using the top 40% of averaged WPLI values among all channel pairs. If the WPLI ij was included in the top 40% of WPLI values, then A ij = 1, otherwise, A ij = 0. We then calculated the node degree, which is the number of connections of each node in the network. We tested several thresholds: 30, 40, 50, 60, and 70%; the analysis results are robust below 50% (See Supplementary Fig. S1).

EEG network configuration for ES conditions

First, we calculated the Spearman correlation between node degree and frequency, which has been identified as one of the network conditions for ES18. For each 10-s window, we calculated the power spectral density (PSD) for all channels (“pwelch.m” in MATLAB, with 2-s Hamming windows, 1-s overlaps, and frequency resolution = 0.2 Hz). At each window, we obtained the frequency that has median power within a frequency band. We took the average of the median frequency over all windows to calculate the correlation between node degree and frequency.

Next, the frequency difference Y ij and the frequency assortativity A f between connected nodes were calculated to investigate the functional network conditions for ES19,21. The frequency difference Y ij is defined as

$${Y}_{ij}=\frac{|{Y}_{i}-{Y}_{j}|}{{Y}_{i}+{Y}_{j}}$$ (2)

where Y i and Y j are the median frequencies of the signal i and j, respectively. We used the average of frequency difference among connected nodes, <Y ij > to evaluate the frequency difference of one individual. The frequency assortativity A f is a Spearman correlation between the node frequencies and the average frequencies of connected neighbors in a network. If A f < 0, the higher frequency nodes tend to link with lower frequency nodes, or vice versa, which indicates that the frequency relationship between nodes in the network is diassortative. Furthermore, Spearman correlations were calculated to investigate the relationship between clinical pain scores and network sensitivities represented by the degree-frequency correlation, frequency difference and frequency assortativity for all subjects.

Lastly, we calculated the Spearman correlation between the pain score and the degree/median frequency of each channel for all patients. This was performed to determine which brain areas were associated with the pain score. We calculated z-scores of degree/median frequency over all channels for each subject. With the z-scores of each channel for all patients, we calculated the correlation between the pain scores and those z-scores. We considered a p-value less than 0.05 as statistically significant.

Human brain network models

A simple coupled oscillator, the Kuramoto model, was used to simulate the interaction between brain areas in the human brain network constructed with DTI of 82 nodes including cortical and subcortical areas27.

$$\frac{{\rm{d}}{\theta }_{i}}{{\rm{d}}t\,}={\omega }_{i}+{\sum }_{j=1}^{N}{\lambda }_{ij}{A}_{ij}\,\sin ({\theta }_{j}-{\theta }_{i})\,$$ (3)

where N is the total number of oscillators, θ i and ω i are the phase and natural frequency of the oscillator i, respectively, i = 1, 2, …, N at time t. For simplicity, we assumed that the coupling strength between oscillators i and j, λ ij is constant and uniform, i.e., λ ij = λ. A ij is the connection matrix of the human DTI anatomical network. A ij = 1 if oscillators i and j are connected, and A ij = 0 if they are not. Each oscillator i is assigned a random initial phase θ i , uniformly distributed between [−π, π] and a random initial frequency ω i , drawn from some arbitrary distribution function g(ω). The oscillators become spontaneously synchronized if λ is larger than a certain critical value.

We designed the natural frequency configurations g(ω) to embody the network conditions that were consistent with the empirical results of patients with varying pain scores. Specifically, empirical data were used to generate ES and non-ES networks as follows:

i) We applied the perturbations to the highest degree nodes (1, 5, 10, 15, and 20). Here we only present the significant perturbation results of the 15 highest degree nodes (~18% of all nodes), which includes brain regions such as the putamen, insula, precuneus, thalamus, superior frontal and superior parietal regions27. The perturbation results and the names of the brain regions for the other highest degree nodes are presented in Supplementary Fig. S2. ii) The values of ω i that were not involved in the 15 highest degree nodes were obtained from the Gaussian distribution function \(g(\omega )={G}_{\bar{\omega },{\sigma }_{\omega }}(\omega )\) with the mean \(\bar{\omega }=10\) Hz, deviation σ ω = 0.2 Hz, mimicking the alpha rhythms of human EEG activity28. iii) The 15 highest degree nodes were divided by two subgroups α and β in which oscillators have higher and lower frequency ranges than average – i.e., \({G}_{\overline{{\omega }_{\alpha }},{\sigma }_{{\omega }_{\alpha }}}({\omega }_{\alpha })\) with \(\overline{{\omega }_{\alpha }}=10.2\) Hz, \({{\rm{\sigma }}}_{{\omega }_{\alpha }}=0.02\) Hz and \({G}_{\overline{{\omega }_{\beta }},{\sigma }_{{\omega }_{\beta }}}({\omega }_{\beta })\) with \(\overline{{\omega }_{\beta }}=9.8\) Hz, \({{\rm{\sigma }}}_{{\omega }_{\beta }}=0.02\) Hz for each subgroup, respectively − causing large frequency differences as well as forming a V-shape relationship between the node degree and frequency in the network, which facilitates ES19. iv) With the above conditions, we carried out the simulation of 100 configurations with positive degree-frequency correlation, frequency difference >0.01, and frequency assortativity <−0.25 as the high pain score condition, determined by empirical observation of the participants who had a pain score >42. We refer to those network conditions as the ES conditions. v) We also generated the frequency configurations that have Gaussian distribution function \(g(\omega )={G}_{\bar{\omega },{\sigma }_{\omega }}(\omega )\) with the mean \(\bar{\omega }=10\) Hz, deviation σ ω = 0.2 Hz. The 100 configurations with the frequency assortativity >−0.25 were considered to be the low pain score conditions. We deem these as non-ES conditions.

The collective dynamics of the ensemble of the oscillators were measured by the order parameter,

$$z(t)=r(t){e}^{-i{\psi }(t)}\equiv \frac{1}{N}{\sum }_{j=1}^{N}{e}^{-i{\theta }_{j}(t)}$$ (4)

where ψ(t) is the average global phase. The modulus r(t) = |z(t)| so-called order parameter represents the degree of synchrony, being equal to 0 when the oscillators’ phases are uniformly distributed in [0, 2π) and 1 when they all have the same phase. The level of phase synchronization is determined by a time average of the order parameter after a transient period T Δ = 5000, R linked = 〈r(t)〉 T , with the whole time period T = 10000. We observe the state of the network by increasing the coupling strength λ by δλ = 0.002, from λ = 0.

Network sensitivity analysis

We introduced a frequency perturbation into the Kuramoto model slightly below the critical point to simulate external stimuli. We used the pair correlation function C p , which has been used to quantify the susceptibility of statistical physics models29,30,31 to measure network sensitivity at a global level as well as to find the critical point of the network. The pair correlation function in the Kuramoto model is defined as,

$${C}_{p}=N\{{\langle R{e}^{2}[z(t)]\rangle }_{t}-{\langle Re[z(t)]\rangle }_{t}^{2}\},$$ (5)

where N and z(t) are the total number of oscillators and complex order parameter, respectively. Re[z(t)] is the real part of z(t) in Eq. (4).

The pain-related brain regions such as insula, precuneus, superior frontal cortices, parietal cortices and the thalamus4,5,32,33 were perturbed as the target sites of the human brain network. The median frequency of the network frequency configuration was given to the target sites to facilitate ES.

The network sensitivity Δ(C p ) was defined as the absolute difference between the pair correlation functions at the critical coupling strength λ c , before and after the frequency perturbation as follows,

$${\rm{\Delta }}({C}_{p})\equiv {|{C}_{p}({s}_{p})-{C}_{p}({s}_{0})|}_{\lambda ={\lambda }_{c}}$$ (6)