System and method for non-invasively recording spinal sensorimotor networks

Abstract

A method is described herein comprising receiving an electroencephalography (EEG) signal from a plurality of EEG sensors during performance of a plurality of tasks, receiving an electrospinography (ESG) signal from a plurality of ESG sensors during the performance of a plurality of tasks, receiving an electromyography data (EMG) signal from a plurality of EMG sensors during the performance of the plurality of tasks, filtering the EEG signal, the ESG signal, and the EMG signal, identifying independent components in the EEG signal for each participant using independent component analysis, wherein the independent components correspond to dipole activity associated with participants' performance of the plurality of tasks, and analyzing the independent components and corresponding ESG signal to determine associations among participant independent components and between each participants independent components and respective ESG signal.

Description

[0001] UNITED STATES PATENT APPLICATION FOR:

[0002] SYSTEM AND METHOD FOR NON-INVASIVELY RECORDING SPINAL

[0003] SENSORIMOTOR NETWORKS

[0004] RELATED APPLICATIONS

[0005] This application claims the benefit of United States Application No. 63 / 611,068, filed December 15, 2023.

[0006] This application claims the benefit of United States Application No. 63 / 614,469, filed December 22, 2023.

[0007] TECHNICAL FIELD

[0008] The disclosure herein involves a system and method for noninvasively recording biopotentials over the surface of the spinal cord and by applying an adaptive filtering algorithm to remove the cardiac artifact, the underlying data can be used to query the state of the spinal cord. This can be used for diagnostic and prognostic applications after injury / disease, and for a spinal cord computer interface (SCCI) for assistive or therapeutical applications.

[0009] INCORPORATION BY REFERENCE

[0010] Each patent, patent application, and / or publication mentioned in this specification is herein incorporated by reference in its entirety to the same extent as if each individual patent, patent application, and / or publication was specifically and individually indicated to be incorporated by reference.

[0011] BACKGROUND

[0012] With the technology and signal processing advancements in recent years, mobile brainbody imaging (MoBI) systems are transformed into much more ambulatory devices, which opens up a potential for the advancement of the ecological validity of brain imaging research and more practical solutions for in-home medical monitoring and brain-computer interface (BCI) applications, as well in consumer electronics applications. With an increasing interest and demand in applying EEG scans in real-world environments, MoBI systems are developed to record brain dynamics during different tasks in the medical and non-medical fields. Even though there is an uptrend of developing commercial headsets in BCI-related research, consumer-like user-friendly headsets are still rare. Most of the commercially available portable systems are relatively expensive, require proprietary software to function, and lack flexibility or modularity. Ergonomically, headsets are not designed to be truly easy and intuitive to use. They often require trained technicians to help to put on the headset and operate the system. There is a growing need for a low-cost MoBI headset that can be set up with user-friendly ergonomics that is easy to operate and performs a quality scan and data collection consistently. Studies show most headsets on the market often do not fit as well as soft EEG caps. Headsets with a poor fit to the user will likely lose scanning signals due to the unstable sensor-skin contact and shifting position while in use. To date, traditional EEG caps are still the best in terms of accommodating both size and shape variation. There is a need for an easy to use one-hand operated headset that provides a custom fit for all users.

[0013] INCORPORATION BY REFERENCE

[0014] Each patent, patent application, and / or publication mentioned in this specification is herein incorporated by reference in its entirety to the same extent as if each individual patent, patent application, and / or publication was specifically and individually indicated to be incorporated by reference.

[0015] BRIEF DESCRIPTION OF THE DRAWINGS

[0016] Figure 1 shows filtering of cardiac artifacts, under an embodiment. Figure 2 shows an example of processed data, under an embodiment. Figure 3 shows an example of processed data, under an embodiment. Figure 4 shows an example of processed data, under an embodiment. Figure 5 shows an example of processed data, under an embodiment. Figure 6 shows an example of processed data, under an embodiment. Figure 7 shows normalized ESG data, under an embodiment. Figure 8 shows power spectral density (PSD) changes during knee flexion, under an embodiment.

[0017] Figure 9 shows power spectral density (PSD) changes during knee flexion, under an embodiment.

[0018] Figure 10 show power spectral density (PSD) changes during plantarflexion, under an embodiment.

[0019] Figure 11 shows power spectral density (PSD) changes during plantarflexion, under an embodiment.

[0020] Figure 12 shows power spectral density (PSD) changes during knee flexion, under an embodiment.

[0021] Figure 13 shows an experimental protocol, under an embodiment.

[0022] Figure 14 shows electroencephalograph (EEG) and electrospinography (ESG) data denoising and analysis, under an embodiment.

[0023] Figure 15 shows filtering of cardiac artifacts, under an embodiment.

[0024] Figure 16 shows Clusters of equivalent dipole sources fit to independent components, under an embodiment.

[0025] Figure 17 shows power spectral density (PSD) changes during knee flexion, under an embodiment.

[0026] Figure 18 shows functional connectivity changes over time, under an embodiment.

[0027] DETAILED DESCRIPTION

[0028] Example 1

[0029] 1. Electrospinography Data

[0030] This section contains example electrospinography (ESG) data from randomly selected participants and demonstrates data prior to and after processing in the time domain. It should be noted that the plotted data showing raw (e.g., unprocessed) data were high-pass filtered at 0. 1Hz to remove baseline drift and improve readability, but data are otherwise untouched.

[0031] Figure 1 demonstrates the data before the step, denoted in blue, and after the Hrj. adaptive algorithm has removed the cardiac artifact, shown in red, across the array for an exemplary participant (NIS004) and a randomly selected five-second time interval. Mid-line channels are highlighted using a dark red color, as seen in the channel guide at the top of the figure. The data demonstrate the difficulty of working with ESG as the cardiac artifact is substantially larger than the underlying data and highlights the effect of the filtering process over the entire array.

[0032] Figures 2-6 demonstrate ESG data during different conditions and different participants. As with the previous figure, mid-line channels are highlighted using a dark red to orient the reader to midline of the array. With the exception of the control (no movement), first movement starts at time t = 5 for each example plot, but sections were randomly selected (i.e., the block for each plotted movement was not fixed).

[0033] Figure 2 shows processed data from a randomly selected 25 second window of nomovement (control). Figure 3 demonstrates ESG data before and during right plantarflexion for an exemplary participant (NIS002). Figure 4 highlights the ESG data during right knee flexion for an example participant (NIS006). Figure 5 depicts changes in the ESG array before and during left plantarflexion. Figure 6 demonstrates changes before and during left knee flexion.

[0034] Figure 7 highlights changes in the ESG recording during the 25 trials, utilizing one second prior to movement onset and one second after movement initiation. For the control (left), two seconds of randomly selected data were chosen to represent each trial. The data were then normalized for that trial, where the normalized amplitude is represented using red shades for positive values, white for zero, and blue shades for negative values. Below each corresponding plot, the average event related potential (ERP) of the 25 trials is shown, where time t = 0 denotes movement onset as found using electromyography data from the soleus muscle for left and right plantarflexion. Movement initiation was defined as the first increase in amplitude larger than two standard deviations from the mean as outlined in the main text. While the experimental protocol was not optimized for ERP analysis, which typically requires 50-100+ trials, the data demonstrate some commonalities around time t — 0, where left plantarflexion tended to cause a positive increase from the mean and right plantarflexion tended toward a negative increase from the mean. Data shown were bandpass filtered between 0.1-20 Hz and are based on bipolar recordings from sensor 14 and sensor 11 (L3 / L2 — L2 / L1).

