The chiropractor employed typical manual or instrument-assisted interventions that are common in the chiropractic profession (Fig. 1). The chiropractor performed either manual high-velocity, low-amplitude (HVLA) adjustments or instrument-assisted HVLA adjustments using an Activator instrument™ or a drop-piece table (Fig. 1) [41]. The Activator instrument™ (Activator Methods International, Ltd., Phoenix, AZ, USA) used was a hand-held device that delivers an HVLA thrust, which can be set at various pre-determined force levels, directed at dysfunctional spine or pelvic joints [26]. The activator instrument has previously been shown to activate a similar, but smaller, neural response as compared to a manual HVLA thrust, and is capable of also altering the median nerve N30 SEP peak amplitude in younger subclinical pain populations [24]. The drop-piece table, which is a purpose-built adjustment table, also enables a rapid movement to be induced at specific spinal segments. The patient is positioned prone, and the appropriate drop pieces are lifted where it sits elevated, until the thrust is delivered, allowing the drop piece to drop down, aiding in the HVLA thrust to allow rapid movement between a motion segment, such as between the sacral base and L5 motion segment [26]. The chiropractor selected the site(s) for the spinal adjustment based on routine clinical indicators [17], which included tenderness to palpation of the relevant joints, restricted intersegmental range of movement on manual palpation, palpable asymmetric intervertebral muscle tension, and any unusual or blocked joint play and end-feel of the joints. The chiropractor adjusted multiple levels of the spine and pelvis depending on each patient’s clinical findings. The segments adjusted for each participant are noted in Table 1.

The control intervention acted as a physiological control for possible changes occurring due to the cutaneous, muscular, or vestibular input that would have occurred with the movements involved in preparing a patient for chiropractic spinal adjustments. The chiropractor applied a sham spinal adjustment by passively moving the subject’s head, spine, and body into positions approaching those used in the chiropractic adjustment intervention group. However, the chiropractor took care to not provide a manipulative thrust/impulse or to take a spinal motion segment to the end-range tension.

The resting-state EEG was recorded for five minutes. For artifact detection and correction, the preprocessed EEG was high-pass filtered using a Kaiser- windowed finite impulse response (FIR) filter ($\beta$ = 5.653) with a cutoff frequency of 1 Hz and an order of 4948 equivalent to the transition bandwidth of 1.5 Hz. The data from subjects with at least two minutes of clean EEG per session were used for further analysis. The analyses were performed on the data downsampled to 256 Hz to reduce the computational time and remove redundant frequency spectrum not required in the analysis. The EEG was segmented into two-second long epochs to encompass two cycles of the lowest frequency of interest (1 Hz). For any session with clean data longer than two minutes, 60 epochs were randomly selected to equalize the length of data across sessions and subjects.

SEPs are electrical responses generated by the nervous system in response to sensory stimuli, typically involving the stimulation of peripheral nerves. These potentials are recorded through electrodes placed on the scalp or over specific regions of the body. SEPs provide valuable information about the integrity and functioning of somatosensory pathways in the nervous system. When a sensory stimulus, such as electrical or tactile stimulation, is applied to a specific part of the body, the resulting electrical signals can be measured and analyzed. The primary components of SEPs include N20, P25, N30, and P45, each corresponding to different stages of sensory processing along the neural pathway. Various peripheral nerves, including the median, ulnar, radial, tibial, peroneal, femoral, sural, and superficial radial nerves, can be stimulated to evoke somatosensory evoked potentials, enabling a comprehensive assessment of sensory pathways in both upper and lower limbs.

The raw EEG was truncated to keep up to an additional 30 s of data at the beginning and end of the recordings to reduce filter artifacts. The early-stage EEG processing pipeline (PREP) [57] identified bad channels, removed line noise, and provided the average referenced data. For both resting-state EEG and SEPs, an epoch-based cleaning method was employed. Epochs for SEPs were extracted from 100 ms pre-stimulus to 150 ms post-stimulus. Baseline correction was applied using the pre-stimulus period. We also divided resting-state EEG into 0.5 s epochs. We marked epochs bad if, for any channel: (i) the amplitude was greater than 100 µV, (ii) peak-to-peak amplitude was more than 150 µV, (iii) the amplitude was greater than 100 µV in a step-function with a sliding window 200 ms wide with a step size of 50 ms, (iv) sample-to-sample difference was more than 50 µV, or (v) the amplitude was less than 1 µV for 150 ms (i.e., flat-lined data). All the epochs were also verified manually by visualisation to classify them as good or bad epochs. The step-like artifacts in the frontal channels, which were related to eye blinks and movements, were not removed [58].

