METHODS OF DETECTING AND MAPPING SPINAL CORD EPIDURALLY-EVOKED POTENTIALS

Abstract

Embodiments of the present invention relate to a novel approach to automatically detect the occurrence of evoked potentials, quantify the attributes of the signal and visualize the effect across a high number of spinal cord epidural stimulation parameters. This new method is designed to automate the current process for performing this task that has been accomplished manually by data analysts through observation of the raw EMG signals, which is laborious and time-consuming as well as being prone to human errors. The proposed method provides fast and accurate framework for activation detection and visualization of the results within five main algorithms.

Claims

1. A method for detection and characterization of evoked potentials in muscles, comprising: recording electrical signals at least one of a patient's muscles; dividing the electrical signals into sequential segments; calculating, for each segment, a highest statistical difference Si between the electrical signal in that segment and background noise; calculating a dynamic threshold h; and detecting an evoked potential in one of the patient's muscles if Si>h in at least 50% of the segments.

2. The method of claim 1, further comprising delivering a plurality of epidural stimulations to the patient's spinal cord at a first stimulation voltage, each of the plurality of epidural stimulations separated by a time interval.

3. The method of claim 2, wherein the delivering and the recording occur concurrently.

4. The method of claim 2, wherein each of the sequential segments has a duration equal to the time interval.

5. The method of claim 2, wherein each segment temporally overlaps with delivery of one of the plurality of epidural stimulations.

6. The method of claim 1, wherein the electrical signals are electromyography signals.

7. The method of claim 1, further comprising reducing the noise of the electrical signals.

8. The method of claim 7, wherein reducing the noise comprises: converting each electrical signal into a 2-D image; and applying an image smoothing method to the 2-D image to reduce signal noise.

9. The method of claim 1, wherein Si is determined by $S_i={.Math.}_{k=j}^{N_i}.Math.s_k=\left(N_i-j+1\right).Math.\ln \left(\frac{{}_0}{{}_i}\right)+\frac{\left(N_i-j\right)}{2}.Math.\left(\frac{{}_i^2}{{}_0^2}-1\right).$

10. The method of claim 1, wherein his determined by
h=S.sub.max+.sub.S.sub.i.

11. The method of claim 1, further comprising delivering a plurality of epidural stimulations to the patient's spinal cord at a first stimulation voltage, followed by delivering a plurality of epidural stimulations to the patient's spinal cord at a second stimulation voltage, wherein the second stimulation voltage is unequal to the first stimulation voltage.

12. The method of claim 11, wherein the second stimulation voltage is higher than the first stimulation voltage.

13. The method of claim 1, further comprising extracting at least one feature of the evoked potential.

14. The method of claim 13, wherein the feature is one of a peak-to-peak interval, a min-max interval, an activation latency, and an integrated EMG.

15. The method of claim 13, further comprising displaying the at least one feature in a visually perceptible format.

16. The method of claim 15, wherein the visually perceptible format is a colormap image.

17. The method of claim 1, wherein recording electrical signals at least one of the patient's muscles includes recording electrical signals at a plurality of the patient's muscles, and wherein detecting an evoked potential in one of the patient's muscles includes associating the evoked potential with one of the patient's muscles in the plurality of the patient's muscles.

18. A method for detection and characterization of evoked potentials in muscles, comprising: delivering a plurality of epidural stimulations to a patient at a stimulation voltage, each of the plurality of epidural stimulations separated by a time interval; recording electrical signals at least one of the patient's muscles, wherein the recording is concurrent with the delivering; dividing the electrical signals into sequential segments; calculating, for each segment, a highest statistical difference Si between the electrical signal in that segment and background noise; calculating a dynamic threshold h; determining whether Si>h in at least 50% of the segments; and increasing the stimulation voltage until an evoked potential is detected; wherein the evoked potential is detected if Si>h in at least 50% of the segments.

19. The method of claim 18, further comprising extracting at least one feature of the evoked potential.

20. The method of claim 19, wherein the feature is one of a peak-to-peak interval, a min-max interval, an activation latency, and an integrated EMG.

21. The method of claim 18, further comprising reducing the noise of the electrical signals.

22. The method of claim 21, wherein reducing the noise comprises: converting each electrical signal into a 2-D image; and applying an image smoothing method to the 2-D image to reduce signal noise.

23. A method for determining the location a muscle activation potential evoked from epidural stimulation and minimum voltage required to evoke activation potential, comprising: delivering a plurality of epidural stimulations to a patient at a stimulation voltage; recording electrical signals at a plurality of the patient's muscles; dividing the electrical signals into sequential segments; reducing the noise of the electrical signals; calculating, for each segment, a highest statistical difference Si between the electrical signal in that segment and background noise; calculating a dynamic threshold h; determining whether Si>h in at least 50% of the segments; and increasing the stimulation voltage until an activation potential is detected in one of the plurality of muscles, wherein the activation potential is detected if Si>h in at least 50% of the segments.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0011] A better understanding of the present invention will be had upon reference to the following description in conjunction with the accompanying drawings.

[0012] FIG. 1 depicts a block diagram of the proposed framework for visualization and onset detection of muscle activations induced by epidural spinal cord stimulation in individuals with spinal cord injury.

[0013] FIGS. 2A-D depict the steps for converting raw EMG signals into 2-D and 3-D images. FIG. 2A is a graph depicting a raw EMG signal. FIG. 2B is a graph illustrating segmentation of the raw EMG signal using stimulation time intervals. FIG. 2C is a 3-D graph depicting the overlaid segments, wherein the X-axis is the evoked potentials duration (s), the Y-axis is the stimulation voltage (V), and the Z-axis is the amplitude of the signals (V). FIG. 2D is a 2-D Colormap image depicting the data of the 3-D graph in FIG. 2C, wherein the X-axis is the evoked potentials duration (s), the Y-axis is the stimulation voltage (V), and the amplitude of the signals (V) is indicated by the color of the image.

