Estimation of a single motor unit's threshold and activation range, a study on patients with muscular disorders.

BACKGROUND: In clinical neurophysiology, threshold tracking studies are used to evaluate the functionality of a muscle through studying the functionality of its motor units (MUs) that govern the muscle. The functionality of an MU can be quantified by estimation of its excitability properties via MU's stimulus-response curve. In this study, we aim to develop a model-based approach to estimate MU's threshold mean and its activation range as indications of MU's excitability. This is a different approach from routine strategies in neurophysiology, which are mostly subjective.
METHODS: To assess the excitability of a single MU, needle electromyography examination was used to obtain the axonal activity of that MU. To improve estimation, the examination was repeated several times on individuals. Replication of experiment introduces serial correlation between observations. We account for this correlation by using a mixed-effects model. We investigate the appropriateness of classical logistic mixed-effects model and its Bayesian formulation for estimation purpose.
RESULTS: Both classical and Bayesian models can obtain a reliable estimation of MU's threshold. However, we found Bayesian approach to provide a better estimate of MU's activation range. Moreover, if data contain outliers both classical and Bayesian methods are vulnerable to some extent.
CONCLUSIONS: Compared to the classical approach, Bayesian method is more flexible in dealing with overdispersion and provides more robust estimation of MU's parameters.

INTRODUCTION

Needle electromyography (NEMG) examination is a routine clinical technique to examine the functionality of muscles. The contraction of muscles releases a microvolt level electrical signals that can be captured through a needle electrode inserted into the target muscle of the subject's tissue. By quantification of these bioelectrical signals, NEMG technique enables clinicians to study the nerve and muscular disorders.

A key concept in clinical EMG is the so-called motor units (MUs). A MU is formed by an alpha motor neuron in the spinal cord, its long protruding axon, and the tens to thousands of muscle fibers that this neuron innervates.[12] In patients with some extend of muscular disability, the muscle activity is limited due to nonfunctionality of a group of MUs. In testing muscle activities, the aim is to gather information from MUs’ activity using conventional methods to understand the extent of loss (for future therapies and evaluation of disease progression). Of particular interest is to obtain a full description of MU's functionality in this evaluation.[3]

In the laboratory, in the event of muscular disorder, the interest is to isolate a single MU to study its excitability. The excitability of a MU is reflected in the intensity of the electrical stimulus that is necessary to evoke an action potential (electrical signal) in MU's axon. Therefore, in the so-called threshold tracking studies, a full description of a MU is obtained by observing MU's responses to a range of stimulus intensities over which it displays a stochastic behavior (MU activation range).[4] The activation range is spanned around the unknown MU's threshold so that it starts with a low stimulus intensity that is not large enough in the magnitude to innervate the nerve that supplies the MU, and ends with a stimulus intensity at which the nerve is always activated. MU's threshold is usually defined as the stimulus intensity that elicits a response to 50% of the stimuli.[5] Generally, a MU fires (response) if the applied stimulus exceeds its threshold, otherwise it remains off (no response).[6] In this sense, the MU responses are binary observations.

As an indication of the MU's excitability behavior, a threshold tracking study leads to MU's stimulus-response curve (a sigmoid curve). Figure 1 displays an example of stimulus-response curve obtained from one NEMG experimentation. The horizontal scale represents the stimulus intensity values, varied from 59.9 mA to 69 mA. For each stimulus intensity, the vertical scale represents the number of successful MU's responses (firing) out of 50 attempts.

Figure 1

A typical stimulus-response curve showing the general excitability behavior of a single motor unit. The motor unit (MU) starts with no response at lower stimulus intensities and ends with always active at high stimulus intensities. The mean threshold is located somewhere close to the middle of x-axis where MU response to an applied stimulus is totally unpredictable. Activation range is the range of stimulus intensities on x-axis

The s-shape curve in Figure 1 is known as MU's stimulus-response curve. In practice, MU's threshold (mean of the curve) and its activation range (slope of the curve) are calculated from the stimulus-response curve and interpreted as MU's excitability properties. Thus, one important question is how one can obtain efficient estimates of these parameters. For estimation of MU's threshold, various strategies are available.[5] However, the judgments are mainly subjective, and statistical methods are seldom used.[6]

In this paper, our goal is to investigate how statistical methods can be merely employed for estimation purpose of both MU's threshold and activation range parameters. We use a real NEMG data to apply our proposed models. NEMG data are noisy thus a great care is needed to model the heterogeneity of data properly. This issue is discussed in Section 3 and our strategy to tackle this problem is introduced. Our models include classical approaches, which are readily available from standard statistical packages, and a more sophisticated approach that it's implementing requires experience and deeper statistical skills (Bayesian models).

The paper is organized as follows: NEMG experiment is introduced in Section 2. At Section 3, our strategy for estimating MU's threshold and activation range parameters is introduced. Results are given in Section 4. The paper ends with a discussion in Section 5.

A brief introduction to needle electromyography experiment

The NEMG is a routine technique in diagnostics of muscular disorders. In this technique, experimenter inserts a needle electrode into the target muscle of the subject's tissue to acquire micro-wavelet electrical signals. In the NEMG experiment we used in this study, the needle was manipulated to isolate a single MU. A potential was considered acceptable when it demonstrated marked “all-or-none” behavior and could be clearly separated out from other potentials throughout the recording session. When visual inspection indicated that such a stable potential was obtained, stimulus intensity was increased until no probabilistic inactivity occurred during 50 stimuli. From that level, which marks the end of the MU's activation curve, stimulus value was decreased in the smallest steps allowed by the equipment until the muscle fiber responded to none of 50 stimuli. At each step, the number of MU's responses to the 50 stimuli was recorded. The experiment was repeated several times at couple of stimulus voltages.

