Electrospinography for non-invasively recording spinal sensorimotor networks in humans

Ten neurologically intact participants (Height: 165 ± 15 cm, Weight: 72 ± 8 kg, Age: 30 ± 12 years, Sex: 4 women) were recruited to participate in this study. Written informed consent to the experimental procedures, which were approved by the University of Houston institutional review board, was obtained from each participant. This research was conducted in accordance with the principles embodied in the Declaration of Helsinki and in accordance with local statutory requirements [17]. Figure 1 provides an overview of the experiment setup ((A) and (B)) and outlines the experimental protocol (C) used in this study.

Zoom In Zoom Out Reset image size

EEG and ESG data were recorded with (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).

Thirty-two channels of EEG were placed following a modified 10–20 international standard arrangement. Channels TP9, TP10, PO9, and PO10 were removed from the cap and utilized for electrooculography (EOG) to measure eye movements and eye blinks. The EOG channels were arranged so the sensors TP9 and TP10 were placed on the left and right temples, respectively, and PO9 and PO10 were positioned superior and inferior to the right eye, respectively. Ground and reference were placed on the left and right earlobes, respectively.

Nineteen channels of ESG were placed between the intervertebral spaces, starting at T10/T11 through L2/L3, found via palpation. Where three sensors were placed at each intervertebral level, with the middle sensor placed along the midline of the spinal column and the left and right sensors placed directly next to the midline sensor. Four additional sensors were placed bilaterally, approximately 5 cm from the midline, just above and below the array. These sensors were utilized to measure EMG, cardiac (EKG), and other artifacts. The electrode impedance was maintained below 25 kΩ which was verified prior to the start of and after the completion of the experiment. To prevent excessive flexion of the neck while the participant was in the prone position, a face-down foam support was used.

Trigno Avanti wireless surface EMG electrodes (Delsys Inc. USA; common-mode rejection ratio $\lt$80 dB; size: 27 mm × 37 mm × 13 mm; input impedance: $\gt$1015 Ω/0.2 pF) were placed at eight different sites: bilaterally over the rectus femoris (RF), the medial hamstring (MH), the tibialis anterior (TA), and the lateral aspect of the soleus (SOL) muscles. EMG data were amplified by a Trigno Avanti amplifier (Delsys Inc.; gain: 909; bandwidth: 20 to 450 Hz). Each sensor contains an integrated 9-degree of freedom inertial measurement unit (IMU) (3-axis accelerometer; sampling resolution: 0.016 ± 0.001 g/bit; noise: $\lt$0.01 g RMS; resolution: 16 bit; 3-axis gyroscope; range: ± 2000^∘, and 3-axis magnetometer; range: ±4900 uT). EMG and IMU data were sampled at 1925 Hz and 148 Hz, respectively. In addition to the adhesive backing used to attach the sensor to the skin, an elastic underwrap athletic foam tape was used to hold the sensors in place to minimize movement artifacts.

Participants were asked to perform one of three tasks: (1) no movement (control), (2) plantarflexion of the left (LPF) or right (RPF) foot, and (3) knee flexion of the left (LKF) or right (RKF) leg, with all participants performing both left and right movements. Prior to the start of each task, participants were given a demonstration of the movement they were to perform, instructed to remain relaxed between repetitions, and to move only the targeted joint to the best of their ability. Movement tasks were performed in five blocks of five repetitions for a total of 25 receptions per movement. Movements were self-paced, and the participant was told to start the block 1–2 s after a verbal cue signifying the start of the block. The time between conditions was no less than 2 min, and blocks within a given condition occurred no less than 1 min apart. The control task (no movement) was performed for 2 min each time preceding and following the two movement conditions to account for possible changes in resting-state cortical and spinal activity over the course of the experiment [25].

Signal pre-processing and statistical analysis were performed off-line utilizing custom software written in MATLAB R2022b (MathWorks, MA) and included functions from the open-access toolbox EEGLAB [19]. A flowchart outlining the signal pre-processing pipeline and analysis is shown in figure 2 and follows a previously utilized pre-processing methodology [24].

Zoom In Zoom Out Reset image size

EOG data (4 channels) were used to filter the EEG signals (28 channels) via $H\infty$ an adaptive algorithm that removes artifacts such as eye blinks or eye movements (q = 10⁻⁹, γ = 1.15) [18]. The data were then zero-phase band-pass (0.1 Hz–100 Hz) filtered by applying a fourth-order ButterWorth filter.

