Participants and clinical characteristics assessment
In total, this study included 43 patients with CTN matched for age, gender, education with 45 healthy controls (HCs). All patients were recruited from Lanzhou University Second Hospital and diagnosed according to the International Classification of Headache Disorders (ICDH-III)[2] by two experienced neurologists. NVC were demonstrated for all patients on either MRI or during surgery. Inclusion criteria were: (1) the duration of disease was more than 1 years; (2) unilateral pain in the area innervated by one or more branches of the trigeminal nerve; (3) paroxysmal pain, described as electric shock-like, shooting or stabbing experience, was activated by trigger factors or in trigger areas; and (4) without obvious sensory deficits. Exclusion criteria were: (1) CTN with concomitant continuous pain; (2) had a previous surgical history, especially microvascular decompression (MVD) for CTN, or a history of head trauma; (3) suffered any other pain disorders or neuropsychiatric disease, (4) with metal implanted in the body, particularly metallic fixed dentures; and (5) abnormal MR manifestation (including severe white matter lesions with Fazekas grade 3 and evidence of multiple sclerosis or space-occupying lesions that indicate secondary TN)[9]. CTN patients all received medical treatment, of which the most is carbamazepine, the rest are mecobalamin or other medicine. No control measure about patients’ treatment were taken. The research was approved by the Ethics Committee of Lanzhou University Second Hospital. According to the Declaration of Helsinki, written Consent was obtained from every participant after introduction about project for details.
The pain intensity of CTN patients was recorded with visual analogue scale (VAS) ranging from 0 (no pain) to 10 (extreme amount of pain). Patients were required to rate their pain intensity for the last 7 days by marking on the 100mm line of questionnaire, and then averaged score for a week. All questionnaires assessment were performed under the supervision of experimenters.
Magnetic resonance imaging data acquisition
Structural MRI and rs-fMRI data were acquired on a 3.0T Siemens Verio MRI system (Siemens Medical System, Erlangen, Germany) with an 8-channel head coil. During scanning, participants were instructed to stay awake and relaxed but to keep their eyes closed, with earplugs and foam padding used to attenuate noise and reduce head motion. High-resolution three-dimensional structural images were acquired using sagittal Magnetization Prepared Rapid Gradient echo sequence (field of view: 256×256 mm; matrix: 256 × 256; time of repetition = 1900 ms; time of echo = 2.93 ms; resolution = 1 × 1 mm; flip angle = 9ᵒ). The rs-fMRI images were acquired via an echo-planar imaging (EPI) sequence (180 volumes; 64 contiguous slices/volume; FOV: 192 × 192 mm; matrix: 64 × 64; spatial resolution = 3 × 3 × 3 mm; TR = 2000 ms; TE = 30 ms; flip angle = m90ᵒ). An experienced radiologist inspected the previous MR images of these participants to make sure that each patient was free with abnormalities as described in above exclusion criteria (5).
fMRI data preprocessing and head motion analysis
The rs-fMRI data were preprocessed with the toolbox for Data Processing & Analysis of Brain Imaging (DPABI, http://rfmri.org/dpabi)[29]. The first ten volumes of the functional images were discarded. The remaining volumes underwent slice-time correction, then were realigned to correct the motion between time points. Wherein, head motion parameters were computed by estimating the translation in each direction and the angular rotation on each axis for each volume. As a result, the participants with mean framewise displacement (FD) (Jenkinson) larger than 0.2 mm or head displacement more than 1.5 mm, maximum rotation greater than 1.5° were excluded from the analysis. According to this exclusion criterion, two subjects of HC and two subjects of CTN were excluded. No significant intergroup differences were found in FD (t = 1.09, p = 0.28). The individuals fMRI data were co-registered to their structural images, followed by segmentation of the gray matter (GM), white matter (WM), cerebrospinal fluid (CSF) and normalization to the Montreal Neurological Institute (MNI) space. The normalized images were spatially smoothed with a 6 mm full-width at half-maximum Gaussian kernel.
GICA analysis and identification of independent components
After preprocessing, fMRI images with 170 volumes underwent group independent component analysis (GICA) to be decomposed into different RSNs by using the Group ICA of fMRI Toolbox software (version 4.0b; mialab.mrn.org/software/gift/)[30]. Two data reduction steps were performed in the principal component analysis[14]. First, we reduce the individuals data into 120 principal components, which preserved more than 99% of the variance. Next, we concatenated the reduced data of all participants across time and further reduced the data to 100 principal components using an expectation maximization algorithm[31]. The reliability and stability of the infomax ICA algorithm[32] was ensured by iterating 20 times in ICASSO implemented in GIFT[33] and using the most central run to reconstruct subject-specific time courses and spatial maps of each IC using the GICA back reconstruction algorithm[34]. The group ICs of the 20 runs were clustered to estimate their reliability, with values more than 0.8 were selected[35]. Through one sample t test across all subjects and for each IC, the t-map of ICs was obtained with a threshold of t > mean (µ) + 4SD (σ)[36]. Details about labels and spatial maps of each IC are presented in Figure S2, and the peak coordinates of ICs are shown in Table S1.
We identified 59 ICs from 100 ICs based on the following evaluation criteria: (1) IC should exhibit peak activations in grey matter; (2) low spatial overlap with known vascular, ventricular, motion, and susceptibility artifacts; (3) should have time courses dominated by low-frequency fluctuations (ratio of powers below 0.1 Hz to 0.15–0.25 Hz in spectrum)[37]. All 59 ICs were then sorted into different nine RSNs according to the spatial correlation values between their spatial maps and atlas used in previous studies[14] (Fig. 1A). Afterwards, additional postprocessing was applied to time courses of 59 ICs as described in Allen et al[14], including detrending, despiking using AFNI’s 3dDespike algorithm, filtering using a fifth order Butterworth filter with a 0.15 Hz high frequencies cut-off, and finally regressing out the movement parameters.
dFNC estimation
We computed dFNC between the time courses (170 time points) of ICs using a sliding window approach, which was performed using the DFNC network toolbox in GIFT. A window size with 20-TR (40s) was chosen because previous studies have suggested that FC fluctuations at resting-state would be captured with windows of 30 ~ 60s[38]. We used a tapered window in steps of 1 TR, which was obtained by convolving a rectangle with a Gaussian (σ = 3) to localize the dataset at each time point. Finally, a total of 150 windows were obtained and 59 × 59 pairwise FC matrices by regularized precision matrix (inverse covariance matrix)[39] were computed in every window. The L1 norm penalty was imposed in the Graphical LASSO framework with 100 repetitions to promote sparsity in estimations[40]. With Fisher's z-transformation, the correlation values of pairwise functional matrices were converted to z-values to improve normality and comparability and then residualized with nuisance variables, including age and gender[36].
State clustering analysis
To assess the dFNC patterns that reoccur over time, k-means clustering was performed on all FC matrices for all participants. The k-means clustering algorithm was iterated 100 times with L1 distance (Manhattan distance) function to estimate the similarity between matrices[41]. Later, the analysis for cluster number validity was made and determined optimal number of k as 2, based on the silhouette criterion[42], which was computed as a ratio of the similarity between windows in the same cluster compared to similarity with windows in a different cluster. In the next 100 clustering iterations, k = 2 was kept. Eventually, we obtained two reoccurring FC states, of which cluster centroids were determined as the median of all matrices allocated to that state over time. The subject-specific centroids of each state were calculated as median value similarly. Further, the subject-specific centroids belong to each group were averaged to get group-specific centroids for better visualization of group comparison patterns[43].
In order to describe the characteristics of the two cluster states, we mainly focused on the degree of global modularity. Modularity is a valuable measurement from graph theory in interpretation of dFNC states, because it evaluates both functional integration and segregation of networks[44]. Thus, we calculated modularity index Q for each state by using a normal Louvain community detection algorithm in Brain Connectivity Toolbox (www.brain-connectivity-toolbox.net/). Wherein, A larger Q represents a higher tendence of assigning ICs into different modules[45].
We calculated four temporal properties: fractional windows, mean dwell time, number of transitions and transition likelihood. The fractional window is calculated as the proportion of time spent in each state as measured by percentage. The mean dwell time represents average duration of time intervals an individual spent in each state, which was calculated by averaging the number of consecutive windows belonging to one state before switching to another. The number of transitions represents the switching times between states, which estimates the flexibility of brain. And the last one, transition likelihood, represents the percentage of switching probability between states. For between-group comparison of different properties, nonparametric permutation tests (10,000 repetition) were used to assess differences of all those temporal properties mentioned above, treating age and sex as covariates. False discovery rate (FDR) correction were applied for fractional windows and mean dwell time.
For the purpose of evaluating the consistency and validity of the k-means clustering at different window sizes, we repeated the dFNC states analysis with 16-TR (32s) and 24-TR (48s). By calculating Pearson's correlation coefficients between the cluster centroids under different window sizes to represent similarity and find the states consistent with the primary analysis[45].
Dynamic topologic analysis
We applied a graph theory approach to obtain topological metrics across all sliding windows using GRETNA software (www.nitrc.org/projects/gretna), so as to observe the variability of topological organization of the functional connectivity network. Based on the framework of graph theory, we defined the 59 ICs as functionally independent nodes with FC between pairs of ICs as edges. At first, FC matrices of all windows were binarized with a series of sparsity thresholds, where edges larger than threshold was designated as 1, and 0 when it was smaller than threshold; only positive FC values were considered. With regard to sparsity, it is defined as the ratio of the number of existing edges divided by the maximum possible number of edges in a network. Referring to previous studies[23, 28], we determined thresholds ranged from 0.10 to 0.35 (with an interval of 0.01) for further analyses.
Next, we calculated both global and regional network properties in a series of adjacent matrix for all participants. The former included: (1) measures of global (Eg) and local network efficiency (Eloc); and (2) small-world global metrics of clustering coefficient (C), characteristic path length (L), small-worldness (σ), normalized clustering coefficient (γ), and normalized characteristic path length (λ); and the later was nodal efficiency. As it has been widely used in previous studies, an AUC approach was chosen to avoid the specific selection of a threshold[23]. The detailed interpretation of topological properties is listed in Table S2. For better characterizing the temporal variation of those measurements, we also computed the coefficient of variation (CV) of AUC of network parameters as what was did in Luo et al[28], where CV was calculated as SD/Mean across all sliding windows.
The nonparametric permutation approach (10,000 iterations) was used again to test for dynamic topological property differences in the AUC of each metric with age and sex as covariates. When group differences in the larger number of nodal properties, a false discovery rate (FDR) correction was used to control the false-positive rate to one per analysis.
Correlational analyses
Given the fact that the dynamic measures obtained in our study were non-normal distribution, we performed Spearman's partial correlation analyses to investigate possible relationships between abnormal properties and clinical data (including illness duration, VAS and attack frequency). Demographics (age, sex, education) and head motion (FD Jenkinson) were regressed out and p < 0.05 was set as statistical significance threshold.