Quantification of Muscle Coordination and Clonic-like Motor Firing Patterns
in Spinal Cord Injured Subjects
During Stepping

B. C. Davis2, T. D. Johnson1, J. A. Beres2, K. J. Sullivan2, S. J. Harkema2
2Neurology, 1Biomathematics, University of California, Los Angeles






 
Table of Contents

::Home
::Introduction
::Methods
-Computation of Coordination
-Coordination--Examples
-Computation of Clonicity
-Clonicity--Examples
::Conclusions
::References

::Introduction
Our objective was to develop a method to quantify muscle activation patterns observed during human locomotion, specifically coordination among different muscles and clonic-like firing patterns. Previous quantitative analyses of muscle coordination have been limited to agonist/antagonist co-contraction relationships, and studies of clonus have been only qualitative. We developed a Bayesian change-point method to model step-cycle EMG data, via reversible jump Markov chain Monte Carlo simulation. This approach allows us to quantify the amount of clonicity in a given muscle, and the coordination of clonic and non-clonic activation patterns among multiple muscles, as well as make statistically valid inferences about these quantities and differences between them.

 

::Methods
In order to detect motor events (specific muscle activation patterns) in an EMG recording, we begin with a suitable model of the EMG data. Our primary assumption is that changes in muscular activity are reflected in changes in the variance of the (high-pass filtered) EMG signal. Hence, we wish to determine where the variance changes in the EMG pattern. The number and location of these change-points, as well as the EMG signal variance between change points, are thus our model parameters, represented here by the vector . In a Bayesian analysis, both model parameters and observed data (Y) are treated as random variables. Our goal is to estimate the joint posterior distribution of the model parameters, given the observed data.

Given the large number of parameters needed to model complex EMG data, the posterior distribution (the above equation) does not have a closed form solution, and we must numerically estimate it. One such estimation method is Markov chain Monte Carlo (MCMC) simulation. The Metropolis-Hastings algorithm (Hastings, 1970) is the most general algorithm for performing MCMC simulation for problems with a fixed number of parameters. However, the number of change-points, and thus the number of model parameters, changes during simulation of our EMG data. We therefore use a reversible jump Markov chain Monte Carlo (RJMCMC, [Green, 1995]) algorithm, which allows the number of parameters to vary. Thus, we use RJMCMC to estimate the posterior distribution of EMG signal variance changes, from which we can further define specific muscle activation patterns.

Our method can be summarized by these five steps:

  1. RJMCMC to select the number and location of change-points along with the inter-change-point variances and the variance of the background EMG signal.
  2. Calculation of the SNR based on the estimated posterior distributions of the variances from the EMG signal and background.
  3. Determination of EMG activity/inactivity based on the SNR and appropriate thresholds.
  4. Computation of quantities of interest-coordination between muscles and clonicity.
  5. Statistical inference about these quantities.

RJMCMC
At each iteration of the simulation we propose, at random, one of the following moves—each with a probability of 1/3.

  1. addition of a change-point
  2. deletion of a change-point
  3. position change of a change-point

The RJMCMC simulation is run for 300,000 iterations. The initial 200,000 iterations are considered a burn-in until convergence of the Markov chain is reached. By convergence, we mean that the Markov chain is now sampling from the posterior distribution. During development, most cases had converged by 150,000 iterations. We run the simulation for an extra 50,000 iterations, to assure convergence. These first 200,000 burn-in iterations are discarded.

We then save every 50th iteration, so estimates of the posterior distribution are based on a sample of 2,000 simulated draws from the joint posterior distribution of the parameters.

These Markov chain samples, however, tend to be autocorrelated (correlation between successive values), which potentially increases the bias of distribution estimates. To reduce the autocorrelation we save every nth iteration, where n is chosen to minimize autocorrelation. For muscle coordination, we found that saving every 50th iteration reduced the autocorrelation to an acceptable level.

The mean and each variance of the inter-change-point, stationary Gaussian distributions are updated at each iteration after convergence is reached, before which they are updated on every 25th iteration.


click for full size

Figure 1 represents 4 typical iterations. Acceptance or rejection of each proposed move is based on comparing a uniform random variate, u ~ U[0,1], to a ratio of the proposal density ratio over the posterior density ratio (see the Addendum).
Panel A) A deletion is proposed. The red change-point was randomly selected as the candidate for deletion. In the present case the change-point was deleted.
Panel B) A location change is proposed to the green change-point. We randomly select a new location from a truncated normal distribution centered on the change-point between its two adjacent neighbors (within the green shaded background). The new location (11.7) was accepted and the change-point moved.
Panel C) Insertion of a change-point is proposed, at the position depicted in blue. Here the change-point was rejected.
Panel D) Another change-point is proposed at a new location, depicted in blue.


click for full size
Signal-to-noise Ratio Calculation
Each RJMCMC iteration is a sample from the posterior distribution of the parameters. The SNR for each sample is calculated, at each time point t, by the following:

