From May 2020 to August 2023, 40 patients with NDPH and 43 healthy controls (HCs) were recruited from the Headache Department, Neurology Centre, Beijing Tiantan Hospital, Capital Medical University. Each recruited patient needed to be diagnosed with NDPH by two specialist neurologists. The inclusion criteria of NDPH group: (1) Satisfied the diagnostic criteria of NDPH according to ICHD-3 criteria; (2) Ages 18 to 70 years; (3) None of the patients enrolled had been prophylactically treated for NDPH for at least 3 months. The exclusion criteria of NDPH group: (1) Combined with other types of primary headache or major systemic diseases; (2) Inability to complete MEG and MRI (e.g., claustrophobia or metal implants in the body); (3) Poor data quality; (4) Pregnancy or breastfeeding. The same exclusion criteria were used for the age- and gender-matched HCs, who had no history of headache and were free of physical and psychiatric disorders. Headache information (headache history, headache frequency, etc.) and clinical scales were collected before MEG acquisition. The clinical scale included Headache Impact Test -6 (HIT-6), Patient Health Questionnaire-9 (PHQ-9), Generalized Anxiety Disorder-7 (GAD-7), Pittsburgh Sleep Quality Index (PSQI), Visual Analogue Scale (VAS) and Montreal Cognitive Assessment (MoCA). The above scales assessed the intensity of headache impact, anxiety and depression symptoms, sleep quality, pain degree, and cognitive level of patients.

The study protocol was approved by the Institutional Review Committee of Beijing Tiantan Hospital of Capital Medical University (KY2022-044), which was registered on the https://www.clinicaltrials.gov (unique identifier: NCT05334927). All participants provided informed written consent according to the Declaration of Helsinki.

All participants were imaged with a 3.0 Tesla MR scanner (GE Healthcare, Milwaukee, WI, USA) at the Nuclear Medicine of Beijing Tiantan Hospital. Examinations were performed by a neuroradiologist who was unaware of the participant's diagnosis. Participants were asked to keep their heads and neck still, stay awake, and close their eyes, with tools to reduce noise and head movements. After checking the images, exclude images with quality problems. T1-weighted volumetric images were obtained by the 3D BRAVO sequence (coronal acquisition, the field of view (FOV) = 256 mm, acquisition matrix = 256, slice number = 192, flip angle = 15°, TR = 850 ms, TE = 320 ms, voxel size = 1 × 1 × 1.5 mm3).

The Elekta Neuromag 306-channel scanner (Elekta TRIUX ®) was used in this study to record neural activity at 2000 Hz with a low-pass filter set to 660 Hz. The Elekta Neuromag scanner with 306 channels (102 magnetometers and 204 gradiometers) was used in this study. The head position of the participants is detected by four (head position indicator) HPI coils, and re-recorded if the head movement was excessive during the scan. Then, the head position is digitally marked (Polhemus Fastrak®). About 300 points were marked on the nasion, anterior points in front of the ear points and scalp for MRI co-registration. Data acquisition was performed using a 2000 Hz sampling rate and a low-pass filter set to 660 Hz, while the participants’ electrooculogram and electrocardiogram were recorded. Resting-state MEG data were collected for five minutes for each participant. During the scan, participants were instructed to keep their heads and neck still, stay awake, and close their eyes.

After checking and excluding bad channels in the original data, the data were filtered using MaxFilter. The data sampling rate was reduced to 1000 Hz and 50 Hz line noise was removed, and finally a bandpass filter of 1-80 Hz was applied. In the independent component analysis (ICA), ocular and cardiac artifacts were marked and excluded. The cleaned data were then used to construct functional connectivity networks at the source level. Neural activity was filtered into five frequency bands: delta (1-4 Hz), theta (4-8 Hz), alpha (8-13 Hz), beta (13-30 Hz), gamma (30-80 Hz). At the source level, MRI images and MEG data were registered, and calculated the surface-based source space and inverse solution. Finally, dynamic statistical parametric mapping (dSPM) was used for source estimation [14].

Compute envelope correlations of orthogonalized activity as FC using pairwise and symmetric orthogonalization in source space. The procedure for symmetric orthogonalization in is: Extract inverse label data from raw; Symmetric orthogonalization; Band-pass filter; Hilbert transform and absolute value. According to the Desikan-Killiany Atlas, the brain was divided into 68 cortical regions (nodes). This power envelope time course was then correlated between brain region for each individual. These communication links between cortical regions (nodes) correspond to “edges” in a graph theory network model, and this specific approach for connectivity estimation shows greater repeatability than a wide range of other choices (Fig. 1) [15]. Network-based statistic (NBS) analysis was used to investigate which functional connections were significantly different between NDPH and HCs.

The pipeline of neural physiological signal to the construction of the brain network

We computed the network topological parameters with graph theory that summarize the aspects of segregation and integration of a network. In this study, we focused on network nodal and global parameters including nodal clustering coefficient, nodal efficiency, nodal degree, global efficiency, local efficiency, shortest path length of a network. The complex network analyses were performed at a sparsity range from 0.05 to 0.39 with an interval of 0.01, and the area under curve (AUC) values under this range of sparsity were calculated for both global and nodal network parameters for statistical analyses. All network analyses were performed by Gretna (http://www.nitrc.org/projects/gretna/) and visualized by using BrainNet Viewer software (http://www.nitrc.org/projects/bnv/).

The sample size was determined based on the available data and previous literature. Assuming no negative correlation between endpoints, a sample size of 75 cases (40 HC group and 35 NDPH group) would provide 90% power to reject the null hypothesis equal means at a two-sided alpha of 0.05. IBM SPSS 26.0 was used to perform the statistical analysis. Independent sample t-test and chi-square test were used to calculate the statistical difference of clinical data between groups. Continuous variables were tested for normality using the Kruskal–Wallis test; data conforming to a normal distribution were expressed as mean ± standard deviation, and otherwise as median with interquartile range. Categorical variables were expressed as numbers (percentages). Correlations between network parameters and demographic data were calculated using Pearson’s correlation, and false discovery rate (FDR) correction was performed. In NBS analysis and graph theory analysis, age and gender were used as covariables. Finally, the patients was grouped according to the headache location (unilateral and bilateral), the network topology parameters between the groups were compared. All statistical tests were two-tailed tests, with P < 0.05 indicating statistical significance.