Effects of transcutaneous spinal stimulation on spatiotemporal cortical activation patterns: a proof-of-concept EEG study

Objective. Transcutaneous spinal cord stimulation (TSS) has been shown to be a promising non-invasive alternative to epidural spinal cord stimulation for improving outcomes of people with spinal cord injury (SCI). However, studies on the effects of TSS on cortical activation are limited. Our objectives were to evaluate the spatiotemporal effects of TSS on brain activity, and determine changes in functional connectivity under several different stimulation conditions. As a control, we also assessed the effects of functional electrical stimulation (FES) on cortical activity. Approach. Non-invasive scalp electroencephalography (EEG) was recorded during TSS or FES while five neurologically intact participants performed one of three lower-limb tasks while in the supine position: (1) A no contraction control task, (2) a rhythmic contraction task, or (3) a tonic contraction task. After EEG denoising and segmentation, independent components (ICs) were clustered across subjects to characterize sensorimotor networks in the time and frequency domains. ICs of the event related potentials (ERPs) were calculated for each cluster and condition. Next, a Generalized Partial Directed Coherence (gPDC) analysis was performed on each cluster to compare the functional connectivity between conditions and tasks. Main results. IC analysis of EEG during TSS resulted in three clusters identified at Brodmann areas (BA) 9, BA 6, and BA 4, which are areas associated with working memory, planning, and movement control. Lastly, we found significant (p < 0.05, adjusted for multiple comparisons) increases and decreases in functional connectivity of clusters during TSS, but not during FES when compared to the no stimulation conditions. Significance. The findings from this study provide evidence of how TSS recruits cortical networks during tonic and rhythmic lower limb movements. These results have implications for the development of spinal cord-based computer interfaces, and the design of neural stimulation devices for the treatment of pain and sensorimotor deficit.