Table 1 shows an example dataset of this experiment with six replications on a subject suffering from amyotrophic lateral sclerosis. At each replication, the MU's responses have been observed at a few distinct stimulus intensities (s). At each stimulus intensity, the observation (y), is the number of MU's firings out of 50 stimuli.

Table 1

The results of parameter estimations under Bayesian and non-Bayesian (classical) logit mixed-effects model using data from Table 2

Table 2

A MU activity recorded by needle EMG experiment for an individual at six replications

In a typical NEMG experiment such as Table 1, MU's responses can be treated as binomially distributed observations.

Procedure

Motor unit's threshold and activation range parameters

In this section, we develop our statistical model to estimate MU's parameters. Since MU's observations were binary identities, NEMG observations modeled under a logistic regression function. However, there was one major issue with experimentation to tackle; NEMG experiment was carried out over a long period of time (several minutes) where during this period the needle inserted to evoke MU moves slightly altering the stimulus voltage and the firing rate. This introduces a new source of variation into the data, the so-called overdispersion or extra-binomial variation. Because of the problem of overdispersion, to improve estimates experiment was repeated on an individual. In statistics overdispersion on binomially distributed data can be handled using a logistic mixed-effects model.[7]

To introduce the notation, let's denote y as the number of MU's responses for a given stimulus intensity s at replication i = 1, 2…r and time t = 1,2…n so that y❘s ~ Bernoulli(p. Consider the following model;

log it(p❘.) = b (s – m) + ɛ      (1)

Where p is the probability the MU to fire at time t and replication i, m and b are threshold and slope (activation range) parameters of the MU in replication i respectively. Model 1 assumes that both threshold and slope parameters are subject to change when we move between replicated NEMG experiments. That is, b and m are random effects distributed normally with mean 0 and variance σ2 and σ2 respectively. The error terms ɛ ˜ N(0, σ2) are added to the model to account for overdispersion.

The main criticism of the classical model 1, however, is that it underestimates the uncertainty around the estimates. Bayesian approach solves this inadequacy by the formulation of the beliefs available prior to experimentation. This is the main motivation to consider Bayesian methods as an alternative choice with the hope to improve parameter estimations. For the sake of simplicity, we skip introducing mathematical Bayesian formulation for model 1. Readers can refer to[8] for details of maths and other justifications of Bayesian approach.

RESULTS

Weakly-informative priors were used for Bayesian models. Table 1 displays the results of parameter estimations of two models. Note that σ2 and σ2 are indications of variability between MU's threshold and slope parameters at six replicates and σ2 can be used as an indication of overdispersion. The R statistical package[9] was used to fit model 1 under classical model and WinBUGS software[10] for Bayesian paradigm. The result of MU's threshold and its activation range parameters for both models are given by Table 1.

Note that classical and Bayesian estimates are similar but have different interpretations. For instance, under model 1, the value 67.54 is the maximum likelihood estimation of the threshold and 0.38 is its standard deviation whereas under model 2, is the posterior mean of the threshold parameter distribution and 1.08 is the posterior standard errors.

Out of six NEMG replications presented in Table 2, for unknown reasons replication 6 was not completed. To address the effect of this unexpected disturbance on parameter estimations by two statistical methods, this replication was removed from the analysis. The results are given by Table 3.

Table 3

The results of parameter estimations under classical and Bayesian logit mixed-effects models after excluding replication 6 from data in Table 2

DISCUSSION AND CONCLUSION

From the results presented in Table 1, the immediate conclusion is the similarity between classical logit mixed model and its Bayesian estimations. However, a mild difference can be observed in estimation of the activation range and variance components between two models. The larger posterior estimation (or maximum likelihood estimator) of variance between thresholds () can be due to the last replication with threshold mean estimation of 71.52, markedly different from others. Activation range of this replication (and hence its threshold) is significantly different from others [Table 2]. Note that NEMG experiment for this replication is not complete and has been terminated for unknown reason/s. In general, under Bayesian approach, the estimation of variance components ( and ) is larger than that of classical method but estimation of , an indication of overdispersion, is smaller. This means Bayesian model can handle overdispersion in the NEMG data more efficiently. However, both approaches suffer to some extend from overdispersion (=1 means no overdispersion).

As expected, after discarding replication 6 [Table 3], the variation among thresholds was decreased significantly for both models. Overdispersion was improved under Bayesian model but no improvement achieved using classical model. Furthermore, note that for classical logit model, σb ≈ 0. As Punheiro and Bates[11] noted, this variance should not be interpreted as no variation between activation range in different replications, it might be an indication of small variation between them that is not large enough to be captured by classical model. In other words, the variations have to be large in order to be estimated by the classical model. Otherwise it fails to estimate this variation appropriately.

Generally speaking, using a real NEMG experiment, the main conclusion is that both approaches can achieve a robust estimation of threshold parameter, but they dealt with estimation of slope parameters (activation range) in different manners. Bayesian approach appeared to provide better estimate of slope parameter, mainly due to better handling of overdispersion. However, when outlier presented in data (such as replication 6), both classical and Bayesian method were vulnerable though latter model showed better performance. It should be noted to avoid exaggerating our finding in better performance of Bayesian approach for overdispersed data, a comprehensive simulation study is recommended to demonstrate the reliability of both approaches in a broad spectrum of data. This can be a topic for future research.