Figure 2, panels D, E, and F show the mean SNR, along with the 0.025 and 0.975 quantiles, for the Soleus, MG and TA from a single step obtained from a non-disabled volunteer. The estimated mean SNR at time t is the ergotic average SNR--that is averaged over all saved iterations.

EMG Activity/Inactivity
Three thresholds are used to convert the SNR into estimated posterior probabilities of EMG activity. Periods of muscle activation are defined by:

  • SNR is greater than 2.5
  • Minimum activity duration of 32 msec
  • Minimum inactivity duration of 38 msec

Computation of Coordination
We then define coordination between two muscles in terms of the proportion of time they are simultaneously active and simultaneously inactive according to the following formula:

Figure 3 demonstrates the computation of coordination between the SOL and MG and between the SOL and TA. “On” represents activity and “Off” represents inactivity. Time segments marked by:
A= both muscles simultaneously active
B= either is active.
C=both muscles are simultaneously inactive
D= either is inactive


click for full size

Statistical Testing
We then compute the probability that the coordination between any two of the three pairs of muscles is equal.

  • Compute difference in coordination for the 2 pairs for each saved iteration, resulting in a distribution of coordination difference values.
  • If 0 lies between the 0.025 and the 0.975 percentiles, the two coordination values are similar, while if 0 lies in one of the distribution tails, the coordination between the two pairs is (significantly) different.


click for full size

Figure 4, panel A shows the estimated posterior densities of coordination between the three major lower leg muscles . These estimates are smoothed histograms based on the 2,000 samples.
Figure 4, panels B, C, and D demonstrate the statistical test for differences in coordination between the 3 muscle pairs. If 0 lies under a red portion of the estimated density, we say the two pairs differ in coordination.

Coordination--Examples
Figures 5 through 9 display the EMG of the SOL, MG and TA, the SNR, and posterior probability of EMG activity from two steps from ND and SCI subjects. A minimum of at least eight steps each were used in the analysis.
In each figure, panels D, E and F show the mean SNR, the 0.025 and 0.975 quantiles. The estimated posterior probability of EMG activity is shown in panels G, H and I.

for all graphs: click for full size

Figure 10 (above). Posterior density esitmates of coordination between the three lower leg muscles for each subject.

Table 1. Means and standard deviations of the estimated posterior coordination distributions for each of the five subjects studied.

Table 2. The probability of coordination differences between the muscles within each subject (red indicates a sufficiently small probability of equality, <0.05).

Computation of Clonicity
We compute clonicity as the ratio of clonic-like EMG activity to all active EMG in the recording period.
This definition can take values in [0,1]. An absence of clonic-like activity will result in a value of zero, and purely clonic-like activity in a value of one.

Definition of Clonus
We define clonic EMG as follows:

Average Burst/Silent Duration and Frequency
Using our Bayesian methodology, we can also esimate the posterior distributions of average burst duration (ABD), average silent duration (ASD), and frequency.

Clonicity--Examples

Figure 11 shows the EMG (panels A, C, and E) and SNR (panels B, D, and E) of the three lower leg muscles during 3 seconds of clinical clonus. Also shown are the upper and lower 0.025 quantiles of the SNR. Of interest is the coordination of clonic activity of the TA with the Soleus and MG.

Figure 12 shows clinical clonus as well as clonic-like activity during stepping from the Soleus of a patient with a SCI.

Figure 13 (above) shows the estimated posterior densities of Soleus clonicity for the four subjects. The clonicity density from the non-disabled volunteer is shown in panel A. Panels B-F show clonicity densities in increasing order of clonicity.

Table 3. Mean values and standard deviations of the estimated posterior distributions of clonicity, ABD, ASD and frequency during clinical clonus.

Table 4. The probabilities that these quantities are equal among the three lower leg muscles during clinical clonus for SCI-C1.

Figure 14 shows the estimated posterior densities of the difference of clonicity for three selected cases (corresponding to Figure 12). Panel A shows that clonicity is different for the two different speeds from SCI-C1.

 

::Conclusion
This method demonstrates the utility of a Bayesian approach to analysis of EMG data, in its ability to make statistically valid estimates of the probability of muscle activity/inactivity, and thus specific muscle activation patterns. Probablistic estimates of temporal activation patterns (on and off-set times) not only allow statistical comparisons of relevant quantities, but also reduce the amount of unknown error associated with such measurements.

We further show the ability of modern Markov chain Monte Carlo simulation to estimate posterior distributions of parameters in complex models of EMG data.

The application demonstrated here could be used to address many investigational goals. These quantitative measurements could serve as meaningful clinical measures of specific neuropathologies, normal locomotor muscle activation patterns, and gait rehabilitative progress and strategies.

 

::References
Green, P. J. (1995). Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika 82, 711-732.

Hastings, W. K. (1970). Monte Carlo sampling methods using Markov chains and their applications. Biometrika 57, 97-109.

Johnson, T. D., Elashoff, R. M., and Harkema, S. J. (submitted). Bayesian change-point analysis of muscular coordination from electromyographic data. Biometrics.
    URL: http://www.biomath.medsch.ucla.edu/faculty/tjohnson/