[0014] FIGS. 3A-B depict application of the GGMRF method to reduce signal noise. FIG. 3A depicts the 2-D Colormap image of FIG. 2D with an overlaid representation of the application of the GGMRF method. FIG. 3B is a graph depicting the EMG signal before and after application of the GGMRF method.

[0015] FIGS. 4A-D depict calculation steps for the activation detection algorithm. FIG. 4A is a graph depicting an evoked potential in a single segment of the EMG signal. FIG. 4B is a graph depicting the evoked potential of FIG. 4A with a histogram of its estimated Gaussian distribution. FIG. 4C is a graph comparing the Gaussian probability density function (PDF) of evoked potential signal to the PDF of background noise. FIG. 4D is a graph plotting the calculated the log-likelihood ratio (LLR) for each segment of the EMG signal and determining the activation threshold.

[0016] FIGS. 5A-B depict selected feature parameters of the EMG activation signal. FIG. 5A depicts the Colormap image of FIG. 2D with the number of peaks of the evoked potential, activation onset and latency identified by visual inspection. FIG. 5B is a graph depicting features extracted by an embodiment of the automated method of the present invention, namely, the peak-to-peak value (Vpp), activation latency, time interval between minimum and maximum values (Tpp) and integrated EMG (summation of absolute values of all gray areas).

[0017] FIG. 6 depicts the conversion of 14 muscle EMG signals into Colormap images. Section A of FIG. 6 depicts raw EMG signals of 14 proximal and distal leg muscles. Section B of FIG. 6 is a Colormap image depicting the corresponding peak-to-peak values (v) of the signals in each stimulation voltage after stimulation threshold detection. The white area in the images indicate the stimulation voltage at which the experiment began and the gray area indicates the stimulation voltages below the stimulation threshold.

[0018] FIGS. 7A-D depict box-plot representations of performance measurements for comparing the disclosed automated activation detection method with manual detection by analysts. FIG. 7A depicts a comparison of accuracy. FIG. 7B depicts a comparison of sensitivity. FIG. 7C depicts a comparison of specificity. FIG. 7D depicts a comparison of Dice similarity.

[0019] FIG. 8 is a graph depicting a comparison of three activation detection methods, TKEO, AGLR, and AGLR with GGMRF, for accuracy as an activation detection method as a function of SNR (dB).

[0020] FIG. 9A-C depict Colormap images output from the disclosed method. FIG. 9A depicts the peak-to-peak amplitude of detected evoked potentials as a function of muscles vs. stimulation voltage for a given frequency and configuration. FIG. 9B depicts the peak-to-peak amplitude as a function of configuration vs. stimulation voltage for a given muscle and frequency. FIG. 9C depicts the integrated EMG values as a function of muscles vs. stimulation frequency for a given voltage and configuration. In these images, the gray area represents the pre-threshold part of the experiment where no activation was induced. In the graphical representation of the configurations, cathodes are black, anodes are gray and un-assigned electrodes are white.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0021] For the purposes of promoting an understanding of the principles of the invention, reference will now be made to selected embodiments illustrated in the drawings and specific language will be used to describe the same. It will nevertheless be understood that no limitation of the scope of the invention is thereby intended; any alterations and further modifications of the described or illustrated embodiments, and any further applications of the principles of the invention as illustrated herein are contemplated as would normally occur to one skilled in the art to which the invention relates. At least one embodiment of the invention is shown in great detail, although it will be apparent to those skilled in the relevant art that some features or some combinations of features may not be shown for the sake of clarity.

[0022] Any reference to invention within this document is a reference to an embodiment of a family of inventions, with no single embodiment including features that are necessarily included in all embodiments, unless otherwise stated. Furthermore, although there may be references to advantages provided by some embodiments of the present invention, other embodiments may not include those same advantages, or may include different advantages. Any advantages described herein are not to be construed as limiting to any of the claims.

[0023] Specific quantities (spatial dimensions, dimensionless parameters, etc.) may be used explicitly or implicitly herein, such specific quantities are presented as examples only and are approximate values unless otherwise indicated. Discussions pertaining to specific compositions of matter, if present, are presented as examples only and do not limit the applicability of other compositions of matter, especially other compositions of matter with similar properties, unless otherwise indicated.

[0024] In some embodiments, the present invention is a method for detection and characterization of evoked potentials in muscles comprising multiple steps. Referring now to FIG. 1, in an initial step 10, epidural stimulation is applied to a patient and an electrical signal evoked from the patient's muscle in response to the epidural stimulation is detected by EMG or other suitable diagnostic process. In a first step 12, the electrical signal is temporally segmented into sequential segments. In a second step 14, image processing techniques are applied to the signal to reduce signal noise. In a third step 16, evoked potentials are detected within the electric signal. In a fourth step 18, features of the evoked potentials are extracted. In a fifth step 20, the evoked potentials are converted into a graphical form for ease of visualization by analysts. Each of these steps are described in further detail below.

[0025] 2. Materials and Methods

[0026] Participants:

[0027] In this study, five male individuals with motor complete SCI participated. Two of these participants have AIS grade B and three of them have AIS grade A. The average age of these five individuals at the time of the experiments was 29.824.458 years old and the average time since injury was 4.221.553 years.

[0028] Spinal Cord Epidural Stimulation:

[0029] A stimulation unit in combination with a chronical 16-electrode array was surgically implanted at the T11-L1 vertebral levels over the spinal-cord segments L1-S2. It was used to deliver electrical stimulation to the lumbosacral spinal cord of each SCI individual. The scES lower limbs mapping experiments started 2-3 weeks after the surgical implantation.

[0030] 2.1 EMG Data Acquisition

[0031] The supine experiments were performed with the participants relaxed in a supine position. The electrical stimulation waveform had a rectangular, biphasic shape with pulse duration of 450 sec. In each experiment a particular combination of electrode array is selected to activate, which is referred to as the stimulation configuration. A total of 12 different stimulation configurations were examined for each individual. For each configuration, the scES intensity (in volts) was increased while the frequency (in hertz) was fixed at (2 Hz). This allowed a simple evoked potential as well as a reasonable time frame in which to complete the experiments. During each intensity ramp up, the scES intensity started at a pre-threshold value and increased with 0.1 or 0.5 V increments up to 10 V. scES delivered a minimum of five stimulus pulses at each intensity level. During the performance of each experiment, the surface EMG signals were recorded from 14 leg muscles, using bipolar surface electrodes that were placed on the left (L) and right (R) soleus (SOL), medial gastrocnemius (MG), tibialis anterior (TA), vastus lateralis (VL), rectus femoris (RF), medial hamstrings (MH), and gluteus maximus (GL). The recorded EMG signals were digitized with a sampling rate of 2000 samples per second. In addition, two more surfaced electrodes were placed laterally over the paraspinal muscles in order to record the onset of the stimulus.

[0032] 2.2 Methodology

[0033] In this section, a set of five algorithms is proposed for the mapping task in order to convert the raw recorded EMG signal into its significant building blocks, which are evoked potentials induced by ES, and extract several features of these evoked potentials and visualize them in order to effectively represent the desired hidden information in the EMG records to the data analysts. Moreover, the pseudocode of each step of the framework is presented below.

[0034] 2.2.1 2-D Representation of EMG Signal

[0035] From a computational point of view, each single record of the EMG (x_k,k1) is a set of sample observations of a discrete random process (X_k,k1). The EMG signals, which were recorded from leg muscles during supine experiments, consist of evoked potentials (w_1, w_2, . . . , w_NW), which were induced by ES, and the corresponding stimulation onset timings (t_i), are marked using the paraspinal muscles signals. Utilizing this information, the whole EMG signal can be segmented into its building blocks (W set) where each segment is the time interval between two consecutive stimulation pulsations (FIG. 2A, 2B). Subsequently, the first algorithm converts the EMG signal, X_k, to a 3-D image, _k (x,y,z), by overlaying all these EMG portions and then showing their value in 3-D graphs (FIG. 2C). By using the Colormap technique, a 3-D graph can be converted to a 2-D image where the amplitude values are represented by colors (FIG. 2D). Each row in the Colormap image is one segment of the whole EMG signal (Algorithm I, below).

TABLE-US-00001 Algorithm I: 2-D representation of the raw EMG signals Find the time intervals between each two consecutive stimulation pulsations using the onset timing of each stimulation Use the time intervals to segment the EMG signals into corresponding pieces Repeat for all the muscles Save the samples of EMG segments into a m n l matrix where m is the number of stimulation pulsation, n is the maximum time interval between consecutive stimulations, and l is the number of muscles. (Optional) Visualize the m n matrix as Colormap image (use MATLAB imagesc function or equivalent) Repeat for all muscles

[0036] 2.2.2 Noise Reduction

[0037] After signal to image conversion, it is possible to use image-processing techniques for smoothing the images and consequently de-noising the signals. The 2-D generalized Gaussian Markov Random Field (GGMRF) model was applied to the constructed 2-D images so as to reduce the noise level from the EMG signal. This is an image smoothing method, which preserves continuity and removes inconsistencies in the image that in this context is caused by background noise in the original EMG signals. In this method, each pixel's value, which is the evoked potential amplitude in V, is being compared to 8-neighborhood pixel set and recalculated using Equation 2.1.

[00001]$\begin{array}{cc}{\hat{}}_s=\underset{{}_s}{\mathrm{argmin}}.Math.\left\{{.Math.{}_s-{\stackrel{~}{}}_s.Math.}^q+{}^q.Math.{}^p.Math.\underset{r{v}_s}{.Math.}.Math.b_{s,r}.Math.{.Math.{\stackrel{~}{}}_s-{}_r.Math.}^p\right\}& \mathrm{Equation}.Math..Math.2.1\end{array}$

[0038] Where .sub.s, {circumflex over ()}.sub.s and {tilde over ()}.sub.s are original pixel value, its recalculated value, and expected estimate respectively. v.sub.s is the 8-neighborhood pixel set; b.sub.s,r is the GGMRF potential; and a and are scaling factors. The parameter p[1.01, 2.0] controls the smoothing level (e.g., p=2 for smooth versus p=1.01 for relatively abrupt edges). The parameter q{1, 2} determines the Gaussian (q=2) or Laplace (q=1) prior distribution of the estimator. Our simulations were conducted with =1, =5, p=1.01, q=2, and b.sub.s,r={square root over (2)}. An example of the input and output of the GGMRF method and their corresponding muscle activation segments is demonstrated in FIG. 3 (Algorithm II, below).

TABLE-US-00002 Algorithm II: Noise reduction Initialize the GGMRF parameters , , p, q, and b.sub.s,r for Eq. 2.1 Determine the window size for the Eq. 2.1 Design an all ones m n mask matrix and set its four borders rows and columns to zero Use m n segmentation matrix, m n mask matrix, parameters and window size to substitute in Eq. 2.1 and calculate the estimated value for the pixel Repeat for all pixels Repeat for all muscles