Next, artifact subspace reconstruction (ASR) was performed to remove high-amplitude artifacts (e.g. muscle activity, sensor movement, wire tugs, etc). Just prior to the start of the experiment, 2 min of EEG recorded with the participant in the prone position without movement were taken as calibration data for ASR and were not utilized as the control condition during the experiment. Corrupted subspaces were reconstructed utilizing the neighboring channels and a mixing matrix that is computed by the covariance matrix via the calibration data. For this process, a sliding window of one second and a variance threshold of four standard deviations were implemented to identify and rebuild the corrupted subspaces.

Two complementary adaptive filters were applied to remove harmonics of 60 Hz power line noise (CleanLine and ZapLine Plus toolboxes for EEGLAB) [19, 26–28]. The EEG channels were then re-referenced by subtracting their common average. Next, using the standard three-shell boundary element head model included in the DIPFIT toolbox, the equivalent dipole that matched the scalp projection of each independent component (IC) was computed [19]. Each IC scalp projection, the equivalent dipole's location, and the power spectra were then visually inspected, and ICs that related to non-brain artifacts (e.g. sensor movement, muscle artifact, etc) were removed from the analysis.

To establish a common group-level comparison of EEG data at the cortical source level, K-means clustering was performed. Specifically, equivalent dipoles of ICs that explained greater than 80% of the data variance were retained for dipole clustering. The number of clusters was determined by three different algorithms (Calinski–Harabasz, Silhouette, Davies–Bouldin) so that the best fit could be achieved [29–31]. This process concluded that the best fit was k = 5, 5, and 9 for Calinski–Harabasz, Silhouette, and Davies–Bouldin, respectively. Because two of the algorithms returned the same value, the ICs were clustered by applying the k-means algorithm with k = 5 centroids and dipoles were clustered using the squared euclidean distance with respect to equivalent dipole locations only, which was done to avoid circular inference [32]. If an IC with a dipole whose distance from a cluster centroid was more than three standard deviations, it was omitted from further analysis. This resulted in a total of 107 ICs being retained across participants and clusters for further analysis.

Four channels, placed bilaterally 5 cm lateral to the midline and approximately parallel to the top and bottom of the lumbar ESG array, were utilized to filter EMG, EKG, and other sources of noise. These sensors are referred to in the text as 'EKG sensors.' The EKG sensors were utilized by the $H\infty$ adaptive filtering algorithm, similar to what was used prior for eye artifacts in the EEG recordings.

In this case, two values q and γ were adjusted for this application (q = 10⁻⁷, γ = 1.2). Specifically, the $H\infty$ adaptive weight estimation filter with time varying weights can formulated such that

$$\begin{equation} \hat{w}_{i+1} = \hat{w}_i + \frac{P_ir_i}{1+r_i^TP_ir_i}y_i \end{equation} \tag{ 1 }$$

where $\hat{w}_{i+1}$ is the estimated weight per channel at sample i + 1, r_i is a vector of reference measurements, in this case the EKG sensors, y_i is defined such that

$$\begin{equation} y_i = s_i - r_i^T\hat{w}_i \end{equation} \tag{ 2 }$$

where

$$\begin{equation} P_i^{-1} = \tilde{P}_i^{-1} - \gamma^{-2}r_ir_i^T , \end{equation} \tag{ 3 }$$

and $\tilde{P}_i^{-1}$ is the noise covariance matrix which satisfies the equation

$$\begin{equation} \tilde{P}_{i+1} = \left[\tilde{P}_i^{-1}+\left(1-\gamma^-2\right)r_ir_i^T\right]^{-1} + qI \end{equation} \tag{ 4 }$$

which is initialized with $\tilde{P}_0 = \mu I$ where µ is a constant. gamma is the bound on the energy-to-energy gain from the disturbances to the output estimation error and defines the levels of disturbance tolerated by the $H\infty$ filter.

In this formulation the weights are assumed to be time-varying with unknown variation dynamics and when γ < 1 the filter provides a guaranteed robust performance under all levels of disturbances in the system and the formulation is valid for the inequality $\gamma^2 \unicode{x2A7D} 1+ q\bar{r}$ where $\bar{r}$ is the supremum of the measured disturbances and q reflects the a priori information of how rapidly the weights vary over time. Thus, by adjusting both the selected value of γ and q, we can modify the behavior of the $H\infty$ filter for both applications [18].

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. Then ASR was performed using a sliding window of 1 s and a variance threshold of four standard deviations. To remove harmonics of 60 Hz power line noise the approach previously outlined for EEG data cleaning was applied (CleanLine and ZapLine Plus toolboxes for EEGLAB) [19, 26–28]. Lastly, ESG channels were re-referenced by subtracting their common average.

