$$\rightleftharpoonup{xx}$$
$$\longleftharp{xx}$$,
$$\longrightharp{xx}$$,
In terms of recruitment results, participants were primarily recruited via mailing of recruitment letters and follow-up phone calls based on the outlined regulations of the Atlanta VA Healthcare System. The study team recruited a total of 50 participants, proving the effectiveness of the methods used in meeting the recruitment goal (see Figure 2). The use of the new clinical fibromyalgia diagnostic criteria allowed the study team to properly screen out individuals who did not meet fibromyalgia criteria32. Asking possible participants about a diagnosis of fibromyalgia is not as robust a measure as the additional screening and could have led to improperly scheduled baseline visits or study participation. Forty-eight participants were randomized into active and sham groups; two participants were excluded because they did not meet eligibility criteria based on the inclusion testing for the study.

Figure 2: Recruitment flow chart. A report and flow diagram of study recruitment, randomization, and allocation of intervention. Please click here to view a larger version of this figure.
Three outcomes were used in sample size calculations: DVPRS (clinical pain), 30 s chair stand test (function), and DMN-SMN connectivity (rs-fcMRI). All power calculations were based on preliminary data. Clinical pain changes using DVPRS were chosen as the primary outcome of interest. The 30 s chair stand test (30sCST) was chosen as the representative functional outcome since it exhibited the smallest between-group change. DMN-SMN connectivity was chosen as the secondary outcome of interest as a neuroimaging biomarker for clinical pain and treatment response. Sample size analyses were conducted assuming a significance of 1% and 80% power (2-sample, 1-sided).
The seeds for this analysis were chosen based on preliminary data, as well as the literature on fibromyalgia, pain, and CES16,24,26,63,64,65,66,67,68,69. Based on the preliminary data, the mean (SD) of the post-treatment change in left primary sensorimotor cortex (L-S1M1) to left posterior cingulate cortex (L-PCC) connectivity is 0.041 (0.079) for the treatment group and -0.026 (0.049) for the standard treatment group; the observed between-group difference effect size is 1.03 (see Table 216,21,26,62,63,64,70,71,72,73). The study would need 20 subjects in each group to achieve 80% power to detect the difference between the CES group and the standard treatment group in their post-treatment change in connectivity for L-S1M1 to L-PCC at the significance level of 0.01 using a two-sided t-test, assuming the between-group difference effect size is 1.03 as observed in the pilot data. Though the pilot study observed a 17% attrition rate for 12 subjects in our prior study of auricular neuromodulation (all completed their follow-up MRI, but two were lost to follow-up at 8 weeks and 12 weeks), in order to maintain a conservative estimate for sample size calculations, this study assumed a 20% attrition rate. With an expected 20% attrition at the post-treatment visit, the study needed to recruit 20/0.8 = 25 subjects per group.
| DMN seeds (x,y,z) | SMN seeds (x,y,z) | SN seeds (x,y,z) |
| Medial Prefrontal Cortex62,67 | Right Putamen64,71 | Right dorsolateral prefrontal cortex62 |
| Right PCC70 | Left M170 | Left anterior insula62 |
| Left PCC16,67 | Right M170 | Right anterior insula62 |
| Precuneus71 | Right S1-Hand16,72 | Left posterior insula64 |
| Left S1-Hand16,72 | Right posterior insula63 |
| Thalamus21 | Dorsal anterior cingulate cortex72,73 |
| | Right temporoparietal junction62 |
Table 2: Seeds for analysis. DMN, SMN, and SN seeds chosen for analysis based on a priori hypotheses. Each seed is presented with references to prior literature supporting its testing in pain syndromes.
Based on prior research of CES in civilian subjects with fibromyalgia, a total of 50 subjects (n = 25 sham, and n = 25 true) should achieve 80% power to detect a difference in pain scores between the two groups17 (see Table 3). Sample size calculations for this study were conducted using sealedenvelope.com and were based on preliminary data.
| Control | Intervention | N per group |
| Mean Change | SD | Mean Change | SD |
| Clinical Pain (DVPRS) | 0.375 | 1.493 | -1.833 | 2.229 | 10 |
| Function (30sCST) | -0.250 | 1.500 | 3.000 | 4.980 | 14 |
| rs-fcMRI (S1M1 to PCC) | -0.026 | 0.049 | 0.041 | 0.079 | 20 |
Table 3: Sample size calculation. Calculations related to study sample size.
Imported fMRIPrep functional data were smoothed using spatial convolution with a Gaussian kernel of 8 mm full width half maximum (FWHM) (see Figure 3) shows the functional output from fMRIPrep normalized into MNI152NLin2009cAsym template space (left) and the smoothed functional image from CONN Toolbox (right). This results in an increase to the signal-to-noise ratio, which in turn improves the detection of blood-oxygen-level-dependent (BOLD) signals.