[0039] The advantage of spatial smoothing of EMG signals, compared to the traditional signal low pass and high pass filtering techniques, is that this method offers the option of comparing each evoked potential with the previous and the next activation, and it reduces the signal irregularities that are caused by different sources of noise. Unlike the filtering methods, the image smoothing method offers a de-noised signal without any significant changes in the position and original shape of muscle activations.

[0040] 2.2.3 Activation Detection

[0041] The activation detection algorithm is designed to determine the presence or absence of scES induced evoked potentials in each segment of the EMG signals and consequently the corresponding stimulation intensity threshold. The pre-assumptions for this task are: 1) The intensity threshold Vs0, which caused the earliest evoked potential, is an unknown random value with unknown distribution; 2) The amplitude of the first visible evoked potential is also an unknown value. These are valid assumptions because of the non-stationary nature of EMG signals and the fact that no a priori information about the distribution of the changing parameters is present. It is also notable that these onset values are highly variable and alter based on the choice of stimulation configuration, frequency, muscle type, and also depending on the day-to-day and pre-post training variability. As mentioned above, scES delivers the minimum of five pulse stimulations at each stimulus voltage referred to as an event, (w.sub.i+j, 0jN.sub.i) where min(N.sub.i)=5. Therefore, it is assumed that if at least 50 percent of the EMG segments corresponding to the same stimulation intensity show epidurally evoked potentials, that intensity will be considered as the intensity threshold V.sub.s0.

[0042] The general technique implemented in this step is known as the statistically optimal decision (SOD) method. One derivation of SOD is the approximated generalized likelihood-ratio (AGLR) detector. Here, this method is modified and adapted to the present activation detection application. The modified method has three main phases. In the first phase, each segmented and de-noised portion of the signal (w_iW) is modeled by a Gaussian probability density function (pdf) and the model parameters .sub.i and .sub.i are estimated using the maximum likelihood estimation (MLE) method (p.sub.(w.sub.i)). FIG. 4B shows the histograms of the one segment of the EMG and its estimated Gaussian distribution. From the statistical point of view, the activation detection method represents a binary selection between the null hypothesis HO that says there is no significant change in the pdf p.sub..sub.o(w.sub.i) of the i.sup.th segment of the signal and the alternate hypothesis that says there is a significant change in the parameters of pdf p.sub..sub.1(w.sub.i). Therefore, in the second phase, the Gaussian model of all the EMG segments will be compared to the Gaussian model of the background noise by the log-likelihood ratio (LLR) measure. Eq. 2.2 shows the general formulation of the LLR.

[00002]$\begin{array}{cc}{s}_k=\ln \left(\frac{p_{{}_1}\left(y_{i_k}\right)}{p_{{}_0}\left(y_{i_k}\right)}\right)& \mathrm{Equation}.Math..Math.2.2\end{array}$

[0043] Where y.sub.i.sub.k is the k.sup.th sampled-value of w.sub.i segment. In order to reduce the high computational cost of this equation, it is presupposed that the occurrence of the scES induced muscle activation does not change the mean of the Gaussian pdf and the most significant changes happen in the standard deviation .sub.i of the pdf as shown in FIG. 4C. Therefore, the Eq. 2.2 is modified and simplified herein to Eq. 2.3.

[00003]$\begin{array}{cc}{s}_k=\ln \left(\frac{{}_0}{{}_i}\right)+\frac{1}{2}.Math.{\left(y_{i_k}-\right)}^2.Math.\left(\frac{1}{{}_0^2}-\frac{1}{{}_i^2}\right)& \mathrm{Equation}.Math..Math.2.3\end{array}$

[0044] The sum of all the s.sub.k values over one segment is referred to as CUSUM value Si and is calculated based on Eq. 2.4.

[00004]$\begin{array}{cc}{S}_i=\underset{k=j}{\overset{N_i}{.Math.}}.Math.s_k=\left(N_i-j+1\right).Math.\ln \left(\frac{{}_0}{{}_i}\right)+\frac{\left(N_i-j\right)}{2}.Math.\left(\frac{{}_i^2}{{}_0^2}-1\right)& \mathrm{Equation}.Math..Math.2.4\end{array}$

[0045] Eq. 2.4 is used to calculate one value for each segment of the signal that represents the highest statistical difference between that segment and the background noise (FIG. 4D).

[0046] In the third phase of the algorithm, a dynamic threshold h is calculated for each EMG signal in order to find the first segment of the signal that includes evoked potential and its corresponding stimulation intensity V.sub.s0. In comparison, traditional AGLR uses a constant threshold. Based on the experimental observations, the first event corresponding to the lowest stimulation voltage does not usually trigger any muscle activation and can be considered as baseline. Consequently, the threshold value h can be computed as the summation of maximum and standard deviation of the set of S.sub.i that belongs to the baseline (Eq. 2.5).

h=S.sub.max+.sub.S.sub.iEquation 2.5

[0047] All the steps of the activation detection method are demonstrated in FIG. 4 (Algorithm III, below).