[0035] Figure 1. Filtering the cardiac artifact via adaptive filtering algorithm across channels: A exemplary five-second section of ESG data pre-H^ filtering shown in blue and post-7 / co filtering shown in red. Channels 2, 5, 8, 11, and 14 are highlighted using a dark red and denote the midline channels in the array as seen in the channel guide (top left). Figure 2. Processed Data Example: A randomly selected 25 second time window of ESG data during the no movement control for participant NIS001 . Mid-line channels are highlighted using red numbers and are denoted in the guide (right).

[0036] Figure 3. Processed Data Example: A 25 second time window of ESG data during right plantarflexion for participant NIS002. First movement occurs at t = 5 seconds. Mid-line channels are highlighted using red numbers and are denoted in the guide (right).

[0037] Figure 4. Processed Data Example: A 25 second time window of ESG data during right knee flexion for participant NIS006. First movement occurs at t = 5 seconds. Mid-line channels are highlighted using red numbers and are denoted in the guide (right).

[0038] Figure 5. Processed Data Example: A 25 second time window of ESG data during left plantarflexion for participant NIS008. First movement occurs at t = 5 seconds. Mid-line channels are highlighted using red numbers and are denoted in the guide (right).

[0039] Figure 6. Processed Data Example: A 25 second time window of ESG data during left knee flexion for participant NISO 10. First movement occurs at t = 5 seconds. Mid-line channels are highlighted using red numbers and are denoted in the guide (right).

[0040] Figure 7. Event Related Potential: Normalized ESG data during two seconds of randomly selected rest (left), one second prior to movement onset and one second postmovement initiation during left plantarflexion (middle) and right plantarflexion (right) are shown at the top portion of the figure, the corresponding average of the trials are shown directly below where the light pink denotes the 95% confidence interval. Individual trials are normalized for that repetition and there are a total of 25 trials for each condition including the control. Data were bandpass filtered between 0.1-20 Hz to highlight low frequency changes and are based on bipolar recordings from sensor 14 and sensor 11 (L3 / L2 — L2 / L1).

[0041] 2. Power Spectrum Changes

[0042] This section contains examples of the power spectral density changes (PSD) changes for several randomly selected participants for left or right knee flexion and left or right pl antarfl exion. Figures 8 and 9 depict examples of PSD changes during knee flexion when compared to the control (rest). Figures 10 through 12 demonstrate exemplary changes in PSD during plantarflexion compared to the control (rest). PSD was calculated using the same methods covered in the methods section of the main text. While a bipolar configuration is utilized for the midline sensors to remove common noise shared between the pair, the EKG sensors - placed at the four comers of the array - are shown using a monopolar arrangement to preserve the spatial distribution of the changes in PSD. The light shaded area represents the 95% confidence interval for that condition.

[0043] Specifically, Figure 8A demonstrates PSD changes during knee flexion, which had similar changes to those shown in the main text for participant NIS007, with decreases in lower frequencies (0.1-20 Hz) for the midline sensors with participant NIS004, while using the bipolar configuration shown by the connecting lines. Fig 8B shows the PSD changes for the EKG sensors placed around the array during the same time periods. In this case, high frequency (30 Hz) increases were noted at all four sensors.

[0044] Figure 9A depicts changes during knee flexion. In this case, low frequency increases were found at the two lower sensor pairs of the array (L1 / L2 — T12 / L1 and L2 / L3 — L1 / L2) with high frequency (20 Hz) changes at the bottom sensor pair (L2 / L3 — L1 / L2) for both left and right knee flexion. Interestingly during the control (rest) the superior sensor pairs (T11 / T12 — T10 / T11 and T12 / L1 — T11 / T12) an increase centered around approximately 20 Hz at rest was noted; however, this was eliminated during knee flexion. The PSD of the EKG sensors placed around the array (Figure 9B) did not demonstrate a significant change when compared to rest.

[0045] Figure 10A shows the PSD changes during pl antarfl exion for participant NIS005. This participants midline sensor pairs (denoted by the lines connecting the sensors) demonstrated a more complex set of changes during movement when compared to control (Rest). For example, left plantarflexion (LPF) PSD calculated at the inferior sensor pairs (L2 / L3 — L1 / L2 and L1 / L2 — T12 / L1) demonstrated broad spectrum decreases while the PSD calculated at the superior-most sensor pair (T11 / T12 — TlO / Tll) increased at frequencies over 5 Hz during right plantarflexion (RPF) only, but this did not occur anywhere else within the array. Instead, there were broad spectrum decreases at sensors T12 / L1 — T11 / T12 during LPF only. The EKG sensors placed around the array (Fig 10B) did not demonstrate any significant changes, except for sensor 2, placed at the superior right of the array (top right) where the PSD decreased during LPF for low-frequencies (0.1-4 Hz).

[0046] Figure 11 A demonstrates the PSD changes during plantarflexion for participant NIS006. Similar to other examples, a low frequency (0.1-20 Hz) drop occurs during movement attempts, in this case the largest differences with respect to the control (Rest) were found at the superior and inferior-most sensor pairs, which showed side- specific decreases with left plantarflexion (LPF) causing a greater decrease in the PSD than right plantarflexion (RPF). There were no changes in the PSD for the “EKG sensors” around the array (Fig 1 IB) during the same time periods except for sensor 2 located at the superior right of the array (top right) where PSD at high frequencies (35 Hz) dropped during RPF only.

[0047] Figure 12A depicts changes in PSD during plantarflexion for participant NIS009.

[0048] Figure 8. PSD Changes during knee flexion: Changes in PSD (A) compared to rest during left knee flexion (LKF) or right knee flexion (RKF) for participant NIS004 midline sensors utilizing a bipolar configuration as denoted by the lines connecting the two sensors demonstrate low frequency (0.1-20 Hz) decreases. In this case, the sensor pair T12 / LI — T11 / T12 (top right) demonstrate side specific changes. When compared to the sensors placed around the array (B) for the same time periods there are no changes at low frequencies, but some high frequency changes (30 Hz) during movement.

[0049] For this participant a broad-spectrum decrease was observed during plantarflexion at all sensor pairs with the exception of L1 / L2 — T12 / L1 (bottom left) which the PSD was found to increase during right plantarflexion (RPF) at frequencies ranging from 0.1-4 Hz. While the EKG sensors only demonstrated a decrease in PSD during left plantarflexion (LPF) at low-frequencies (0.1-5 Hz) at the inferior right sensor (bottom right).

[0050] Figure 9. PSD Changes during knee flexion: Changes in PSD (A) compared to rest during left knee flexion (LKF) or right knee flexion (RKF) for participant NIS008 midline sensors utilizing a bipolar configuration as denoted by the lines connecting the two sensors. In this case, there was a broad-spectrum increase in power when compared to rest. Interestingly, at rest a increase in PSD centered around approximately 20 Hz was eliminated during movement attempts. When compared to the sensors placed around the array (B) for the same time periods there are no changes across the shown frequencies.