The high-pass filtered EEG was down sampled to 512 Hz and epoched in a similar way as described above. The bad epochs and bad channels (from PREP) were removed from the data, and adaptive mixture independent component analysis was used to decompose EEG signals into independent spatial and temporal components [59]. This algorithm was chosen because it outperforms other independent component analysis (ICA) algorithms [60]. The ICA weights obtained were applied to the EEG data obtained after the PREP pipeline. A Kaiser-windowed FIR filter ($\beta$ = 5.653) with an order of 7420 and a transition bandwidth of 1 Hz was used for resting-state EEG. For SEPs, we used a 2nd order Butterworth filter. After performing ICA, the resulting independent components (ICs) were manually categorized into brain or non-brain components based on whether they were associated with muscle activity, eye activity, or noise from the channels or mains. The ICs were classified according to spatial distribution using scalp topography, time course, spectrograms, event-related potential (ERP) images, and current dipole models, considering recommended procedures [30, 36] and prior studies [61, 62], and guided by a website (https://labeling.ucsd.edu/). The marked ICs were removed, and the noisy channels were interpolated to provide a cleaned dataset for further processing.

The data were excluded if: (i) the semi-automated epoch classification approach mentioned above failed due to high levels of noise, making it difficult for the automatic methods to mark epochs as bad (e.g., all epochs marked bad), and it was impossible to visually verify whether the epochs marked are correct or not, (ii) the number of good ICs was zero after ICA classification.

Non-parametric cluster-based permutation tests [63] were used to evaluate the differences between interventions based on the global power spectrum of the resting state using a cluster threshold of 0.05. The clusters were identified as two or more channel-frequency pairs, each having a p-value of less than 0.05 from the paired t-test (two-tailed). The maximum of cluster-level statistics was determined by adding the t-values within each cluster. If the Monte Carlo probability for each tail exceeded the threshold of 0.025 compared to the reference distribution, which was approximated by the Monte Carlo method with 5000 permutations, the cluster was considered significant.

The LORETA-KEY software was used to perform a statistical procedure for source localization of resting-state EEG. Non-parametric mapping [64] was used, which utilized Fisher’s random permutation test with 5000 randomizations to control for the multiple comparison problem. A paired two-tailed t-test was used to find differences in the current sources across different frequency bands. The tests were used to identify differences between the pre-sessions, the post- and pre-chiropractic sessions, and the post- and pre-control sessions.

In this study, we employed GraphVar (Charité Universitätsmedizin, Berlin, Germany), a toolbox by [65], to analyse 14 $\times$ 14 PLI connectivity matrices. We examined the connections between each pair of nodes within the matrix. To handle the challenge of multiple comparisons, GraphVar organized significant links into Graph-Components, which can be considered sub-networks. These components were measured in a manner similar to how clusters are identified in functional magnetic resonance imaging (fMRI) [66]. To assess whether the size of a Graph-Component was non-random, we compared it to randomly generated data within GraphVar. We computed a p-value for each non-random component. This allowed us to pinpoint significant connectivity patterns. In the statistical section of GraphVar, we used a within-subject design, where data from subjects were collected across multiple sessions (pre- and post). GraphVar calculated the mean PLI difference between two sessions simultaneously (e.g., Post - Pre). Importantly, these calculations only considered significant non-random graph components. The results highlighted the brain connections where the mean PLI significantly differed between sessions, helping us identify changes in connectivity patterns over time. The same analysis was run for both the chiropractic and control groups.

A mixed-effects model was used to identify the effects of the intervention on the SEPs, with intervention (control and chiropractic) and session (pre and post) as fixed effects and subjects as a random effect. The akaike information criterion corrected (AICc) values were used to evaluate the choice of the link function (identity or log). Based on the AICc values, we used a gamma distribution-based model with a log link for Alzheimer’s data and a linear model with an identity link for Parkinson’s data. The statistical procedure was performed in R software (version 4.0.4, Vienna, Austria) using lme4 package version 1.1.26 [67]. The contrasts were obtained using the emmeans package version 1.5.4 (https://CRAN.R-project.org/package=emmeans) [68], adjusted for multiple comparisons using Tukey’s honestly significant difference test (HSD).

Two-way repeated-measures analysis of variances (ANOVAs) were conducted, using intervention (control and chiropractic) and session (pre and post) as factors, to determine the effects of the treatment on brain source strengths. Pairwise comparisons were performed using Tukey’s HSD. The statistical analysis was conducted using MATLAB 2015b (The MathWorks, The MathWorks, Inc., Natick, MA, USA).

Alzheimer’s Disease: For all frequency bands, no significant differences in the power were seen between the pre-intervention sessions (one positive cluster p = 0.1170) (Fig. 3A), post- and pre-control intervention (no cluster) (Fig. 3B), and post- and pre-chiropractic spinal interventions (one negative cluster p = 0.1702) (Fig. 3C). Although not statistically significant, there was an increase in the grand-averaged absolute power across all frequency bands after the chiropractic intervention. After the control intervention, the absolute power was similar or lower in delta, beta, and gamma bands, and higher in theta and alpha bands, compared to baseline values. The difference between the two grand averages is shown in Fig. 3D.

Fig. 3.

Resting-state frequency-domain analysis (Alzheimer’s Disease). Scalp topographies of the difference in grand-averaged power comparing (A) the pre-chiropractic intervention to the pre-control, (B) the post-control to the pre-control session, (C) the post-chiropractic to the pre-chiropractic session, and (D) the difference between interventions [(C) minus (B)]. Although the cluster-based permutation tests showed no significant differences in comparisons (A), (B), and (C), the absolute power in all frequency bands was found to be higher after the chiropractic intervention.

Parkinson’s Disease: For all frequency bands, no significant differences in power were seen between the pre-intervention sessions (no clusters) (Fig. 4A), post- and pre-control intervention (no cluster) (Fig. 4B), and post- and pre-chiropractic spinal interventions (no clusters) (Fig. 4C). Although not statistically significant, there was an increase in the grand-averaged absolute power across all frequency bands after the chiropractic intervention, while after the control intervention, the absolute power was similar or lower compared to baseline values. The difference between the two grand averages is shown in Fig. 4D.

Fig. 4.

Resting-state frequency-domain analysis (Parkinson’s Disease). Scalp topographies of the difference in grand-averaged power comparing (A) the pre-chiropractic intervention to the pre-control, (B) the post-control to the pre-control session, (C) the post-chiropractic to the pre-chiropractic session, and (D) the difference between interventions [(C) minus (B)]. Although the cluster-based permutation tests showed no significant differences in comparisons (A), (B), and (C), the absolute power in all frequency bands was found to be higher after the chiropractic intervention.

In individuals with Alzheimer’s and Parkinson’s disease, the sLORETA analysis showed that the pre-sessions were similar, and neither of the chiropractic and control interventions had any significant effect on brain activity in any frequency band (all p $>$ 0.05).

Alzheimer’s Disease: The mixed model (Table 5) showed a significant effect of the chiropractic intervention for the N30 SEP peak amplitude (p = 0.005). The post-hoc analysis showed that the N30 amplitude was decreased by 15% after chiropractic spinal adjustment (post: 2.13 $\pm$ 1.89 µV, pre: 2.76 $\pm$ 1.92 µV, pre/post: p = 0.027; 95% confidence interval (95% CI) = [1.03, 1.75]). There was no change after the control intervention (post: 2.92 $\pm$ 1.82 µV, pre: 3.00 $\pm$ 1.75 µV, pre/post: p = 0.7876; 95% CI = [0.798, 1.35]). Fig. 6 (Ref. [69]) shows the distribution of N30 amplitude across the four sessions.

Fig. 6.

N30 SEPs amplitude (Alzheimer’s Disease). Dots represent the N30 amplitudes from all the analyzed subjects, while the boxplots display the median, 25th, and 75th percentiles. The error bars indicate the mean $\pm$ 95% CI. The distribution plots show the density distribution estimated using a Gaussian kernel with a SD of 1.5. The N30 amplitude was similar for the pre-intervention sessions. However, after chiropractic spinal adjustment (represented by the dashed black line), the N30 amplitude significantly decreased, but it remained unchanged after the control intervention (represented by the solid black line). The figure was created using a modified version of the code provided by Allen et al. [69]. 95% CI, 95% confidence interval; SD, Standard Deviation.

Parkinson’s Disease: Data from eight participants were analyzed as the N30 peak for one participant could not be identified. The mixed model (Table 6) showed no significant intervention effect and no main effects of intervention or session for the N30 SEP peak amplitude (p $>$ 0.05). The mean amplitudes for chiropractic group were: post = 3.19 $\pm$ 1.71 µV and pre = 3.03 $\pm$ 1.06 µV. The mean amplitudes for the control group were: post = 2.51 $\pm$ 1.41 µV and pre = 2.38 $\pm$ 1.21 µV. Fig. 7 (Ref. [69]) shows the distribution of N30 amplitude across the four sessions.

Fig. 7.

N30 SEPs amplitude (Parkinson’s Disease). Dots represent the N30 amplitudes from all the analyzed subjects, while the boxplots display the median, 25th, and 75th percentiles. The error bars indicate the mean $\pm$ 95% CI. The distribution plots show the density distribution estimated using a Gaussian kernel with an SD of 1.5. The N30 amplitude was similar for the pre-intervention sessions and did not show a significant difference following either the chiropractic spinal adjustment (represented by the dashed black line) or the control intervention (represented by the solid black line). The figure was created using a modified version of the code provided by Allen et al. [69].

The results of the LORETA analysis showed the presence of five separate regions during the 20–60 ms time frame following the stimulus: the primary somatosensory cortex (SI) on the opposite side, prefrontal cortex, cingulate, and both secondary somatosensory cortices (SII). Based on these findings, it was assumed that there were five sources and the dipoles were placed accordingly in these regions.

Alzheimer’s Disease: Data from 13 participants were included as residual variance for one participant was higher than the threshold (10%). The ANOVA showed no interactions between the interventions and sessions for source strengths (all p $>$ 0.05).

Parkinson’s Disease: Data from eight participants were analyzed as N30 peak of one participant could not be identified. The ANOVA showed no interactions between the interventions and sessions for source strengths (all p $>$ 0.05).