TABLE-US-00003 Algorithm III: Activation detection using AGLR method Select the background noise as the first non-zero segment of the de-noised EMG signal Estimate the Gaussian distribution parameters .sub.o and .sub.o (use MATLAB mle function or equivalent) For each segment of the de-noised signal use mle function to estimate the Gaussian distribution parameters .sub.i and .sub.i Calculate S.sub.i using Eq. 2.4 for all segments i of the de-noised signal Select the baseline as the first event of the de-noised EMG signal Calculate S.sub.max and .sub.S.sub.i for the selected baseline to find the activation threshold h using Eq. 2.5 Compare each S.sub.i with h: If S.sub.i > h, detect the evoked potential in the i.sup.th segment and represent it with 1, else there is no evoked potential and represent it with 0 If 50% of the segments inside one event is active call the whole event active, else call it inactive Repeat for all events Find the corresponding stimulation intensity to the first active event and assign it as the voltage threshold Repeat for all muscles

[0048] 2.2.4 Feature Extraction

[0049] The objective of this step is to represent each potential evoked by scES with a set of features that offers the essential characteristics of that evoked potential. As a result, several parameters are selected for the EMG fragments that were segmented in the previous step. These parameters can be calculated automatically or observed visually using the 2-D representation of the EMG signal. Parameters such as activation latency, overall onset of muscle activations and number of peaks of the evoked potentials are observable parameters from 2-D image demonstration of the EMG signal (FIG. 5A). The calculated parameters in this framework are: Peak-to-peak value, which is the absolute value of the difference between the highest and the lowest peaks in the signal and its normalized value based on the highest peak between left and right muscle; Activation latency, the time interval between the stimulus timing and the onset of muscle activation; Time interval between minimum and maximum values, the time interval between the highest and lowest peak; Integrated EMG value, the area under the motor unit curvature after rectifying (taking the absolute value of) the EMG signal; and binary 0/1 values: an indication of the absence or presence of evoked potentials in each segment of the signal (Algorithm IV, below). The aforementioned features are illustrated in FIG. 5B.

TABLE-US-00004 Algorithm IV: Feature extraction Peak-to-peak and min-max interval: Divide each segment into 8 pieces and find the maximum A.sub.max and minimum A.sub.min amplitude values inside the first piece and their corresponding timings T.sub.max and T.sub.min. Calculate the peak-to-peak value (A.sub.pp = A.sub.max A.sub.min) and min-max interval (T.sub.mm = |T.sub.max T.sub.min|) for each segment i Repeat for all segments Take the average over all A.sub.pp and T.sub.mm values inside one event Repeat for all events Repeat for all muscles Activation latency: Divide each segment, which are detected as active by Algorithm III, into 20 pieces and take the first five pieces as baseline Apply the same method in Algorithm III to detect the onset timing of the evoked potential in each segment L.sub.i Repeat for all segments Take the average over all L.sub.i inside one event L Repeat for all events Repeat for all muscles Integrated EMG: Rectify each segment of the EMG signal Take the integral of the rectified signal I.sub.EMG Repeat for all segments Take the average over all I.sub.EMG values inside one event Repeat for all events Repeat for all muscles

[0050] 2.2.5 Visualization

[0051] The last step of the proposed framework is to represent the data processing results in informative ways to show the connection between stimulation parameters and generated results from the computer-based method for each muscle in order to make a convenient representation for the examiner to interpret the data and modify the experiments accordingly. FIG. 6 shows the transformation of 14 EMG signals (FIG. 6A) into a single Colormap image that is showing the peak-to-peak values and the binary image that represents the binary values 0/1 that are the output of the onset detection algorithm (FIG. 6B).

[0052] Not every possible stimulation intensity is typically tested in most experimental trials, particularly voltages above the minimum voltage necessary to evoke an activation potential in a muscle. For example, in one set of experiments, testing begins by applying a low voltage and ramping up with 0.1-volt increments. Once all the muscles start to show activations, increments increase to 0.5-volts, and stimulation is applied up to 10-volts, or the highest intensity that the patient feels is comfortable. This increase in increments results in some gaps between data points corresponding to overlooked stimulation voltages. Therefore in order to maintain the continuity and resolution of the Colormap images in the visualization step, the missing values need to be interpolated based on the observed values. This task has been done using cubic smoothing spline method. Assume that there is a given set of coordinates (x.sub.0, y.sub.0), (x.sub.1, y.sub.1), . . . , (x.sub.n, y.sub.n) where the values of xi are in ascending order (stimulation voltages). The objective of the cubic smoothing spline method is to link the gap between adjacent points (x.sub.i, y.sub.i), (x.sub.i+1, y.sub.i+1) using the cubic functions P.sub.i; i=0, . . . , n1; in order to piece together a curve with continuous first and second derivatives. In other words, this technique tends to preserve the values of observed data and interpolated the missing values using the piecewise curve function. The curve function also needs to satisfy the end points matching conditions. The smoothing spline method allows selection of the smoothing parameter. The smoothing parameter determines the relative weight that one would like to place on the contradictory demands of having P.sub.i be smooth (p=0), where P.sub.i is the least-squares straight line that is fitted on the data versus having P.sub.i be close to the data (p=1) so that it is variational cubic spline interpolant. In one embodiment, p=1 in order to preserve a smoothed yet accurately interpolated curves (Algorithm V, below).

TABLE-US-00005 Algorithm V: Visualization Use csaps function (cubic smoothing spline) to interpolate the missing values for A.sub.pp, T.sub.mm, L and I.sub.EMG Assign each feature values to the corresponding stimulation intensity voltage A.sub.pp Repeat for all events Repeat for all the muscles Use imagesc function to visualize each matrix as Colormap image

[0053] 3 Results

[0054] In this section, the performance of the computer-based activation detection algorithm has been evaluated by comparing the output of the algorithm with the output of manually detected evoked potentials by trained data analysts as the ground truth. The proposed method was also compared with two other existing methods to show its advantages. Finally, the run-time of all the algorithms is presented.