[0051] Figure 10. PSD Changes during plantarflexion: Changes in PSD (A) compared to rest during left pl antarfl exion (LPF) or right plantarflexion (RPF) for participant NIS005 midline sensors utilizing a bipolar configuration as denoted by the lines connecting the two sensors. This participant demonstrated several broad spectrum decreases, at the L1 / L2 — T12 / L1 and L2 / L3 — L1 / L2 sensor pairs with an increase during RPF at the Tll / T12 — T10 / T11 sensor pair after approximately 7 Hz. The sensors around the array (B) for the same time periods demonstrated no changes across the shown frequencies with the exception of sensor two (top right) which recorded a low frequency (0.1-4 Hz) drop during LPF only.

[0052] Figure 11. PSD Changes during plantarflexion: Changes in PSD (A) compared to rest during left plantarflexion (LPF) or right pl antarfl exion (RPF) for participant NIS006 midline sensors utilizing a bipolar configuration as denoted by the lines connecting the two sensors. For this participant, during movement (left or right), the largest changes were seen at the extreme ends of the array (711 / 712 — 710 / 711 and L3 / L2 — L1 / L2). The PSD changes at the low frequencies (0.1-12 Hz) show a side specific decrease during movement with the left plantarflexion condition demonstrating a larger decrease than right when compared to the resting state PSD. The sensors placed around the array (B) for the same time periods demonstrated no changes at low frequencies, but did drop during plantarflexion at high frequencies (35 Hz) at sensor 2 (top right).

[0053] Figure 12. PSD Changes during Knee flexion: Changes in PSD (A) compared to rest during left pl antarfl exion (LPF) or right plantarflexion (RPF) for participant NIS009 midline sensors utilizing a bipolar configuration as denoted by the lines connecting the two sensors. In this case, during plantarflexion there were broad-spectrum decreases in PSD with the exception of the 71 / 72 — 712 / 71 sensor pair, which recorded increases during RPF at low frequencies (0.1-4 Hz) and a decrease during LPF from 3 to 40 Hz. The sensors around the array (B) for the same time periods demonstrated no changes at low frequencies with the exception of the inferior right sensor (bottom right) which recorded a decrease in PSD during LPF for frequencies between 0.1-5 Hz.

[0054] Example 2 Abstract.

[0055] Objective. Currently, few non-invasive measures exist for directly measuring spinal sensorimotor networks. Electrospinography (ESG) is one non-invasive method but is primarily used to measure evoked responses or for monitoring the spinal cord during surgery. Our objectives were to evaluate the feasibility of ESG to measure spinal sensorimotor networks by determining spatiotemporal and functional connectivity changes during single-joint movements at the spinal and cortical levels. Approach. We synchronously recorded electroencephalography (EEG), electromyography, and ESG in ten neurologically intact adults while performing one of three lower-limb tasks (no movement, plantar-flexion and knee flexion) in the prone position. A multi-pronged approach was applied for removing artifacts using Hoo filtering, artifact subspace reconstruction and independent component analysis. Next, data were segmented by task and independent components (ICs) of EEG were clustered across participants. Within-participant analysis of ICs and ESG data was conducted, and ESG was characterized in the time and frequency domains. Generalized Partial Directed Coherence (gPDC) analysis was performed within ICs and between ICs and ESG data by participant and task.

[0056] Results. K-means clustering resulted in five clusters of ICs at Brodmann areas (BA) 9, BA 8, BA 39, BA 4, and BA 22. Areas associated with motor planning, working memory, visual processing, movement, and attention, respectively. Time-frequency analysis of ESG data found localized changes during movement execution when compared to no movement. Lastly, we found bidirectional changes in functional connectivity (p < 0.05, adjusted for multiple comparisons) within IC’s and between IC’s and ESG sensors during movement when compared to the no movement condition.

[0057] Significance. To our knowledge this is the first report using high density ESG for characterizing single joint lower limb movements. Our findings provide support that ESG contains information about efferent and afferent signaling in neurologically intact adults and suggests that we can utilize ESG to directly study the spinal cord.

[0058] 1. Introduction

[0059] Spinal cord injury (SCI) is a debilitating and life-long condition affecting approximately 302,000 people living in the United States alone, with an average of 18,000 new cases of SCI each year [1], The costs associated with SCI are substantial, with lifetime cost estimates that can reach $5.8 million (USD) depending on the severity of the injury and age at injury [1], While other areas of healthcare have improved — for example, people with HIV have seen a dramatic increase in life expectancy since the disease was first identified in the 1980s — the same cannot be said for SCI [2], In fact, the average years of life for a person with SCI have not improved since the 1980s and remain significantly lower than that of a person without SCI [1],

[0060] One potential area for improvement is the way SCI is categorized. Currently, the severity of an SCI is ranked using the American Spinal Injury Association Impalement Scale (AIS). This scale uses a binary approach, ranking from AIS ”A”, where a person has no motor or sensory function in the sacral segments, to ”E”, where a person has an SCI but motor and sensory function are normal [3], However, this fails to address the heterogeneity of injury, where each SCI is unique to itself and that even within AIS A patients, different levels of injury (thoracic T2-T10 and thoracolumbar T11-L2) had significantly different motor improvements [4],

[0061] While electroencephalography (EEG) has been used for non-invasive study of the brain and has led to different ways to predict or diagnose disease, seizure, and decoding intent, there are few non-invasive methods to study the spinal cord [5-7], Researchers typically treat the spinal cord as a black box. One method is to apply a known physical or electrical perturbation to the body or spinal cord and relate the response to the stimulus [8, 9], In this case, the output (e.g., movement or changes in movement) is used to infer the underlying motor network behavior. An alternative approach is modeling; for example, high-density electromyography (hdEMG) can be decomposed and related to the activity of spinal motor pools

[0010] , While these approaches can be very effective at characterizing the underlying circuitry of the spinal cord and the behaviors, they are limited in the amount of information provided; namely, they provide information only regarding changes in the spinal cord related to the movement being performed.

[0062] One approach to filling this gap is to use electrospinography (ESG). ESG is a well-studied technique that utilizes one or more surface electrodes placed over different areas of the spinal column and is typically used to record evoked responses from electrical stimulation to peripheral nerves and to monitor the spinal cord during surgery [11, 12], It has been demonstrated that these responses are spinal cord in origin and, while significantly lower in amplitude, have the same characteristics as their invasively recorded counterparts [12, 13], However, little research has been done to determine if ESG can be utilized for other applications such as decoding intent, with only a single study known to the authors that shows a proof of concept

[0014] ,

[0063] In this study, we systematically investigate the use of high-density ESG with neurologically intact participants. First, we recorded EEG, ESG, and electromyography (EMG) simultaneously during single joint movements. We demonstrate that by using an adaptive filtering algorithm, we can effectively remove cardiac artifacts from ESG. Next, we characterize the power spectral density (PSD) changes in the ESG signal during movement and compare them against sensors placed away from the ESG array. Lastly, we look at the changes in functional connectivity between cortical areas and the spinal cord during movement preparation and execution. 2. Methods

[0064] 2. 1. Participants

[0065] 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. Figure 13 provides an overview of the experiment setup (A and B) and outlines the experimental protocol (C) used in this study.