The data were segmented by condition, then by block, where each condition had five blocks consisting of five repetitions. First, PSD was calculated with a 750 ms sliding window using 99% overlap to apply Thomson's multitaper estimate in MATLAB, with 4990 frequency bins (0.1–50 Hz) and a half-bandwidth product of four for the entirety of each block [24, 33].

Next, the PSD data were further segmented by start and stop periods, where a movement that caused EMG changes larger than two standard deviations from the mean was defined as the start of the movement, and the end of the movement was defined as the point where EMG returned below this cutoff. PSD data was averaged over windows that fell completely within the movement period. To compare to the baseline, the same number of windows were taken between 5 and 20 s prior to the first movement within each block. The data were averaged across repetitions and blocks within individual participants. ESG sensors were analyzed using a bipolar configuration to reduce common noise and increase differences between ESG sensor pairs.

Next, to determine the effect noise had on the resultant data, the PSD was also calculated using the four EKG sensors at the same time points and analyzed by sensor.

Figure 2 (data analysis) provides a flowchart of the steps taken to determine generalized partial directed coherence (gPDC) from within participant clustered IC data, where IC's belonging to the participant were selected, and between participant IC's and ESG data, which was modified from a previously utilized methodology for determining functional connectivity [24]. Due to the differences in signal characteristics between the cortical sources and ESG recordings, we expected the model orders that best described the clustered data only and the model orders that best described the entire data (e.g. connectivity between clustered ICs and ESG) would be different. For that reason, two models were utilized to best represent the two pathways.

Using this method, two models were created for each participant, utilizing the ICs in the clusters belonging to that participant only and the participants' ICs and their respective ESG data, which was vertically concatenated with the IC data. Each model was computed from a model order p = 4 to 100 for each condition and each participant, then stored in a matrix where the first row is the minimum BIC value calculated between the ICs for that participant and the second row is the minimum BIC determined for the participant ICs and their respective ESG data. Therefore, each column represents the minimum BIC for that condition for the two models, and each model order matrix represents one participant.

Next, the mean of each row in the model order matrix is calculated, and a final multivariate autoregressive (MVAR) model with a two-second window and 50% overlap is constructed for each participant and condition using the determined model order [22, 34]. The MVAR model can be described by the equation

$$\begin{equation} \textbf{X}\left(n\right) = \sum_{k = 1}^p\textbf{A}_k X\left(n-m\right)+\textbf{E}\left(n\right) \end{equation} \tag{ 5 }$$

in which X(n) is a multivariate M-channel process of length n, A_k are MxM parameter matrices that are comprised of the coefficients $a_{ij}(k)$ that relate the ijth series at lag m and describe the interactions between time series pairs over time. The vector of model innovations E has a zero mean and a covariance matrix $\boldsymbol{\Sigma}_w$, and the model order p is the corresponding model value.

To determine information flow, gPDC was applied because to the MVAR models because it can be utilized to determine directional influences and can detect the interaction of multiple influences (M > 2) [23, 35]. The gPDC equation is given as

$$\begin{equation} \mathrm{gPDC}_{j \rightarrow i}\left(f\right) = \frac{\frac{1}{\sigma_i}\bar{A}_{ij}\left(f\right)}{\sqrt{\sum\limits_{k = 1}^N\frac{1}{\sigma^2_k}\bar{A}_{kj}\left(f\right)\bar{A}_{kj}^*\left(f\right)}} \end{equation} \tag{ 6 }$$

where $\textbf{A}(f)$ is the Fourier transform of MVAR model coefficients A_k , $\bar{A}_{ij}$ is the i jth element of $\textbf{A}(f)$, $\bar{A}^*_{ij}$ signifies the complex conjugate operation. The difference between gPDC and PDC is a normalization factor used to correct for cases where the innovation matrix is unbalanced which occurs, in part due to noise in the model, the normalization factor is σ_k and signifies the kth diagonal element of the covariance matrix $\boldsymbol{\Sigma}_w$. The gPDC output is normalized and therefore bound by the interval [0,1].

Functional connectivity was used to determine information flow within the participants' IC's and between the participants' IC's and ESG data to minimize or eliminate the effect of any remaining muscle artifact due to preparatory or postural back muscle contraction before and during movement, which may cause a false positive when correlating ESG to EMG activation.