[0055] 3.1 Comparison with Manual Activation Detection

[0056] This comparison is based on calculating sensitivity, specificity, Dice similarity and accuracy as evaluation measures. These parameters are calculated based on true positive (TP), true negative (TN), false positive (FP), and false negative (FN) values. The mathematical equations for calculating these parameters and their values from comparing automatic vs. manual activation detection method on 700 EMG signals recorded from 14 leg muscles of five participants during supine experiments are presented in Table 1.

TABLE-US-00006 TABLE 1 Performance measurements for comparing the automated activation detection method with the manual ground truth. The five-number summary of all the calculated values are presented in this table. Total Accuracy % Total Sensitivity Total Specificity Total Dice (TP + TN)/ % % Similarity % 2TP/ (TP + FP + TN + FN) TP/(TP + FN) TN/(TN + FP) (2TP + FN + FP) Maximum 100.00 99.9999 99.9999 99.9999 Upper 100.00 99.9998 99.9989 99.9999 quartile Median 100.00 99.9997 99.9975 99.9998 Lower 97.7273 99.9996 91.6659 98.4614 quartile Minimum 94.3396 99.9994 79.5917 96.2961

[0057] This table shows the five-number summary of all the comparison measurements. These five values are maximum, median, upper and lower quartile, and minimum values of all the measurements (outliers are excluded). Moreover, the box-plots of all the performance measurement parameters for each individual are shown in FIG. 7. The mean and root mean square error (RMSE) of all the measurements (including outliers) are presented separated based on subjects in Table 2, and separated by muscles in Table 3.

TABLE-US-00007 TABLE 2 Mean and RMSE of the performance measurements for comparing the automated activation detection method with the manual ground truth as a function of subjects. Subject # 1 2 3 4 5 Sensitivity Mean 99.8226 99.4128 99.3792 99.8518 98.8426 % RMSE 1.2273 1.8114 1.9518 0.6056 8.7196 Specificity Mean 93.4503 97.1311 94.7098 95.3723 92.0207 % RMSE 11.8018 6.4772 9.6330 10.4982 11.1296 Similarity Mean 99.0072 99.1890 97.2539 99.5583 97.6363 % RMSE 2.6263 1.4585 7.5508 0.9379 8.6058 Accuracy Mean 98.5775 98.8611 97.0499 99.2641 97.6272 % RMSE 3.4938 1.9659 6.0939 1.5476 3.1851

TABLE-US-00008 TABLE 3 Mean and RMSE of the performance measurements for comparing the automated activation detection method with the manual ground truth as a function of muscles. Accuracy % (TP + TN)/ Dice Similarity % Specificity % Sensitivity % (FP + FN + TP + TN) 2TP/(2TP + FP + FN) TN/(TN + FP) TP/(TP + FN) Mean RMSE Mean RMSE Mean RMSE Mean RMSE LSOL 99.19 1.80 97.61 13.50 97.91 6.50 97.56 13.93 LMG 97.33 6.97 98.39 5.53 95.56 10.02 99.04 2.93 LTA 98.03 3.45 97.80 6.25 95.08 10.43 99.43 1.83 LMH 99.01 1.46 99.38 0.90 94.46 11.40 99.78 0.74 LVL 98.78 1.70 99.17 1.18 93.06 10.59 99.91 0.62 LRF 98.82 1.60 99.20 1.14 94.49 9.79 99.71 0.86 LGL 96.63 5.94 97.75 4.37 87.63 16.06 99.75 1.02 RSOL 98.19 4.18 97.88 7.52 98.30 4.94 99.08 2.33 RMG 97.78 6.18 98.21 6.46 96.09 9.05 99.55 1.15 RTA 98.24 2.47 98.43 2.86 94.40 8.67 99.49 1.47 RMH 98.61 2.69 99.04 1.92 94.40 9.73 99.93 0.43 RVL 98.97 1.67 99.27 1.23 96.56 6.42 99.71 1.08 RRF 98.49 2.38 98.84 1.99 95.00 7.42 99.95 0.31 RGL 97.90 2.45 98.64 1.63 91.23 12.64 99.68 1.20

[0058] 3.2 Comparison with Other Activation Detection Methods

[0059] The performance of proposed method was compared with two other methods: the Teager-Kaiser Energy Operation (TKEO) method and AGLR without GGMRF smoothing. The comparison is based on recorded EMG signals as well as simulated EMG signal. The simulated EMG signal was designed using an activation shape signal, X(k) that is convolved with a train of Dirac delta impulses .sub.n=1.sup.k0.01n (tn), where the amplitudes are linearly increasing. The additive white Gaussian noise, n(t), will then be added to the signal based on the desired signal to noise ratio (SNR) to generate the final simulated signal. This signal has 50 segments, and the first 20 segments do not include any activation, and the first 10 segments are considered as baseline. One embodiment of the pseudocode of the TKEO method implementation is presented in Algorithm VI, below.

TABLE-US-00009 Algorithm VI: Activation detection using TKEO method Segment the signal based on stimulation timings (similar to Algorithm I) Calculate the TEKO value for each sample x.sub.ijk inside the segment (TKEO.sub.ijk = x.sub.ijk.sup.2 (x.sub.i1jk .Math. x.sub.i+1jk)) Repeat for all segments j Repeat for all muscles k Take the first event after TKEO operation as baseline Calculate the baseline maximum and standard deviation and find the activation threshold h = TKEO.sub.base.sub.max + TKEO.sub.sd Calculate the maximum value for the segment TKEO.sub.max If TKEO.sub.max > h, detect the evoked potential in the current segment and represent it with 1, else there is no evoked potential and represent it with 0 Repeat for all segments If 50% of the segments inside one event is active call the whole event active, else call it inactive Repeat for all events Find the corresponding stimulation intensity to the first active event and assign it to be the activation intensity threshold Repeat for all the muscles