[0066] Figure 13. Experimental protocol: The experiment was performed with (A) the participant in the prone position while measurements were taken with 28-channel electroencephalography (EEG), 4-channel electrooculography (EOG), 15-channel electrospinography (ESG) with 4 channels for noise around the ESG array, and 8-channel electromyography (EMG) / inertial measurement unit (IMU). ESG was recorded at the intervertebral spaces starting at T10 / T11 and ending at L2 / L3 (B) which corresponds to the approximate location of the lumbar enlargement. The participant performed one of three movement tasks (C), a no movement control condition, knee flexion of the left or right leg, and plantarflexion of the left or right foot. These tasks were repeated for 5 blocks of 5 repetitions per block for a total of 25 repetitions per movement per leg, with at least two minutes between conditions and one minute between repetitions. The control condition (no movement) was performed twice for two minutes at the start and end of the movement conditions. Spinal segment chart adapted from

[0015] and utilizes data from

[0016]

[0067] 2.2. Experimental Setup

[0068] 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 1,000 Hz through BrainVision Recorder software (Brain Products GmbH).

[0069] 32 channels of EEG were placed following a modified 10-20 international standard arrangement. Channels TP9, TP 10, 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 TP 10 were placed on the left and right temples, respectively, and P09 and POlO were positioned superior and inferior to the right eye, respectively. Ground and reference were placed on the left and right earlobes, respectively.

[0070] 19 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 AD 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 facedown foam support was used.

[0071] Trigno Avanti wireless surface electromyography (EMG) electrodes (Delsys Inc., USA; common-mode rejection ratio < 80 dB; size: 27 mm 37 mm x 13 mm; input impedance: > 1015 Q / 0.2 pF) were placed at eight different sites: bilaterally over the vastus lateralis (VL), 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 (DOF) inertial measurement unit (IMU) (3-axis accelerometer; sampling resolution: 0.016 ± 0.001 g / bit; noise: < 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 1,925 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.

[0072] 2.3. Experimental Procedure

[0073] Participants were asked to perform one of three tasks: 1) no movement (control), 2) pl antarfl exion 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 j oint 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 seconds after a verbal cue signifying the start of the block. The time between conditions was no less than two minutes, and blocks within a given condition occurred no less than one minute apart. The control task (no movement) was performed for two minutes 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.

[0074] 2.4. Signal Pre-Processing

[0075] Signal pre-processing and statistical analysis were performed off-line utilizing custom software written in MATLAB R2022b (MathWorks, MA) and included functions from the openaccess toolbox EEGLAB

[0018] , A flowchart outlining the signal pre-processing pipeline and analysis is shown in Figure 2 and follows a previously utilized pre-processing methodology [8],

[0076] 2.5. EEG Pre-processing

[0077] EOG data (4 channels) were used to filter the EEG signals (28 channels) via Hco an adaptive algorithm that removes artifacts such as eye blinks or eye movements (q = 10-9, y = 1.15)

[0017] , The data were then zero-phase band-pass (0.1 Hz - 100 Hz) filtered by applying a fourth-order ButterWorth filter.

[0078] 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, two minutes 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.

[0079] Two complementary adaptive filters were applied to remove harmonics of 60 Hz power line noise (CleanLine and ZapLine Plus toolboxes for EEGLAB) [18, 23-25], 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

[0018] , 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. Independent components are latent hidden variables that represent brain components estimated by ‘unmixing’ the EEG recorded on the surface of the scalp under some assumptions. Once you have the ICs, you can use an MRI (corresponding to how gray and white matter, skull / bone, etc are distributed within the head) to estimate the equivalent dipoles corresponding to / or associated with locations of electrical activity causes by brain, eye, muscle movement, etc. Once independent components of significance are identified (and artifactual ICs removed to clean the signals), then IC data can subsequently be analyzed as explained below.

[0080] 2.6. K-means Clustering

[0081] To establish a common group-level comparison of EEG data at the cortical source level, K-means clustering was performed. Specifically, equivalent dipoles of independent components (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 [26-28], 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

[0029] , 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.

[0082] 2. 7. ESG Pre-processing

[0083] Four channels, placed bilaterally 5cm 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 Hco adaptive filtering algorithm, similar to what was used prior for eye artifacts in the EEG recordings.

[0084] In this case, two values q and y were adjusted for this application (q = 10-7, y = 1.2). Specifically, the Hco adaptive weight estimation filter with time varying weights can formulated such that wi+i= wt+PriiP'llrl I. yt (i)

[0085] Where wz+1is the estimated weight per channel at sample i + 1, rtis a vector of reference measurements, in this case the EKG sensors, ytis defined such that where

[0086] P1= P1~ y-2rir , (3) and P1is the noise covariance matrix which satisfies the equation which is initialized with Po= \il where ,z / 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 Hco filter.

[0087] In this formulation the weights are assumed to be time-varying with unknown variation dynamics and when y < 1 the filter provides a guaranteed robust performance under all levels of disturbances in the system and the formulation is valid for the inequality y2< 1 + qf where f 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 y and q, we can modify the behavior of the Hco filter for both applications

[0017] ,

[0088] 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 second 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) [18, 23-25], Lastly, ESG channels were re-referenced by subtracting their common average.

[0089] 2.8. Power Spectral Density

[0090] The data were segmented by condition, then by block, where each condition had five blocks consisting of five repetitions. First, power spectral density (PSD) was calculated with a 750 ms sliding window using 99% overlap to apply Thomson’s multitaper estimate in MATLAB, with 4,990 frequency bins (0.1-50 Hz) and a half-bandwidth product of four for the entirety of each block [8, 30], 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 seconds 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.

[0091] 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.

[0092] 2.9. Functional Connectivity