The 95% confidence interval of the PSD data were determined via the Signal Processing Toolbox available for MATLAB (R2022b) through the pmtm function which uses the Chi-squared approach.

Permutation testing was used to determine gPDC significance because it requires no knowledge of how the test statistic of interest is distributed. The permutation process empirically 'generates' the distribution under the null hypothesis to determine if a difference exists [36, 37].

To perform permutation testing, the gPDC is first separated into the frequency bands δ (0.1–4 Hz), θ (4–8 Hz), α (8–12 Hz), β (12–30 Hz), and γ (30–50 Hz) by condition. Next, the test statistic $\textbf{T}_\mathrm{obs}$, which is the difference in the observed value for the band and task (i.e. control, RPF, LPF, RKF, LKF) being permuted is calculated as the movement condition minus the control (no movement) condition.

The null hypothesis distribution $\textbf{T}_\mathrm{null}$ is generated separately by randomly shuffling the task (movement) and control (no movement) labels, then subtracting the two variables, and repeating this process for L = 10 000 times. If the null hypothesis is correct, then there are no differences between the two conditions, therefore by shuffling the labels we can determine if the observed values ($\textbf{T}_\mathrm{obs}$) and the null distribution ($\textbf{T}_\mathrm{null}$) are different.

Significance was determined by summing the instances when values of $\textbf{T}_\mathrm{null}(k)$ are larger than or equal to $\textbf{T}_\mathrm{obs}(k)$ and is divided by the total number of permutations [24].

As multiple comparisons were made, a Bonferroni correction was applied and the corrected α, $\alpha_\mathrm{critical}$ was set at 0.0025. For the gPDC between clusters and ESG sensors, there were k = 40 comparisons performed and $\alpha_\mathrm{critical}$ was set at 0.00 125.

Figure 3 demonstrates the results from the $H_\infty$ filtering algorithm for cardiac artifact removal and shows ESG data prior to the $H_\infty$ filtering step (black) and ESG data after the $H_\infty$ filter has been applied (blue) for a (top left) 50 s time window, a (top right) two second time window, and a (bottom) 200 ms time window at decreasing amplitude scales, ± 100, 50, 30 µV, respectively. The data from the top right and bottom middle plots are denoted by the black square.

Zoom In Zoom Out Reset image size

Figure 4 shows the results of the k-means clustering algorithm, which resulted in five clusters and 107 total dipoles. Figure 4 also shows the number of participants that make up each cluster, the number of IC sources per cluster, the Talairach coordinates, and the Brodmann area (BA) for the corresponding cluster. Talairach coordinates and BAs were identified from the Talairach atlas [38]. The BAs were searched within ± 5 mm cube ranges around each cluster centroid and the maximum distance from the centroid to the BA area was 3 mm.

Zoom In Zoom Out Reset image size

Figure 5(A) demonstrates PSD changes during left and right knee flexion when compared to rest employing a bipolar configuration denoted by the line connecting sensors for a representative participant (NIS007). The data demonstrate a significant (p < 0.05) decrease in PSD during left and right knee flexion when compared to resting state. Changes are localized to the top half of the array (top left and right PSD plots) while the bottom half of the array had minimal changes. It should be noted that while several participants, including the one shown, had decreases in PSD during movement while compared to rest, others (n = 3) had increased PSD compared to rest. Further, in some cases the level of PSD change found (increase or decrease) was side specific (i.e. left vs. right). Additional examples of changes in PSD during knee flexion and plantarflexion can be found in the supplementary materials.

Zoom In Zoom Out Reset image size

Analysis of the EKG sensors as seen in figure 5(B) for the top and bottom left and right of the array during the same time periods used for the analysis found minimal change in PSD during left or right knee flexion when compared to rest. Similar changes were found across participants with changes in the PSD during movement more prevalent during knee flexion than plantarflexion, however none of the changes found in the EKG sensors matched the pattern of changes seen in the ESG array.

Figure 6 demonstrates the functional connectivity changes within the participant IC data and between the cluster data and ESG data using a two second window of time just prior to right plantarflexion execution figure 6(A) and a two second window of time during movement preparation and execution (B). Normalized amplitude data is shown for the clusters (bottom left) and ESG data (bottom middle), while the time window is shown using a gray box overlaid on the EMG data for the right soleus (bottom right). Data are for a representative participant and videos presenting gPDC over time can be found in the supplementary materials.

Zoom In Zoom Out Reset image size