[0060] The results of comparing the proposed method with regular AGLR and TKEO based on the total accuracy in the recorded EMG signals from five patients are presented in Table 4.

TABLE-US-00010 TABLE 4 Comparison of methods based on five-number-summary of the total accuracy. AGLR + GGMRF AGLR TKEO Maximum 100.00 100.00 100.00 Upper 100.00 100.00 100.00 quartile Median 100.00 98.9583 98.4251 Lower 97.7273 97.4359 97.2222 quartile Minimum 94.3396 93.6508 93.1818

[0061] As indicated in Table 4, AGLR and TKEO showed lower accuracy rates compared to the AGLR+GGMRF method disclosed herein. This shows that adding the GGMRF smoothing technique to the pre-processing step increases the effectiveness of the activation detection method.

[0062] In order to show the sensitivity of the proposed framework to SNR and to compare it with other methods, the SNR of the stimulated signal was increased from 20 to 9 dB and the accuracy measurement versus SNR plotted. The results are shown in FIG. 8. It is clear from this plot that the AGLR with GGMRF significantly outperforms other methods in signals with lower SNR values.

[0063] 3.3 Calculating Total Run-Time of the Algorithms

[0064] The hardware that was used to process the EMG signals and calculate the run-time is a Dell computer, Optiplex 9020, Intel Corei7-4790 CPU @ 3.60 GHz, 16.0 GB RAM, 64-bit Operating System. The software is MATLAB R2011a. The run-time of the proposed framework is directly dependent on the length of the experiments. Table 5 presents the average run-time of each of five algorithms and the total time for 10 intensity ramp-up experiments with 5 to 7 minutes length.

TABLE-US-00011 TABLE 5 The run-time mean and RMSE of each steps of the proposed framework Converting signal to Noise Activation Feature image reduction detection extraction Visualization Total Run-time Mean 3.9342 6.6107 1.5433 0.0506 0.4090 12.7167 (s) RMSE 1.3343 1.1896 0.2480 0.0069 0.0063 2.3310

[0065] The average and standard deviation of total execution time for processing the data through all five algorithms is approximately 12.712.33 seconds, which compared to the previous manual methods that used to take up to several days, the disclosed automatic method is outstandingly time-efficient. Visualization of the extracted features of the evoked potentials in a visually perceptible format, specifically, colormap images, has been found to greatly aid experimentalists in their efforts to understand the effects of scES and its various parameters (voltage, frequency and electrode configuration) on the motor function of the paralyzed muscles in SCI patients. FIG. 9(A-C) displays a series of colormap images which allow an experimentalist to easily view how different muscles respond to changes in stimulation voltage, frequency and configuration.

[0066] 4 Discussion

[0067] Proposed herein is a novel framework consisting of five algorithms for activation detection and visualization of EMG signals recorded from multiple leg muscles of individuals with spinal cord injuries who have epidural stimulator implantation. Using this framework, a raw EMG signal is converted into a two/three dimensional image, the signal is de-noised using GGMRF image smoothing techniques, the occurrence of scES induced muscle activation is automatically detected, features of the muscle activation are extracted, and visual output is generated for interpretation. Each of these five novel algorithms has several advantages. For instance, the first algorithm that converts signals into images has several benefits such as clearly showing the latencies of all activations as well as the overall intensity onset of scES induced motor responses. It is also useful for the next step that is GGMRF image smoothing to get a de-noised signal without affecting the shape and latency of muscle activations. As it is shown in the results section, adding this step to the framework is also noticeably advantageous for the performance of the next step, which is the activation detection. In the activation detection algorithm, a statistically optimal decision method was developed by applying maximum likelihood estimator (MLE) together with comparing the probability density functions of the muscle activations to the background noise utilizing a log-likelihood ratio and calculating the dynamic activation threshold. In comparing the automatic method for activation threshold detection vs. manual detection (ground truth) on 700 EMG signals, the new automated approach developed here demonstrated an overall accuracy of 98.28% based on the errors of combined false positive and false negative data. It is important to note that the proposed method does not need any a priori information for the statistical model, which makes it a fully automatic method that uses only the current and previous EMG signal values to build the statistical model. Finally, the combination of modified SOD method and GGMRF has proven to have minimum sensitivity to the changes in the signal to noise ratio compared to the other well-known EMG onset detection methods, i.e. SOD method without smoothing and TKEO method using both simulated EMG signal and real EMG signals. Comparing the three methods on the recorded EMG signals indicates the robustness in the accuracy of the proposed activation method in the situation where no information about the noise level in the signal is known. In addition, the feature extraction and visualization steps of the framework help analysts to make accurate and quick connections between the desired EMG features and the scES parameters like intensity voltage, configuration and frequency. In conclusion, this method clearly demonstrates the advantages of implementing a set of algorithms for improving the accuracy and speed in complex EMG data analysis.

[0068] Various aspects of different embodiments of the present invention are expressed in paragraphs X1 and X2 as follows:

[0069] X1. One aspect of the present invention pertains to a method for detection and characterization of evoked potentials in muscles, comprising recording electrical signals at least one of a patient's muscles; dividing the electrical signals into sequential segments; calculating, for each segment, a highest statistical difference Si between the electrical signal in that segment and background noise; calculating a dynamic threshold h; and detecting an evoked potential in one of the patient's muscles if Si>h in at least 50% of the segments.