Introduction
Approximately 300 000 Americans are currently living with spinal cord injury (SCI), with roughly 17 810 new cases added each year [1]. SCI is defined as any damage that impacts the spinal cord and the more rostral the location of the injury, the more bodily function is impacted [2][3][4]. Damage to the spinal cord can impact sensory, motor, and reflex capabilities below the injury site [5,6]. This can lead to secondary complications such as, muscle atrophy, pressure ulcers, cardiovascular, respiratory, urinary and bowel, neuropathic pain, or osteoporosis and bone fractures [3,[7][8][9]. It is estimated that the lifetime direct costs for a person with SCI who is 25 years old is between 2.1 and 5.4 million dollars [10]. Cost of care for the U.S. veterans administration alone exceeded $14.47 million or $21 450 per patient with SCI in 2005 [11].
Treatments for SCI currently focus on physical therapy, which is done to improve long-term outcomes, such as mobility, recovery of motor function, and independence [12]. Due in part to the limited success with physical therapy, several invasive surgical treatments have been explored, ranging from peripheral nerve grafts and more recently, electrical neuromodulation using techniques such as epidural spinal cord stimulation (ESS) [13,14]. ESS has been shown to be effective at restoring motor and autonomic functions in people with SCI [14][15][16][17][18][19]. While long-term studies with ESS for people with SCI are still in progress, the results are promising. Unfortunately, ESS is an invasive treatment, which carries high risks, such as infection and lead migration or other hardware problems [20]. Transcutaneous spinal cord stimulation (TSS) is a promising non-invasive alternative to ESS. Applying TSS at either the cervical or lumbar level of the spinal cord has been shown to improve upper and lower limb function in SCI populations [21][22][23][24][25][26].
Functional electrical stimulation (FES) artificially induces muscle contractions by applying electrical stimulation through electrodes placed on the skin over muscles directly, typically to restore lost function [27][28][29]. The strength of muscle contraction from the stimulation is a function of total electrical charge transferred to the muscle and depends on pulse waveform, duration, amplitude, and frequency [30] FES devices are currently being explored for restoring function after neurological injury such as stroke or SCI [30,31].
A key limitation to electrical neuromodulation is the identification of optimal stimulation parameters to induce recovery of volitional motor function. A more complete understanding of the ascending effects of TSS could help find better parameters for treatment of individuals with neurological injuries and disorders, such as SCI or stroke, by using stimulation characteristics to target particular regions and function of supraspinal neural circuitries. Once quantified, stimulation parameters could then be altered to better match cortical activity in different clinical populations, which may lead to improved outcomes [32]. Once quantified, stimulation parameters could then be altered to better match cortical activity in neurologically intact participants, which may lead to improved outcomes. Additionally, understanding the short-term ascending effects of spinal cord stimulation could answer several open questions, such as the mechanisms behind efficacy of TSS for sensorimotor disorders.
Current studies have focused on functional magnetic resonance imaging (fMRI) to understand the effects of spinal stimulation on the brain [33][34][35]. While fMRI provides high spatial resolution it relies on the blood-oxygen-level-dependent response and therefore has relatively low temporal resolution [36]. While several alternative methods for recording cortical activity exist, because of its safety, high scalp coverage, and high temporal resolution, electroencephalography (EEG) is one of the preferred non-invasive methods for recording brain activity [37,38].
In this study, we hypothesize that TSS causes measurable changes to cortical activation and functional connectivity via ascending connections of the spinal cord. To test this, we recorded EEG activity during different tasks while undergoing TSS or bilateral FES of the plantarflexors as a control condition. While TSS and FES utilize similar methods of non-invasive electrical stimulation, TSS is transmitted through the dorsal roots of the spinal cord and excites multiple motor pools at the same time along with the spinal circuitry [39,40], while FES excites individual muscles by activating both sensory and motor nerves [30]. We hypothesize cortical activation patterns will differ between the two modalities. We characterize cortical changes in the presence of stimulation using power spectral density in the timefrequency domain, the evoked responses in the timeamplitude domain, and the functional connectivity changes over time to provide a robust quantitative assessment of cortical changes during tonic contraction and alternating rhythmic contraction lower limb tasks in the presence of TSS or FES stimulation.

Participants
Five neurologically intact participants (height: 170 ± 18 cm, weight: 75 ± 13 kg, age: 28 ± 6 years, sex: 3 male) were recruited to participate in this study. Written informed consent to the experimental procedures, which were approved by the Houston Methodist Research Institute and University of Houston institutional review boards, was obtained from each participant. Figure 1 shows the experiment setup (A) and outlines the protocol (B) used in this study.

Experimental Setup
TSS was delivered to the skin over the lumbosacral spinal enlargement. FES was delivered bilaterally to the plantarflexors. Stimulation was administered with a constant current stimulator DS8R (Digitimer Ltd, UK) utilizing conductive self-adhesive electrodes (PALS, Axelgaard Manufacturing Co., Ltd, USA). For TSS, the circular cathode (size: 5 cm diameter) was placed at midline between the T12 and L1 spinous processes, and two oval anodes (size: 7.5 cm × 13 cm) were symmetrically placed on the abdomen approximately 5 cm lateral from the umbilicus [23,41]. For FES stimulation, the cathode and anode (size: 5 cm × 9 cm) were positioned over the plantarflexors bilaterally with the anode placed distal to the cathode.
Stimulation was delivered by a 500 µs biphasic rectangular pulse. Prior to the start of the experiment stimulation was delivered at a rate of 0.2 Hz to determine motor threshold and verify central placement of the electrode over the lumbar enlargement. TSS stimulation began at 30 mA and increased in increments of 5 mA until motor responses were observed in the lower limb muscles. FES stimulation began at 5 mA and was increased in increments of 5 mA until above motor threshold. During the experiment, stimulation frequency was set to either 15 or 30 Hz per the intervention protocol ( figure 1(B)). Stimulation amplitude for TSS was set at the participants' tolerance level (54 ± 9 mA) which ranged from 70% to 90% of motor threshold. FES stimulation was delivered at the same ratio to motor threshold (6 ± 2 mA). Both stimulation modalities did not induce muscle contractions in any of the participants.
Trigno Avanti wireless surface electromyography (EMG) electrodes (Delsys Inc., USA; common-mode rejection ratio <80 dB; size: 27 mm × 37 mm × 13 mm; input impedance: >1015 Ω/0.2 pF) were placed at placed at ten different sites: longitudinally over the left and right vastus lateralis, the semitendinosus portion of the medial hamstrings, tibialis anterior, gastrocnemius, and lateral aspect of soleus muscles. All EMG electrodes were placed over the muscle belly. EMG data were amplified with a Trigno Avanti amplifier (Delsys Inc.; gain: 909; bandwidth: 20-450 Hz) and recorded at a sampling frequency of 2000 Hz through a PowerLab data acquisition system (ADInstruments, Australia).
EEG data were recorded with 64-channel (Brain Products GmbH, Morrisville, NC) active Ag/AgCl electrodes amplified via a Brain Products BrainAmp DC amplifier (Brain Products GmbH) at a sampling frequency of 1000 Hz through BrainVision Recorder software (Brain Products GmbH). Electrodes were placed by following a modified 10-20 international standard arrangement. Channels AFz and FCz were used as the ground and reference, respectively. Channels TP9, TP10, PO9, and PO10 were removed from the cap and utilized for electrooculography (EOG) to measure eye-movements. The EOG channels were arranged such that TP9 and TP10 were placed on the left and right temples on the sides of the face, respectively and PO9 and PO10 were positioned superior and inferior to the right eye, respectively. Electrode impedance was maintained below 25 kΩ which was verified prior to the start of and after the completion of the experiment. To prevent the head from applying pressure on the sensors while the participant was in the supine position, a towel was used to support the head via the neck, which reduced the likelihood of motion artifacts.
Force output was measured with Exolab (Antex Lab LLC, Russia), a device incorporating four calibrated, two-sided precision load cells with a capacity of 50 kg (Futek, USA), designed to measure isometric forces generated by voluntary or stimulation-induced leg muscle activation while the participant's hip, knee, and ankle joints were supported at 155, 130, and 90 degrees, respectively. Measurements were recorded at a sampling frequency of 2000 Hz through a PowerLab data acquisition system (ADInstruments, Australia).

Experimental procedure
Participants were asked to perform one of three tasks: (1) no contraction (Control), (2) rhythmic isometric muscle contraction of the lower limbs mimicking walking (RC), or (3) tonic isometric contraction . EEG data denoising and analysis: Raw EEG data are cleaned by applying an H∞ filter which uses EOG recordings [33] to identify and remove ocular artifacts, signal drifts, and other shared sources of noise. The data are high pass filtered, line noise is filtered, then artifact subspace reconstruction is done to remove other types of artifacts. The data are common average referenced, then ICA is utilized to further identify and remove any remaining artifacts. The algorithm DIPFIT determines the equivalent dipole locations for each IC. Raw EMG is time corrected for the 60 ms delay introduced by the amplifiers, downsampled to 1000 Hz, then bandpass filtered, while force data is downsampled to 1000 Hz. Analysis of the cleaned data included calculating ERP, PSD, and gPDC.
to mimic standing (TC). During the RC task participants were instructed on the pace and amount of force to produce prior to the start the task. During the TC participants were instructed on the amount of force to produce at the start of the task. Each task was performed eight times for 30 s with one of three interventions: (1) no stimulus, (2) TSS, or (3) FES, with a minimum of 30 s between blocks. The no contraction and RC were performed at 30 Hz stimulation rates, which has been shown to activate locomotor specific spinal networks [15,39,42]. The TC task interventions were performed at 15 Hz, which has been shown to activate postural specific networks [21,43,44]. The experiment was performed three times over the course of one week with EEG recordings done during the first and last session to account for changes in cortical activity and functional connectivity due to task learning [45,46].

Signal pre-processing
EEG signal pre-processing and statistical analysis were performed off-line utilizing custom software written in MATLAB R2020b (MathWorks, MA) and functions from the open-access toolbox EEGLAB [47]. A flow-chart outlining the EEG signal pre-processing pipeline is shown in figure 2. EOG channels (4 channels) were utilized for adaptive filtering of the EEG signals (60 Channels) through the H∞ algorithm to remove eye artifacts, correct baseline drifts, and remove other shared sources of noise [48]. The data were then zero-phase high pass filtered by applying a fourth order ButterWorth filter with a cut off frequency of 0.1 Hz. An adaptive filter was then applied to remove harmonics of 60 Hz power line noise (CleanLine toolbox for EEGLAB) [47,49].
Artifact subspace reconstruction (ASR) was then applied to remove high amplitude artifacts (e.g. muscle bursts, sensor movement, etc). ASR utilizes a sliding window of EEG data with principal component analysis to identify channels of high variance by statistical comparison with a window of EEG data containing minimal artifacts. Due to volume conduction, a single channel capturing true cortical signals should not principally account for a large amount of variance within the given window. Channels which have a variance above a threshold when compared to the calibration data are identified as corrupted. In this study, two minutes of EEG recorded with the participant in the supine position without movement was taken as calibration data for ASR. Corrupted channels (or the subspaces of multiple channels) were reconstructed utilizing neighboring channels and a mixing matrix that is computed by the covariance matrix of the calibration data. A sliding window of 1 s and a variance threshold of five standard deviations were implemented to identify corrupted subspaces. The EEG channels were then re-referenced by subtracting to their common average.
Next, the equivalent current dipole that matched to the scalp projection of each independent component (IC) source was computed by a standard threeshell boundary element head model included in the DIPFIT toolbox [47]. Each IC scalp projection, its equivalent dipole's location, and its power spectra were then visually inspected and ICs that related to non-brain artifacts (e.g. sensor movement, muscle artifact, stimulation artifact) were removed. Force data were downsampled to 1000 Hz.

K-means clustering
For a group-level comparison of EEG data at the source level (ICs) K-means clustering was performed. Only the ICs which the equivalent dipoles explained greater than 80% of variance of the IC scalp projections were retained for clustering and the optimum number of clusters was determined by three different algorithms (Calinski-Harabasz, Silhouette, Davies-Bouldin) this was done to create a model that best fits the data [50][51][52]. Data were clustered by their equivalent dipole locations only to avoid circular inference, which is known to inflate the false positive rate [53,54]. Then the lowest number of clusters found across the three algorithms was selected. ICs were clustered by applying the k-means algorithm with k = 5 cluster centroids to the dipole locations and their respective distance between each other in vector space using squared euclidean distance. ICs with a dipole whose distance was more than three standard deviations from a cluster centroid were categorized into an outlier cluster and omitted from further analysis. The clustering algorithm was repeated 10 000 times to increase the chances of escaping a local minima and utilized randomly initialized centroid positions. This resulted in a total of 51 IC's retained across clusters for further analysis. Figure 3 provides a flowchart of the steps taken to determine generalized partial directed coherence (gPDC) from clustered IC data. By using the multivariate autoregressive representation (MVAR) of EEG signals the direction of the information flow between channels be inferred and direct or indirect influences detected. For a multivariate M-channel process X(n) defined:

Functional connectivity
the MVAR model can be represented as: where A k are M × M MVAR parameter matrices comprised by the coefficients a ij (k) that relate the ij series at lag m which describes the interactions between time series pairs over time, and is the vector of model innovations (zero mean with a covariance matrix Σ w ). The variable p is the order of . Calculating gPDC from clustered IC's: the mean of each cluster is calculated, the data are then segmented by condition. Bayesian information criterion is applied to determine optimal model order for each condition, then the mean model order is taken across conditions. A final model is created for each condition using a 2 s window with 50% overlap. The model coefficients are then used to calculate gPDC for each condition.
the MVAR model, which is determined with Bayesian information criterion (BIC) and is expressed by: where σ e is the variance of the prediction error for the model, σ X is the variance of the EEG signal, and N is the number of samples. Using BIC, the optimal order of an AR model is one that minimizes equation (4) [55,56]. In this case, BIC is computed iteratively from p = 4 to 100 for each of the conditions and the model order that obtained the lowest BIC is stored in a matrix. Next, the mean of the model order matrix is found and a final AR model with a 2 s window with 50% overlap is constructed for each condition using the mean model order value. Partial directed coherence (PDC) and gPDC are just a two MVAR-based measures which have been introduced to determine directional influence in multivariate systems [57,58]. The partial directed coherence measure was formally defined by Baccala and Sameshima in 2001 [57]. PDC can be defined as: where A(f) is the Fourier transform of MVAR model coefficients A k , A ij is the ijth element of A(f), and A * ij denotes the complex conjugate operation. Because the PDC output is normalized the output is bounded from the interval [0, 1]. However, severely unbalanced predictive modelling errors (innovations noise) can cause a reduction in performance. The gPDC measure was proposed to resolve cases where the innovation matrix is unbalanced by applying normalization factor in the denominator [58]. The gPDC equation can be expressed as: where σ k signifies the kth diagonal element of covariance matrix Σ w .

Statistical analysis
Permutation testing is a non-parametric test and requires no knowledge of how the test statistic of interest is distributed under the null hypothesis (e.g. no significant group difference). This is because the permutation process empirically 'generates' the distribution under the null hypothesis from the data by permuting data labels [59]. This type of statistical analysis is increasingly accepted with EEG and across the neuroscience disciplines [60][61][62]. Figure 4 outlines the steps in determining significant differences between the neuromodulation conditions (stim) to the corresponding no stimulus condition for a given task (control, TC, or RC). First the gPDC for each condition are segmented into frequency bands δ (0.1-4 Hz), θ (4-8 Hz), α (8-12 Hz), β (12-30 Hz), and γ (30-50 Hz). Next, the observed difference for each band is calculated where test statistic T obs is the observed value for a given band and a given task (control, TC, RC) and is the difference between the electrical stimulation condition minus the no stimulation condition.
The null hypothesis distribution T null is generated by randomly shuffling the stimulus and no stimulus condition labels, subtracting the two variables, then repeating this process for L = 10 000 times. Under the null hypothesis there is no difference between the two conditions, thus the labels are artificial. By shuffling the labels and recalculating the value we can determine if the observed values (T obs ) are significantly different from the test distribution (T null ).
Significance was determined by summing the instances where values of T null (k) are larger than or equal to T obs (k) which is then divided by the total number of permutations such that: where Θ is a step-function and, T null (k) and T obs (k) are a M x M matrices, X(i, j) is the ijth element of the difference between T null and T obs , Steps to determine significance: gPDC data are segmented into frequency bands, the difference of the stimulus condition during a given task (no contraction, RC, or TC) is subtracted from the corresponding condition with no electrical stimulus (e.g. TC with TSS minus TC with no stimulus). Separately, the stimulus and no stimulus data labels are shuffled 10 000 times and the difference between the shuffled classes is taken, then the mean and standard deviation are calculated and the difference data is then tested under the null hypothesis. A Bonferroni correction is applied to account for multiple comparisons (p adj ). and p(i, j, k) is the level of significance for the ijkth element.
Significance threshold was corrected for multiple comparisons applying a Bonferroni correction defined as: where k = 20 is the number of comparisons performed, α = 0.05 is the threshold value to determine if a test statistic is statistically significant and α critical was set at 0.00256.

Independent components of event related potentials
Each cluster's ICs of the ERP was determined by lowpass filtering the data with a zero phase fourth order ButterWorth filter with a cutoff frequency of 300 Hz to minimize distortion [63]. Data are divided by tasks where stimulation was applied and the data were segmented by the stimulus window where t = 0 is stimulus onset and window length was 30 and 60 ms for 30 and 15 Hz stimulus conditions respectively. The 30 Hz conditions use two window lengths to display consecutive responses. Latency was determined by the peak of the first response occurring later than either 10 ms for TSS or 25 ms for FES. These times were selected to account for conduction velocity, where the typical latency caused by median nerve stimulation

Power spectral density
The difference for each task was taken by subtracting the mean power spectral density (PSD) of the stimulus condition (e.g. TSS, FES) from the no stimulus condition from the same task. First, the forward projection of the cortical sources (clusters) to scalp sensors is done. Then PSD was calculated with a 1 s sliding window with 95% overlap applying Thomson's multitaper estimate in MATLAB, with 4096 frequency bins (1-50 Hz) and half-bandwidth product of 5. The Thomson multitaper method was chosen to help reduce the variance and spectral leakage associated with estimating the PSD with periodograms in a non-parametric manner [66]. The mean PSD was determined by averaging over blocks of repetitions and segmented into frequency bands. Figure 5 shows the results of the k-means clustering algorithm, which resulted in five clusters and 51 total dipoles included in the analysis. Figure 5 also shows the number of participants that make up each cluster, the number of IC sources, the Talairach coordinates, and the Brodmann area (BA) for the corresponding cluster. Talairach coordinates and BA were identified from the Talairach atlas [67]. The BA were searched within ±5 mm cube ranges around each cluster centroid.

Results
Changes in functional connectivity over time are seen in figure 6 for the rhythmic contraction (RC) condition for time t = 1, 10, 20, and 30 s broken into band power. Only significant differences (p adj < 0.00256) are shown comparing the RC with TSS condition to the RC without TSS condition. We found differences greater than 30% in functional connectivity during TSS over time across bands. Overall there was a trend of suppression between motor areas of the brain (clusters 2 and 3) and an overall increase in connectivity in clusters 1, 4, and 5. TSS during the rhythmic movement condition caused changes in functional connectivity up to 33.4% difference to the rhythmic movement without stimulus condition. Tonic contraction in the presence of TSS had differences up to 27.4% when compared with the tonic contraction and no stimulus condition. No significant changes in functional connectivity were found for movement conditions where FES was utilized. Videos demonstrating functional connectivity over time for differences between the conditions with TSS intervention and the control conditions can be found in the supplementary materials.
The ICs of the evoked responses caused by stimulation are shown in figure 7 and are organized by condition (column) and cluster (row). The dark purple or dark pink line is the mean of either 3592 responses with two consecutive stimuli shown for 30 Hz stimulation conditions or 3592 responses for 15 Hz stimulation conditions and the 95% confidence interval is demarcated by the light purple or light pink. We found that 30 Hz TSS caused a clear response peaking at approximately t = 20 ms after stimulus onset for both the control and RC conditions across clusters, but appears to be instantaneous, while 30 Hz FES caused no significant response in either condition. During 15 Hz TSS the response during TSS is between t = 20-25 ms while 15 Hz FES caused a minor response in clusters 2 and 3 at approximately 50 ms after stimulus was applied. Figure 8 shows the spatial distribution differences in band power for each cluster between the RC with stimulus and the RC without stimulus conditions by frequency band. In most cases only modest changes in band power were seen with maximum changes across clusters and bands of ±3 dB. For Cluster 1 delta band had the largest differences with a absolute maximum difference of 2.64 dB in the δ band. Clusters 2 and 3 were found to have smaller changes with cluster 2 in the α band with an absolute maximum difference of 2.53 dB while for cluster 3 the γ band was even lower at 2.18 dB. The largest changes were found in cluster 4 for the γ band which was found to have an absolute maximum of 2.97 dB. Cluster 5 had the largest changes in the α band which had an absolute maximum of 2.77 dB. Two bands, θ and β, had minimal changes across clusters while δ, α, and γ were found to have at least one cluster with large changes.

Discussion
To better understand the effects of TSS on cortical activity, we investigated several different tasks performed in the presence of either TSS or FES. We clustered dipoles across subjects into five common Due to the high frequency, the pink response (40-60 ms) is most likely due to the previous pink stimulus (0-5 ms) as conduction velocity is not fast enough to cause responses with close to zero latency. This is demonstrated utilizing two different colors to differentiate the two stimulation periods (current and prior stimulation) within the analysis window during the 30 Hz TSS conditions. Figure 8. Scalp maps of power spectral density contribution difference by cluster: the difference in PSD contribution of each cluster (row) for the RC condition by frequency band (columns). Difference in PSD was determined by forward projecting the cortical source (cluster) to the scalp sensors, then subtracting the PSD of the RC with stimulus condition from the PSD for the RC task with no stimulus condition. areas of activation. We then demonstrated distinct changes in functional connectivity between these clusters during TSS when compared to the no stimulus conditions, further these changes were not found during FES. We also observed increases and decreases in functional connectivity across all frequency bands in the presence of stimulation, highlighting their contributions in both rhythmic and tonic tasks. These results help better characterize the effects of TSS on cortical activity and establish a baseline for comparison with clinical populations.
Several of the clusters found are associated with movement planning and execution (figure 5). For example, BA 9 (Cluster 1) of the right hemisphere has been shown to be involved in planning and working memory [68][69][70]. BA 4 (Cluster 2) is the designation for the primary motor cortex which is associated with executing movements and sends direct projections to the spinal cord [71,72]. BA 6 (Cluster 3) is composed of the premotor cortex and, medially, the supplementary motor area which is associated with sensory guidance of movement, control of proximal and trunk muscles. This area also contributes to the planning of complex and coordinated motor movements [73][74][75].
BA 18 (Cluster 4), also referred to as the right visual association cortex (V2), is the region of the brain associated with visual selective attention tasks when stimuli are present [76][77][78]. Finding IC's that make up this cluster across all participants was unexpected because none of the tasks involved visual cues or required visual memorization.
One potential explanation for the recruitment of the V2 cluster, proposed by authors JLC-V and AGS, is that spinal stimulation activates neurons in the spinomesencephalic tract, which synapses at the periaqueductal gray (PAG). The PAG is located around the cerebral aqueduct within the tegmentum of the midbrain and is associated with the release of endogenous opioid neurotransmitters [79]. Furthermore, connectivity analysis suggests the secondary visual cortex preferentially innervate the dorsolateral PAG [80][81][82]. Direct stimulation of the PAG has been shown to cause analgesia and while further study is needed, stimulation induced activation of the PAG would explain the analgesic effects during ESS and could answer why those effects diminish over time [18,83]. Future studies following this experimental design while utilizing simultaneous EEG and fMRI may help determine if this is the case.
Alternatively, authors GAM and DGS propose that the changes in visual/occipital cortex, as found in the V2 cluster, can be attributed to the experimental setup, where the larger pressure on the electrodes at the back of the head resulted in their increased contact with the scalp. This could affect the amplitude of stimulation artifacts, which in turn could skew the magnitude of activities observed in the V2 cluster. Note that authors JLC-V and AGS do not agree with this explanation because if this were the case, the effects would be generalized across the frequency bands in the functional connectivity analysis, would be seen across the scalp, and would not be specific to right visual cortex. Performing this experiment in a seated position with the head selfsupported would determine if these findings were caused due to increased contact with the scalp. Future work should also include eye tracking, which would eliminate changes in involuntary eye movement as a potential cause of V2 activation.
Cluster 5 was found at BA 30 which is known as the retrosplenial cortex [84]. The retrosplenial cortex has projections to several different areas of the brain including: the hippocampus, the prefrontal cortex, the parietal cortex, and the visual cortex [84,85] Functional imaging has suggested that it is involved in memory, visuospatial processing, proprioception, and emotion [86,87]. This area has also been associated with the perception of pain [82,[88][89][90].
We found several temporally distinct changes in functional connectivity within tasks using cortical sources and permutation testing. The largest differences between TSS and no stimulation were greater than ±30% as seen in figure 6. Our findings suggest that TSS induces higher functional connectivity, the effects are immediate, do not require a 'build up' period, and do not increase over time. More importantly, the effects were seen across bands for both TC and RC conditions, which indicates that a wide range of frequencies, from delta band to gamma band, are involved in rhythmic and tonic tasks, a finding that matches similar studies done during walking without TSS [91]. Our results also suggest that TSS tends to increase connectivity between areas of motor planning (Clusters 1 and 5), but reduces areas associated with direct motor control (Clusters 2 and 3).
A similar functional connectivity analysis was performed using the forward projection of the cortical sources to EEG sensor space, but due to the high number of comparisons the statistical power was so low that no significance was found. This sensor based analysis result matches the findings from Goudman et al which found only a single sensor pair that had significant increase in functional connectivity for the FC3-TP9 electrodes in the beta frequency band while the participant underwent high dose ESS [92]. For this reason we believe that performing functional connectivity analysis using cortical sources provides a more accurate representation of information flow. One limitation with this analysis is the window size required to estimate low frequency changes. This creates a trade-off between low frequency estimation and a high time resolution.
No changes in functional connectivity were found during FES, which suggests that sub-motor threshold FES does not cause significant alteration to cortical activity and by extension spinal activity. This finding may explain why TSS can enable voluntary multi-joint lower limb movements in people with SCI [15,93], while FES requires a complex timed stimulation approach to achieve the same result [44,94].
The effect of the stimulus artifact on functional connectivity was an initial concern, but the 17 ms window utilized in the calculation of the gPDC was much shorter than the 33 or 67 ms inter-stimulus intervals used during 30 or 15 Hz stimulus respectively. We also found several decreases in the functional connectivity during stimulation when compared to the no stimulus condition, even in the β and γ bands where the highest effects from the stimulus artifact should be present ( figure 6). This further supports the theory of limited, or even no, meaningful contribution from the artifact on functional connectivity analysis because we would expect increases exclusively if there was a contribution.
The latency of the ERP in IC space (figure 7), measured at first peak across 15 Hz TSS stimulus conditions was found to be 21 ± 3 ms, which matches well with the latency found for ESS stimulation [95]. A limitation to our analysis is that ERP studies are typically done with a stimulation rate of between 1 and 5 Hz [63,64]. The ERP analysis here was done at a frequency of either 30 or 15 Hz and higher frequencies of stimulation have been shown to reduce the amplitude of response [63]. While 0.2 Hz was used initially to determine motor threshold only 10-20 trials were performed per subject. To ensure reproducibility it is recommenced to use at least 500 trials with two separate blocks of trials performed for each subject [64]. During the 30 Hz TSS stimulation conditions, while the IC responses peaks at around 20 ms, it appears that several, marked with two separate colors, have responses starting immediately after stimulus onset. This result could be due to the continuation of the response from the previous stimulus or the ERP response could be caused by the previous stimuli. To demonstrate this the first stimulus is colored in pink with the second response in the same color. Further supporting this theory, when 15 Hz TSS was performed, the responses do not return to baseline before the next stimulus was delivered. We therefore believe a response overlap effect was occurring. The purpose of this study was exploratory and thus did not use the parameters specifically for ERP analysis, which makes interpretation of the results difficult. Further study should be done to investigate the evoked responses from TSS specifically. Analysis of the evoked response could aid in the understanding of the pathways being activated during TSS.
While the pulse width for both TSS and FES was very short, lasting just 500 µs, the artifact is clearly present in several of the cluster ERP's. The artifact was minimized by applying a combination of ICA, H∞, and cluster based analysis, but was not fully eliminated in all clusters/conditions. This demonstrates the difficulties with removing electrical artifacts from neuromodulation interventions like TSS. Further work will need to be done to determine the best ways to minimize or, if possible, eliminate neuromodulation artifacts.
One major limitation of this study is the sample size (n = 5). While participants had EEG recordings taken twice over the course of the experiment, future work should expand the number of participants and explore the effect of other stimulation frequencies and motor tasks on cortical activity. The main purpose of this study was to establish a baseline and to investigate the effects of TSS on cortical networks so that the data can be translated and investigated in different clinical populations. The proposed protocol should be performed in individuals with various degrees of spinal cord injury to assess the effects TSS on cortical activity in this clinical population. Comparison of clinical populations to the cortical responses from the intact spinal cord could help determine adaptive changes after SCI and provide a direct way to quantify severity of injury.
To maximize comfort for the participants of this study, the tasks were performed in the supine position as shown in figure 1(A). This position reduces discomfort by lowering the impedance between the cathode and the surface of the skin and causes the conus medullaris and nerve roots to shift dorsally reducing the response threshold [96,97]. This position limits the movement tasks a participant could reasonably perform and was a potential source of noise due to possible pressure on the EEG electrodes or EMG activation of the neck. Future studies could utilize mobile EEG, which would allow for the investigation of cortical changes during complex walking tasks (e.g. stair ascent and descent) in the presence of TSS. Further study should also be done using a longitudinal methodology to determine the long-term changes in functional connectivity.

Conclusion
To the authors knowledge this is the first study investigating the effects of lumbosacral TSS on cortical activity using EEG. We clustered equivalent dipoles across participants (N = 5) and identified five clusters for analysis. We found significant differences in functional connectivity between clusters during TSS when compared to baseline, but not during FES. Both TSS and FES caused evoked potentials in cortical sources, with TSS mean latencies of 21 ± 3 ms and FES mean latencies of 42 ± 4 ms, respectively. Our findings demonstrate that TSS causes complex and broad changes in cortical networks, changes which were not seen during FES. The findings of this study could aid in determining optimal parameters for treatment of people living with spinal cord injury, assist in the design of a closed-loop system for spinal stimulation, and the design of spinal cord based computer interfaces.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors. ted by the IUCRC BRAIN (National Science Foundation Award # 1650536) to J L C-V, and philanthropic funding from Paula and Rusty Walter and the Walter Oil & Gas Corporation as well as by TIRR Mission Connect Grant No. #019-102-Jerry Johnston Andrew Award to D G S. The funders were not involved in the design of the study, the collection, analysis, and interpretation of the experimental data, the writing of this article, or the decision to submit this article for publication.
The authors have confirmed that any identifiable participants in this study have given their consent for publication

CRediT authorship contribution statement
Alexander G Steele: data acquisition, methodology, software programming, formal analysis, data curation, writing-original draft. Gerome A Manson: methodology, data acquisition, writing-review & editing. Philip J Horner: conceptualization, supervision, writing-review & editing. Dimitry G Sayenko: conceptualization, funding acquisition, methodology, resources, supervision, visualization, and writing-review & editing. Jose L Contreras-Vidal: conceptualization, funding acquisition, methodology, resources, supervision, visualization, EEG analysis strategy, and writing-review & editing. All authors have reviewed and approved the manuscript.