[0093] Figure 14 (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 [8], Due to the differences in signal characteristics between the cortical sources and ESG recordings, we expected the model order 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.

[0094] Figure 14. Electroencephalograph (EEG) and Electrospinography (ESG) data denoising and analysis: Raw EEG data were cleaned by applying an Hco filter which uses electrooculography (EOG) recordings

[0017] to identify and remove ocular artifacts, signal drifts, and other shared sources of noise. The data are then band-pass filtered (0.1-100 Hz) using a 4th order ButterWorth filter, line noise was filtered via two different complementary adaptive filters, then artifact subspace reconstruction (ASR) was done to remove other types of artifacts. The data were common average referenced, then independent component analysis (ICA) was utilized to further identify and remove any remaining artifacts. Finally, the algorithm DIPFIT was applied to determine the equivalent dipole locations for each IC

[0018] ,

[0095] Raw ESG data were cleaned by applying an Hco filter with different filter weights (q = I O7, y = 1.2) and using recordings around the array to identify and remove cardiac artifacts and other sources of shared artifacts. Then data were band-pass filtered and line noise was filtered utilizing the same method that was applied to EEG, then ASR was performed. Lastly, the data were common average referenced. Raw EMG was downsampled to 1000 Hz, then bandpass filtered using a 4th order ButterWorth (BW) filter to determine the envelope of activation.

[0096] Analysis of the cleaned data included using Bayes information criterion (BIC) [19, 20] to determine the best model order, creating a multivariate autoregressive model (MV AR)

[0021] and calculating the generalized partial directed coherence (gPDC)

[0022] within clustered EEG dipoles and between clustered dipoles and ESG sensors. The EEG data pre-processing (denoising) and analysis methodologies were adapted from [8],

[0097] 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 Bayes information criterion (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.

[0098] 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 [21, 31], The MVAR model can be described by the equation

[0099] (5) in which X(n) is a multivariate M -channel process of length n, A.kare MxM parameter matrices that are comprised of the coefficients a-ij(k)' that relate the ij-th 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 and the model order p is the corresponding model value. To determine information flow, generalized partial directed coherence 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) [22, 321. The gPDC equation is given as (6) where A(f) is the Fourier transform of MVAR model coefficients Ak, A, is the ij-th element of A( ), AT signifies the complex conjugate operation. The difference between gPDC and Partial Directed Coherence (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 crkand signifies the A: -th diagonal element of the covariance matrix 2^. The gPDC output is normalized and therefore bound by the interval [0,1],

[0100] Functional connectivity was used to determine information flow within the participants’ independent components (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.

[0101] 2. JO. Statistical Analysis

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

[0103] Permutation testing was used to determine generalized partial directed coherence (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 [33, 34], To perform permutation testing, the gPDC is first separated into the frequency bands 8 (0.1 - 4 Hz), 0 (4 - 8 Hz), a (8-12 Hz), P (12-30 Hz), and y (30-50 Hz) by condition. Next, the test statistic Tobs, 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.

[0104] The null hypothesis distribution TmMis 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 Tobsand the null distribution Tllullare different.

[0105] Significance was determined by summing the instances when values of Tnu(k') are larger than or equal to Tobs( )and is divided by the total number of permutations [8],

[0106] As multiple comparisons were made, aBonferroni correction was applied and the corrected a, acriticaiwas set at 0.0025. For the gPDC between clusters and ESG sensors, there were k = 40 comparisons performed and aCrtticai was setat0.00125.

[0107] 3. Results

[0108] 3.1. Cardiac Artifact altering

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

[0110] Figure 15. Filtering the cardiac artifact via adaptive filtering algorithm: A exemplary section of ESG data pre-Hco filtering shown in black and post-Hco filtering shown in blue at a time scale of (top left) 50 seconds, (top right) 2 seconds, (bottom) 200 ms with an amplitude scale of ± 100, 50, and 30 pV , respectively. The black square denotes where the data for the next plot, designated with the arrow, was derived.

[0111] 3.2. K-means Clustering

[0112] Figure 16 shows the results of the k-means clustering algorithm, which resulted in five clusters and 107 total dipoles. Figure 16 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 Brodmann areas were identified from the Talairach atlas

[0035] , The Brodmann areas 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.

[0113] Figure 16. Clusters of equivalent dipole sources fit to independent components: Five clusters of dipoles were found across participants, trials, and conditions. Brodmann areas are the regions found within a ± 5 mm search range of the cluster centroids and had a maximum distance of 3 cm in this study. Cluster number is for referencing of individual cluster location in the text. The left side of the transverse view is denoted with a L on cluster 1.

[0114] 3.3. ESG PSD Analysis

[0115] Figure 17A demonstrates power spectral density (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.

[0116] Analysis of the EKG sensors as seen in Figure 17B 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.

[0117] Figure 17. Power spectral density (PSD) changes during knee flexion: Average changes in PSD during right (RKF) and left knee flexion (LKF) when compared to control (Rest) (A) utilizing a bipolar configuration where the inferior sensor is subtracted from the superior and is denoted by the lines connecting sensors. Changes in PSD of the EKG sensors, or sensors placed around the array, (B) from the same participant (NIS007) and time period. Light shaded area represents the 95% confidence interval for that condition. Additional examples of PSD changes during knee flexion or plantarflexion can be found in the supplementary materials 3.4. Functional Connectivity

[0118] Figure 18 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 (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 grey 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.

[0119] Figure 18. Functional connectivity changes over time: Significant differences in functional connectivity (gPDC) during right plantarflexion for a representative participant. Data shown are from a single block and organized by band (top) corresponding cluster data (bottom left) and ESG data (middle). The grey shaded area of the EMG data from the right soleus (bottom right) is the two second section of data being shown which is (A) just prior to first movement and (B) during first movement initiation. Cluster number location is shown in the delta band panel (top left). See Supplementary Materials for videos demonstrating changes in functional connectivity over time for both plantarflexion and knee flexion conditions.

[0120] 4. Discussion

[0121] In this study we identified five common clusters of EEG activity among ten neurologically intact participants. Specifically, BA 9 (cluster 1) has been shown to be involved with motor planning and working memory. BA 8 (cluster 2) has been demonstrated to be important for visual attention. BA 19 (cluster 3) is associated with processing visual information. BA 4 (cluster 4) is also known as the primary motor cortex and is responsible for execution of movements. BA 39 has been shown to play a role in spatial cognition, memory retrieval, and attention. Each cluster contained at least one IC from each participant, with the exception of one participant who did not have IC’s that fell within the designated range for clusters 3 and 5.

[0122] This is the first study demonstrating the feasibility of high-density ESG to characterize single joint movements from changes in the electrical activity at the level of the spinal cord during movement. While caution must be taken in interpreting the results, we demonstrated that sensors placed bilaterally approximately 5 cm from midline at the top and bottom of the lumbar ESG array did not show changes in power spectral density during movement that were comparable with the analysis performed using the midline sensors in the array. The experiment was conducted in a controlled environment to minimize the possibility of sensor movement during the two tasks performed and, in most cases, the changes were localized within the array in the superior-inferior direction. Furthermore, it has been demonstrated that peripheral nerve stimulation can evoke a response from the spinal cord that can be recorded utilizing surface electrodes with only one paper known to the authors showing proof of concept for changes in surface recordings during movement using a similar methodology [13, 14],

[0123] Unfortunately, this study ultimately cannot conclude that the recorded signals originated from the spinal cord. Nevertheless, we found that the power spectral changes behave similarly to invasively recordings done using animal models. For example, it has been demonstrated that there is a decrease in 8-30 Hz activity recorded from the epidural space in the skull before and during forelimb movement in rats

[0036] , Another study recording activity from the lumbar spinal cord using micro wires found low frequencies (0.5-30 Hz) and high frequencies ( 100 Hz) conveyed more information about hind limb movement in cats

[0037] , Additionally, gPDC analysis suggests that the ESG recordings are bidirectional, which further reduces the likelihood that the data are from some other physiological source such as EMG which should cause mono-directional changes in gPDC (i.e., from a cluster location to the muscle being activated). Our findings suggest that the information provided by ESG is localized to the ESG array and provides useful information despite unknown origins.

[0124] There are several limitations to the study and one is due to the anatomical differences between individuals where the spinal cord can terminate up to two vertebral levels from the mean location, thus averaging the data was not best practice

[0038] , Future work could use MRI data to standardize locations between participants so that a group average could be performed, however even this may be challenging due to the anatomical differences of the vertebra and the conductivity of the tissues, which may need to be taken into account for this type of analysis. Nevertheless, individual variability in the patterns of ESG measurements are likely to have at least some diagnostic value. Another limitation of the proposed noninvasive technique is that the ESG electrodes are well positioned to record mainly the posterior funiculus, which contains ascending pathways (dorsal column fibers) carrying information concerning touch and limb position from the body to the brain. Thus, the ESG technique seems less ideal to record descending motor information to the muscles. However, we suspect that some of these descending signals may be registered indirectly by the ESG electrodes via volume conduction, particularly prior to movement when the spinal cord may be “quiet”. Interestingly, the gPDC findings (see Figure 18A, alpha band) shows significant changes in functional connectivity in alpha, beta and gamma bands directed from the brain to the spinal cord (that is, descending information) prior to movement. Conversely, gPDC changes directed from the spinal cord to the brain (i.e., ascending information) in the delta, beta and gamma bands are shown during the initial phases of movement (Figure 18B). Future studies should use source analysis and an appropriate spinal cord model and lead field matrix to estimate the sources of these ESG signals within the spinal cord during volitional (self-initiated or attempted) lower- limb movements

[0039] ,

[0125] Due to the low conductivity of fatty tissues, we hypothesize that body composition is one of the largest limitations to the usefulness of ESG in adults and while this was not a factor in the participant selection for this study, the body composition of the participants was well within normal ranges

[0040] ,

[0126] Due to preparatory contraction of back muscles just prior to and during movement / postural adjustments, the experiments were performed in the prone position to minimize or eliminate the influence of muscle artifact contamination on the ESG signal. Further studies should be done to determine if ESG can be recorded in the sitting or standing positions, which would further expand uses for the recording modality. Additionally, while this experiment was performed using recordings over the approximate location of the lumbar enlargement during lower limb movements, future work with recordings over the approximate location of the cervical enlargement may be useful to determine if changes can be detected during arm and / or hand movement(s) and if different fine movements can be discriminated.

[0127] Another limitation is the low signal to noise ratio which was found to be at best similar to EEG and could be made worse due to spinal cord injury, which is thought to decrease the amplitude of spared descending signals

[0041] , However, previous work has demonstrated that transcranial magnetic stimulation (TMS) in conjunction with voluntary effort can create motor evoked potentials that can be recorded on the surface of the skin over the spinal column above, at, and below the lesion in individuals who are motor and sensory complete (ATS A)

[0042] , This suggests that ESG recordings may provide useful information even for people with a SCI where the signal is degraded.

[0128] Future work will be done with people who have a SCI to determine if ESG is sensitive enough to detect changes in descending drive (i.e., attempting movement) or for measuring the resting state of the spinal cord after injury, which may provide valuable insights into adaptations made by the spinal cord and how to treat them. If validated, longitudinal studies could also be done with participants with an SCI while undergoing therapy to determine if ESG can be used as a metric for recovery.

[0129] 5. Conclusion

[0130] In this study we utilized a combination of EEG, ESG, and EMG in neurologically intact participants ( V = 10) while performing knee flexion or plantarflexion in the prone position. To the authors knowledge, this is the first time ESG has been utilized in this method and the first study demonstrating the feasibility of ESG to determine movement in neurologically intact participants. ESG was demonstrated to have localized changes in power spectral density during movement when compared to rest and these changes were not found with sensors placed bilaterally 5 cm from the top and bottom of the array. We then demonstrated bidirectional changes in functional connectivity (increases and decreases) during movement. The findings of this study could aid in uncovering fundamental knowledge about the intact spinal cord and how to better treat SCI.

[0131] Note that the method described above is applicable to various muscles and muscle groups including upper boy muscles.

[0132] A method is described herein comprising receiving an electroencephalography (EEG) signal from a plurality of EEG sensors during performance of a plurality of tasks, wherein the plurality of EEG sensors are located on scalps of participants, receiving an electrospinography (ESG) signal from a plurality of ESG sensors during the performance of a plurality of tasks, wherein the plurality of ESG sensors are attached along a spinal region of each participant, receiving an electromyography data (EMG) signal from a plurality of EMG sensors during performance of the plurality of tasks, wherein the plurality of EMG sensors are positioned at locations on each participant’s body, wherein the receiving the EEG signal, the ESG signal, and the EMG signal is synchronous, filtering the EEG signal, the ESG signal, and the EMG signal;, identifying independent components in the EEG signal for each participant using independent component analysis, wherein the independent components correspond to dipole activity associated with participants’ performance of the plurality of tasks, and analyzing the independent components and corresponding ESG signal to determine associations among participant independent components and between each participants independent components and respective ESG signal.

[0133] In embodiments the plurality of EEG sensors comprise four (4) electrooculogram (EOG) data channels.

[0134] In embodiments the filtering the EEG signal comprises applying an Hco adaptive algorithm to remove artifacts using information of the EOG data channels.

[0135] In embodiments the artifacts include eye blinks and eye movements.

[0136] In embodiments the filtering the EEG signal comprises applying a band pass filter.

[0137] In embodiments the filtering the EEG signal comprises applying artifact subspace reconstruction.

[0138] In embodiments the filtering the EEG signal comprises applying adaptive filters to remove harmonics of power line noise.

[0139] In embodiments the filtering the EEG signal comprises re-referencing the plurality of EEG sensors by subtracting their common average.

[0140] In embodiments the plurality of ESG sensors comprise four (4) electrocardiogram (EKG) data channels.

[0141] In embodiments the filtering the ESG signal comprises applying an Hco adaptive algorithm to remove artifacts using information of the EKG data channels.

[0142] In embodiments the artifacts include cardiac artifacts.

[0143] In embodiments the filtering the ESG signal comprises applying a band pass filter.

[0144] In embodiments the filtering the ESG signal comprises applying artifact subspace reconstruction.

[0145] In embodiments the filtering the ESG signal comprises applying adaptive filters to remove harmonics of power line noise.

[0146] In embodiments the filtering the ESG signal comprises re-referencing the plurality of ESG sensors by subtracting their common average.

[0147] In embodiments the filtering the EMG signal comprises down sampling to 1000Hz. In embodiments the filtering the EMG signal comprises applying a high pass filter.

[0148] In embodiments the filtering the EMG signal comprises applying a low pass filter.

[0149] In embodiments the analyzing comprises clustering each participant’s independent components using K-means clustering.

[0150] In embodiments the analyzing comprises segmenting each participant’s filtered ESG signal and independent component data by task of the plurality of tasks.

[0151] In embodiments the analyzing further comprises segmenting each participant’s filtered ESG signal and independent component data using information of the corresponding EMG signal.

[0152] In embodiments the analyzing comprises applying a generalized partial directed coherence model to within participant segmented independent cluster data.

[0153] In embodiments the analyzing comprises applying a generalized partial directed coherence model to each participant’s segmented independent cluster data and respective segmented ESG data.

[0154] In embodiments the plurality of EEG sensors are positioned according to a modified 10- 20 international standard arrangement.

[0155] In embodiments wherein the spinal region comprises intervertebral spaces, starting at T10 / T 11 through L2 / L3.

[0156] In embodiments the locations comprise sites bilaterally located over the vastus lateralis (VL), the medial hamstring (MH), the tibialis anterior (TA), and the lateral aspect of the soleus (SOL) muscles.

[0157] References:

[0158] [1] NSCISC, “National spinal cord injury statistical center (facts and figures at a glance),” 2023.

[0159] [2] A. Trickey, C. A. Sabin, G. Burkholder, H. Crane, A. d. Monforte, M. Egger, M. J. Gill, S. Grabar, J. L. Guest, I. Jarrin et al., “Life expectancy after 2015 of adults with hiv on longterm antiretroviral therapy in europe and north america: a collaborative analysis of cohort studies,” The Lancet HIV, 2023.

[0160] [3] T. T. Roberts, G. R. Leonard, and D. J. Cepela, “Classifications in brief: American spinal injury association (asia) impairment scale,” 2017.

[0161] [4] M. F. Dvorak, V. K. Noonan, N. Fallah, C. G. Fisher, C. S. Rivers, H. Ahn, E. C. Tsai, A. Linassi, S. D. Christie, N. Attabib etal., “Minimizing errors in acute traumatic spinal cord injury trials by acknowledging the heterogeneity of spinal cord anatomy and injury severity: an observational Canadian cohort analysis,” Journal of neurotrauma, vol. 31, no. 18, pp. 1540-1547, 2014.

[0162] [5] S. Gaubert, F. Raimondo, M. Houot, M.-C. Corsi, L. Naccache, J. Diego Sitt, B. Hermann, D. Oudiette, G. Gagliardi, M.-O. Habert et al., “Eeg evidence of compensatory mechanisms in preclinical alzheimer’s disease,” Brain, vol. 142, no. 7, pp. 2096-2112, 2019.

[0163] [6] A. G. Steele, S. Parekh, H. F. Azgomi, M. B. Ahmadi, A. Craik, S. Pati, J. T. Francis, J. L. Contreras- Vidal, and R. T. Faghih, “A mixed filtering approach for real-time seizure state tracking using multi-channel electroencephalography data,” IEEE Transactions on Neural Systems and Rehabilitation Engineering, vol. 29, pp. 2037-2045, 2021.

[0164] [7] A. Kilicarslan, S. Prasad, R. G. Grossman, and J. L. Contreras- Vidal, “High accuracy decoding of user intentions using eeg to control a lower-body exoskeleton,” in 2013 35th annual international conference of the IEEE Engineering in Medicine and Biology Society (EMBC). IEEE, 2013, pp. 5606-5609.

[0165] [8] A. G. Steele, D. A. Atkinson, B. Varghese, J. Oh, R. L. Markley, and D. G. Sayenko, “Characterization of spinal sensorimotor network using transcutaneous spinal stimulation during voluntary movement preparation and performance,” Journal of Clinical Medicine, vol. 10, no. 24, p. 5958, 2021.

[0166] [9] J. Weiler, P. L. Gribble, and J. A. Pruszynski, “Spinal stretch reflexes support efficient hand control,” Nature neuroscience, vol. 22, no. 4, pp. 529-533, 2019.

[0167]

[0010] R. Ornelas-Kobayashi, A. Gogeascoechea, and M. Sartori, “Person-specific biophysical modeling of alpha- motoneuron pools driven by in vivo decoded neural synaptic input,” IEEE transactions on neural systems and rehabilitation engineering, vol. 31, pp. 1532— 1541, 2023.

[0168]

[0011] B. Ferrero and A. Di Liberto, “Spinal cord injury: role of neurophysiology,” in Spinal Cord Injury (SCI) Repair Strategies. Elsevier, 2020, pp. 13-37.

[0169]

[0012] D. Jeanmonod, M. Sindou, and F. Mauguiere, “The human cervical and lumbo-sacral evoked electrospinogram. data from intra-operative spinal cord surface recordings,” Electroencephalography and Clinical Neurophysiology / Evoked Potentials Section, vol. 80, no. 6, pp. 477-489, 1991.

[0170]

[0013] M. Dimitrijevic, L. Lehmkuhl, E. Sedgwick, A. Sherwood, and W. McKay, “Characteristics of spinal cord- evoked responses in man,” Stereotactic and Functional Neurosurgery, vol. 43, no. 3-5, pp. 118-127, 1980.

[0014] E. Yom-Tov and G. Inbar, “Movement-related potentials in the human spinal cord preceding toe movement,” Clinical neurophysiology, vol. I l l, no. 2, pp. 350-361, 2000.

[0171]

[0015] D. G. Sayenko, D. A. Atkinson, C. J. Dy, K. M. Gurley, V. L. Smith, C. Angeli, S. J.

[0172] Harkema, V. R. Edgerton, and Y. P. Gerasimenko, “Spinal segment-specific transcutaneous stimulation differentially shapes activation pattern among motor pools in humans,” Journal of Applied Physiology, vol. 118, no. 11, pp. 1364- 1374, 2015.

[0173]

[0016] F. P. Kendall, E. K. McCreary, P. G. Provance, M. M. Rodgers, and W. A. Romani, “Muscles: Testing and function, with posture and pain,” Lippincott Williams and Wilkins, Pennsylvania, determination of footedness. J. Phys. Med. Rehabil, vol. 2, pp. 835-841, 1983.

[0174]

[0017] A. Kilicarslan, R. G. Grossman, and J. L. Contreras-Vidal, “A robust adaptive denoising framework for real-time artifact removal in scalp eeg measurements,” Journal of neural engineering, vol. 13, no. 2, 2016.

[0175]

[0018] A. Delorme and S. Makeig, “Eeglab: an open source toolbox for analysis of single-trial eeg dynamics including independent component analysis,” Journal of neuroscience methods, vol. 134, no. 1, pp. 9-21, 2004.

[0176]

[0019] A. A. Neath and J. E. Cavanaugh, “The bayesian information criterion: background, derivation, and applications,” Wiley Interdisciplinary Reviews: Computational Statistics, vol. 4, no. 2, pp. 199-203, 2012.

[0177]

[0020] S. Mousavi, M. Niknazar, and B. V. Vahdat, “Epileptic seizure detection using ar model on eeg signals,” in 2008 Cairo International Biomedical Engineering Conference. IEEE, 2008, pp. 1-4.

[0178]

[0021] T. Schneider and A. Neumaier, “Algorithm 808: Arfit — a matlab package for the estimation of parameters and eigenmodes of multivariate autoregressive models,” ACM Transactions on Mathematical Software (TOMS), vol. 27, no. 1, pp. 58-65, 2001.

[0179]

[0022] L. A. Baccala, K. Sameshima, and D. Y. Takahashi, “Generalized partial directed coherence,” in 200715th International conference on digital signal processing. IEEE, 2007, pp. 163-166.

[0180]

[0023] T. Mullen, “Cleanline eeglab plugin,” San Diego, CA: Neuroimaging Informatics Toolsand Resources Clearinghouse (NITRC), 2012.

[0181]

[0024] M. Klug and N. A. Kloosterman, “Zapline-plus: A zapline extension for automatic and adaptive removal of frequency-specific noise artifacts in m / eeg,” Human Brain Mapping, vol. 43, no. 9, pp. 2743-2758, 2022.

[0182]

[0025] M. Miyakoshi, L. M. Schmitt, C. A. Erickson, J. A. Sweeney, and E. V. Pedapati, “Can we push the “quasi-perfect artifact rejection” even closer to perfection?” Frontiers in Neuroinformatics, vol. 14, p. 597079, 2021.

[0183]

[0026] D. L. Davies and D. W. Bouldin, “A cluster separation measure,” IEEE transactions on pattern analysis and machine intelligence, no. 2, pp. 224-227, 1979.

[0184]

[0027] P. J. Rousseeuw, “Silhouettes: a graphical aid to the interpretation and validation of cluster analysis,” Journal of computational and applied mathematics, vol. 20, pp. 53-65, 1987.

[0185]

[0028] T. Calihski and J. Harabasz, “A dendrite method for cluster analysis,” Communications in Statistics-theory and Methods, vol. 3, no. 1, pp. 1-27, 1974.

[0186]

[0029] N. Kriegeskorte, W. K. Simmons, P. S. Bellgowan, and C. I. Baker, “Circular analysis in systems neuroscience: the dangers of double dipping,” Nature neuroscience, vol. 12, no. 5, pp. 535-540, 2009.

[0030] D. J. Thomson, “Spectrum estimation and harmonic analysis,” Proceedings of the IEEE, vol. 70, no. 9, pp. 1055-1096, 1982.

[0187]

[0031] A. Schldgl and G. Supp, “Analyzing event-related eeg data with multivariate autoregressive parameters,” Progress in brain research, vol. 159, pp. 135-147, 2006.

[0188]

[0032] L. A. Baccala and K. Sameshima, “Partial directed coherence: a new concept in neural structure determination,” Biological cybernetics, vol. 84, no. 6, pp. 463-474, 2001.

[0189]

[0033] E. Maris and R. Oostenveld, “Nonparametric statistical testing of eeg-and meg-data,” Journal of neuroscience methods, vol. 164, no. 1, pp. 177-190, 2007.

[0190]

[0034] F. Pesarin and L. Salmaso, “The permutation testing approach: a review,” Statistica, vol. 70, no. 4, pp. 481- 509, 2010.

[0191]

[0035] J. L. Lancaster, M. G. Woldorff, L. M. Parsons, M. Liotti, C. S. Freitas, L. Rainey, P. V. Kochunov, D. Nickerson, S. A. Mikiten, and P. T. Fox, “Automated talairach atlas labels for functional brain mapping,” Human brain mapping, vol. 10, no. 3, pp. 120-131, 2000.

[0192]

[0036] M. W. Slutzky, L. R. Jordan, E. W. Lindberg, K. E. Lindsay, and L. E. Miller, “Decoding the rat forelimb movement direction from epidural and intracortical field potentials,” Journal of neural engineering, vol. 8, no. 3, p. 036013, 2011.

[0193]

[0037] Y. Fathi and A. Erfanian, “Decoding bilateral hindlimb kinematics from cat spinal signals using three- dimensional convolutional neural network,” Frontiers in Neuroscience, vol. 16, 2022.

[0194]

[0038] A. Macdonald, P. Chatrath, T. Spector, and H. Ellis, “Level of termination of the spinal cord and the dural sac: a magnetic resonance study,” Clinical Anatomy, vol. 12, no. 3, pp. 149- 152, 1999.

[0195]

[0039] Y. Adachi, M. Miyamoto, J. Kawai, G. Uehara, S. Tomizawa, and S. Kawabata, “Source analysis of cervical spinal cord evoked field detected by vector squid sensors,” in International Congress Series, vol. 1300. Elsevier, 2007, pp. 557-560.

[0196]

[0040] G. R. Hernandez-Labrado, J. L. Polo, E. L6pez-Dolado, and J. E. Collazos-Castro, “Spinal cord direct current stimulation: finite element analysis of the electric field and current density,” Medical & biological engineering & computing, vol. 49, pp. 417-429, 2011.

[0197]

[0041] H. J. Sfreddo, J. R. Wecht, O. A. Alsalman, Y.-K. Wu, and N. Y. Harel, “Duration and reliability of the silent period in individuals with spinal cord injury,” Spinal Cord, vol. 59, no. 8, pp. 885-893, 2021.

[0198]

[0042] P. H. Ellaway, M. Catley, N. J. Davey, A. Kuppuswamy, P. Strutton, H. L. Frankel, A. Jamous, and G. Savic, “Review of physiological motor outcome measures in spinal cord injury using transcranial magnetic stimulation and spinal reflexes,” Journal of rehabilitation research and development, vol. 44, no. 1, p. 69, 2007.

Claims

CLAIMS1. A method comprising, receiving an electroencephalography (EEG) signal from a plurality of EEG sensors during performance of a plurality of tasks, wherein the plurality of EEG sensors are located on scalps of participants; receiving an electrospinography (ESG) signal from a plurality of ESG sensors during the performance of a plurality of tasks, wherein the plurality of ESG sensors are attached along a spinal region of each participant; receiving an electromyography data (EMG) signal from a plurality of EMG sensors during the performance of the plurality of tasks, wherein the plurality of EMG sensors are positioned at locations on each participant’s body, wherein the receiving the EEG signal, the ESG signal, and the EMG signal is synchronous; filtering the EEG signal, the ESG signal, and the EMG signal; identifying independent components in the EEG signal for each participant using independent component analysis, wherein the independent components correspond to dipole activity associated with participants’ performance of the plurality of tasks; and analyzing the independent components and corresponding ESG signal to determine associations among participant independent components and between each participants independent components and respective ESG signal.

2. The method of claim 1, wherein the plurality of EEG sensors comprise four (4) electrooculogram (EOG) data channels.

3. The method of claim 2, wherein the filtering the EEG signal comprises applying an Hco adaptive algorithm to remove artifacts using information of the EOG data channels.

4. The method of claim 3, wherein the artifacts include eye blinks and eye movements.

5. The method of claim 4, wherein the filtering the EEG signal comprises applying a band pass filter.

6. The method of claim 5, wherein the filtering the EEG signal comprises applying artifact subspace reconstruction.

7. The method of claim 6, wherein the filtering the EEG signal comprises applying adaptive filters to remove harmonics of power line noise.

8. The method of claim 7, wherein the filtering the EEG signal comprises re-referencing the plurality of EEG sensors by subtracting their common average.

9. The method of claim 1, wherein the plurality of ESG sensors comprise four (4) electrocardiogram (EKG) data channels.

10. The method of claim 9, wherein the filtering the ESG signal comprises applying an Hco adaptive algorithm to remove artifacts using information of the EKG data channels.

11. The method of claim 10, wherein the artifacts include cardiac artifacts.

12. The method of claim 11, wherein the filtering the ESG signal comprises applying a band pass filter.

13. The method of claim 12, wherein the filtering the ESG signal comprises applying adaptive filters to remove harmonics of power line noise.

14. The method of claim 13, wherein the filtering the ESG signal comprises applying artifact subspace reconstruction.

15. The method of claim 14, wherein the filtering the ESG signal comprises re-referencing the plurality of ESG sensors by subtracting their common average.

16. The method of claim 1, wherein the filtering the EMG signal comprises down sampling to 1000Hz.

17. The method of claim 16, wherein the filtering the EMG signal comprises applying a high pass filter.

18. The method of claim 17, wherein the filtering the EMG signal comprises applying a low pass filter.

19. The method of claim 1, wherein the analyzing comprises clustering each participant’s independent components using K-means clustering.

20. The method of claim 19, wherein the analyzing comprises segmenting each participant’s filtered ESG signal and independent component data by task of the plurality of tasks.

21. The method of claim 20, wherein the analyzing comprises computing power spectral data of the ESG signal.

22. The method of claim 21, wherein the analyzing further comprises segmenting each participant’s filtered ESG signal and independent component data using information of the corresponding EMG signal.

23. The method of claim 22, wherein the analyzing comprises applying a generalized partial directed coherence model to within participant segmented independent cluster data.

24. The method of claim 23, wherein the analyzing comprises applying a generalized partial directed coherence model to each participant’s segmented independent cluster data and respective segmented ESG data.

25. The method of claim 1, wherein the plurality of EEG sensors are positioned according to a modified 10-20 international standard arrangement.

26. The method of claim 1, wherein the spinal region comprises intervertebral spaces, starting at T10 / T11 through L2 / L3.

27. The method of claim 1, wherein the locations comprise sites bilaterally located over the vastus lateralis (VL), the medial hamstring (MH), the tibialis anterior (TA), and the lateral aspect of the soleus (SOL) muscles.

Images