[0070] X2. Another aspect of the present invention pertains to a method for detection and characterization of evoked potentials in muscles, comprising delivering a plurality of epidural stimulations to a patient at a stimulation voltage, each of the plurality of epidural stimulations separated by a time interval; recording electrical signals at least one of the patient's muscles, wherein the recording is concurrent with the delivering; dividing the electrical signals into sequential segments; calculating, for each segment, a highest statistical difference S.sub.i between the electrical signal in that segment and background noise; calculating a dynamic threshold h; determining whether S.sub.i>h in at least 50% of the segments; and increasing the stimulation voltage until an evoked potential is detected; wherein the evoked potential is detected if S.sub.i>h in at least 50% of the segments.

[0071] X3. A further aspect of the present invention pertains to a method for determining the location a muscle activation potential evoked from epidural stimulation and minimum voltage required to evoke activation potential, comprising delivering a plurality of epidural stimulations to a patient at a stimulation voltage; recording electrical signals at a plurality of the patient's muscles; dividing the electrical signals into sequential segments; reducing the noise of the electrical signals; calculating, for each segment, a highest statistical difference S.sub.i between the electrical signal in that segment and background noise; calculating a dynamic threshold h; determining whether S.sub.i>h in at least 50% of the segments; and increasing the stimulation voltage until an activation potential is detected in one of the plurality of muscles, wherein the activation potential is detected if S.sub.i>h in at least 50% of the segments.

[0072] Yet other embodiments pertain to any of the previous statements X1 or X2 which are combined with one or more of the following other aspects.

[0073] Wherein the method further comprises delivering a plurality of epidural stimulations to the patient's spinal cord at a first stimulation voltage, each of the plurality of epidural stimulations separated by a time interval.

[0074] Wherein the delivering and the recording occur concurrently.

[0075] Wherein each of the sequential segments has a duration equal to the time interval.

[0076] Wherein each segment temporally overlaps with delivery of one of the plurality of epidural stimulations.

[0077] Wherein the electrical signals are electromyography signals.

[0078] Wherein the method further comprises reducing the noise of the electrical signals.

[0079] Wherein reducing the noise comprises converting each electrical signal into a 2-D image; and applying an image smoothing method to the 2-D image to reduce signal noise.

[0080] Wherein the 2-D image comprises a plurality of pixels, each pixel having a value corresponding to an amplitude of the evoked potential, and wherein applying the image smoothing method to the 2-D image comprises recalculating the value of each pixel in the 2-D image using the following equation

[00005]${\hat{}}_s=\underset{{}_s}{\mathrm{argmin}}.Math.\left\{{.Math.{}_s-{\stackrel{~}{}}_s.Math.}^q+{}^q.Math.{}^p.Math.\underset{r{v}_s}{.Math.}.Math.b_{s,r}.Math.{.Math.{\stackrel{~}{}}_s-{}_r.Math.}^p\right\}$

[0081] Wherein the 2-D image comprises a plurality of pixels, each pixel having a value corresponding to an amplitude of the evoked potential, and wherein applying the image smoothing method to the 2-D image comprises recalculating the value of each pixel in the 2-D image using a 2-D generalized Gaussian Markov Random Field model.

[0082] Wherein an event is delivery of a plurality of epidural stimulations at a stimulation voltage.

[0083] Wherein a plurality of events are applied to a patient, each event including a different stimulation voltage.

[0084] Wherein a plurality of events are applied to a patient, each event including increasing stimulation voltages.

[0085] Wherein a baseline is a first event applied to a patient

[0086] Wherein Si is determined by

[00006]$S_i=\underset{k=j}{\overset{N_i}{.Math.}}.Math.s_k=\left(N_i-j+1\right).Math.\ln \left(\frac{{}_0}{{}_i}\right)+\frac{\left(N_i-j\right)}{2}.Math.\left(\frac{{}_i^2}{{}_0^2}-1\right).$

[0087] Wherein .sub.i is an estimated standard deviation.

[0088] Wherein N.sub.i is at least 5.

[0089] Wherein h is determined by

h=S.sub.max+.sub.S.sub.i.

[0090] Wherein S.sub.max is a maximum deviation of the set S.sub.i of the baseline.

[0091] Wherein .sub.Si is a standard deviation of the set S.sub.i of the baseline.

[0092] Wherein the baseline is the lowest stimulation voltage event.

[0093] Wherein the method further comprises delivering a plurality of epidural stimulations to the patient's spinal cord at a first stimulation voltage, followed by delivering a plurality of epidural stimulations to the patient's spinal cord at a second stimulation voltage, wherein the second stimulation voltage is unequal to the first stimulation voltage.

[0094] Wherein the second stimulation voltage is higher than the first stimulation voltage.

[0095] Wherein the method further comprises extracting at least one feature of the evoked potential.

[0096] Wherein the method further comprises extracting at least one feature of the activation potential.

[0097] Wherein the feature is one of a peak-to-peak interval, a min-max interval, an activation latency, and an integrated EMG.

[0098] Wherein the method further comprises displaying the at least one feature in a visually perceptible format.

[0099] Wherein the visually perceptible format is a colormap image.

[0100] Wherein recording electrical signals at least one of the patient's muscles includes recording electrical signals at a plurality of the patient's muscles, and wherein detecting an evoked potential in one of the patient's muscles includes associating the evoked potential with one of the patient's muscles in the plurality of the patient's muscles.

[0101] The foregoing detailed description is given primarily for clearness of understanding and no unnecessary limitations are to be understood therefrom for modifications can be made by those skilled in the art upon reading this disclosure and may be made without departing from the spirit of the invention.