Figure 3: Single subject comparing an unsmoothed functional image in MNI space (left) to its smoothed counterpart at 8 mm FWHM. Please click here to view a larger version of this figure.
The data were then denoised using a standard denoising pipeline74, including the regression of potential confounding effects characterized by white matter time-series (5 CompCor noise components), CSF time-series (5 CompCor noise components), motion parameters and their first order derivatives (12 factors)75, outlier scans (below 295 factors)48, and linear trends (2 factors) within each functional run, followed by bandpass frequency filtering of the BOLD timeseries76 between 0.008 Hz and 0.09 Hz. CompCor49,77 noise components within white matter and CSF were estimated by computing the average BOLD signal as well as the largest principal components orthogonal to the BOLD average, motion parameters, and outlier scans within each subject's eroded segmentation masks (see Figure 4). From the number of noise terms included in this denoising strategy, the effective degrees of freedom of the BOLD signal after denoising were estimated to range from 33 to 240.6 (average 173.4) across all subjects. The denoising resulted in a reduction of physiological and other extraneous noise from the data that could have resulted in confounding effects.

Figure 4: Quality checks. Quality checks graphs from CONN Toolbox displaying the effects of denoising on functional connectivity (FC), mean global signal and maximum motion. Note (A) upper most graph displays data for a single subject and session, (B,C) graphs B and C are the results of group-level denoising. Please click here to view a larger version of this figure.
Group-level analyses were performed using a General Linear Model (GLM)78. For each individual voxel, a separate GLM was estimated, with first-level connectivity measures at this voxel as dependent variables (one independent sample per subject and one measurement per task or experimental condition, if applicable), and groups or other subject-level identifiers as independent variables. Voxel-level hypotheses were evaluated using multivariate parametric statistics with random effects across subjects and sample covariance estimation across multiple measurements. Inferences were performed at the level of individual clusters (groups of contiguous voxels). Cluster-level inferences were based on parametric statistics from Gaussian Random Field theory62,79. Results were thresholded using a combination of a cluster-forming p < 0.005 voxel-level threshold and a familywise error corrected p < 0.0016 cluster-size threshold80. Following these steps results in group-level functional connectivity values comparing the CES and sham conditions based on regions of interest (ROI). These results can be visualized in a multitude of ways through the results explorer GUI within CONN Toolbox. To see a volume display visualization of an example group-level results with the red area indicating regions with greater positive connectivity with the ROIs and the blue are regions with greater negative connectivity with the ROI, see Figure 5.

Figure 5: Volume display of post cingulate seed. The image displays greater positive connectivity with the anterior cingulate gyrus and (red) and greater negative connectivity with the precuneus (blue) for the True condition relative to the Sham, pFWEc < 0.05. This figure represents partial group-level data (n = 34). Please click here to view a larger version of this figure.
Supplementary File 1: CES device order instructions Please click here to download this file.
Supplementary File 2: CES device log Please click here to download this file.
Supplementary File 3: fMRIPrep Boilerplate Please click here to download this file.
Supplementary File 4: CES CONN instructions Please click here to download this file.
Supplementary File 5: CES R code plots Please click here to download this file.
Supplementary File 6: R Code CES eddy-qc Anova Please click here to download this file.
Supplementary Figure 1: CONN workbook.(A) Second-level covariates. (B) Experimental conditions. (C) Denoising. (D) Results. (E) Structural data. Please click here to download this figure.