IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING, VOL. 23, NO. 2, MARCH 2015

297

Inter-Discharge Interval Distribution of Motor Unit Firing Patterns With Detection Errors Javier Navallas, Javier Rodriguez-Falces, and Armando Malanda

Abstract—Inter-discharge interval (IDI) distribution analysis of motor unit firing patterns is a valuable tool in EMG decomposition and analysis. However, the firing pattern obtained by EMG decomposition may have detection errors: false positives (incorrectly classified firings) and false negatives (missed firings). In this paper, the mathematical derivation of an IDI distribution model that accommodates false positives and false negatives of the detection process is presented. An approximation of the general model to adapt to specific EMG decomposition conditions is also presented. To illustrate the usefulness of the model, the obtained distribution is used to derive the maximum likelihood estimates of the statistics of motor unit firing patterns, the IDI mean and standard deviation, and estimates of the false negative and false positive ratios. Results obtained from simulation experiments and tests with real motor unit firing patterns show an enhanced estimation performance when compared to previously available algorithms. Goodness-of-fit tests applied to estimations for real data corrupted with false positives showed that the model-driven estimations fitted the uncorrupted data better than EFE estimations: 82% versus 52% not rejectable, respectively, when false positives were about 10% of IDIs. With about 5% false positives, the not rejectable estimations were 85% versus 70%. Index Terms—Electromyography (EMG), inter-discharge interval (IDI), motor unit firing pattern, motor unit potential train.

I. INTRODUCTION

M

OTOR unit (MU) firing pattern is defined as the set of firing instants of a given MU. Inter-discharge intervals (IDIs) can be calculated as the time difference between consecutive discharges within an MU firing pattern. MU firing patterns can be obtained by applying EMG decomposition techniques that identify different motor unit potential (MUP) waveforms and determine the firing instants of each MUP train from an EMG signal [1]. In automatic EMG decomposition algorithms, information about the IDI statistics [IDI mean and standard deviation (SD)] of MU firing patterns is widely employed in order to both refine the decomposition [2]–[4] and to validate the decomposition [5]–[8]. IDI statistics are also widely employed in the study of muscle physiology [9], [10]. Two different kinds of detection errors may appear when determining the firing instants within a MU firing pattern during the EMG decomposition process: false positives, which are discharges in the detected firing pattern that do not correspond to

Manuscript received May 05, 2014; reised October 05, 2014. Date of publication October 16, 2014; date of current version March 05, 2015. The authors are with the Department of Electrical and Electronic Engineering, Public University of Navarra, Pamplona, Navarra 31006, Spain (e-mail: javier. [email protected]). Digital Object Identifier 10.1109/TNSRE.2014.2363133

true discharges from the MU; and false negatives, which are discharges from the MU that are missed and therefore not included in the detected MU firing pattern. The impact of these two kinds of errors on the resulting IDIs is also distinct. In the case of false positives, calculated IDIs will be shorter than the actual IDIs, due to the presence of artefact intermediate firings. In the case of false negatives, calculated IDIs will be longer, because they will be the summation of the two or more actual IDIs. Mathematically, a MU firing pattern can be regarded as a realization of a point process. In the absence of detection errors, during short duration ( 10 ) constant force isometric contractions [9], the firing pattern can be modelled as a renewal point process [10]. Within this model, all the IDIs within the MU firing pattern are independent and identically distributed following a certain probability density function (PDF), and the process is completely determined by an adequate statistical description of the IDI PDF [10]–[12]. Experimental observations support the hypothesis that the IDI distribution is Gaussian or Gamma with low skewness [12]. Such IDI distribution models as these, which ignore detection errors, have been used, through their derived statistics (IDI mean and SD), to refine and validate EMG decomposition methods [4]–[6], [8]. False negatives of the detection process can be incorporated in the IDI PDF model by including a certain probability of not detecting individual firings during the decomposition process [13]. Essentially, if some firings are missed, some of the detected IDIs become longer: they represent the summation of several actual IDIs [Fig. 1(b)]. The probability of consecutively missing more than one firing drops as the number of missing firings increases. Hence, the IDI distribution can be modelled as the summation of Gaussian distributions with increasing mean and SD but with lower and lower probability. This model has been successfully employed to calculate the MU firing pattern certainty in EMG decomposition [8] and to estimate IDI statistics of MU firing patterns [14]. To the best of our knowledge, there is no published description of an IDI distribution model that takes into account false positives of the detection process. Given that classification errors are common to all automatic decomposition algorithms [15]–[19], the lack of such a model is a gap in the set of available analysis tools. The aim of this paper is to present a mathematical model for the IDI distribution when both false positives and false negatives are considered. In order to demonstrate the applicability of the model, the distribution model is used to estimate the IDI statistics (IDI mean and SD) as well as the probability of missing firings (false negatives) and the probability of including misclassified firings (false positives), using maximum likelihood estimation (MLE). The

1534-4320 © 2014 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

298

IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING, VOL. 23, NO. 2, MARCH 2015

Fig. 1. (a) Original MU firing pattern given by the set of firings , with ; (b) Thinned MU normally distributed IDIs, firing pattern after individually and independently detecting each firing with probability , and the corresponding set of IDI observations (in the example two firings are missed, generating two increased IDIs); (c) MU firing pattern with additional firings due to classification errors, occurring in proportion during the EMG decomposition. Each additional firing splits an IDI into two , there is new erroneous IDIs. Three different cases are depicted: when no false positive firing (F-firing) between the two adjacent true positive firings IDI is measured; when , F-firings lie in the (T-firings), and the IDIs (leading and closing IDIs in the interval), and TT interval, and 2 IDIs are measured.

MLE performance when estimating the IDI mean and SD is compared to that of error-filtered estimation (EFE) [20], maximum likelihood estimation [14] using an IDI PDF model that does not take into account false positives [13] ( ), and trimmed median estimation (MnT) [4]. Comparisons are made both in simulation experiments and with real MU firing patterns. In addition, the estimator performance when estimating false negative and positive ratios is compared with two parameters previously developed to estimate these quantities: the identification rate and the inconsistent percent of IDIs [5]. II. DERIVATION OF THE MODEL A. Derivation of IDI PDF Model In constant force isometric contractions and if duration is short enough, firings from a MU firing pattern , can be modeled as a renewal point process [10]. The IDIs, , calculated as the temporal differences between consecutive firings, can then be modeled as independent and identically distributed random variables following a distribution [Fig. 1(a)]: (1) We refer to such an IDI as a true-true (TT) IDI, because it is the interval between two true positive firings (T-firings), i.e., two firings from the MU firing pattern under analysis. The detection process of the individual firings of a MU firing pattern, e.g., EMG decomposition, may miss some of the firings (false negatives of the detection process) with probability , where is the firing detection probability [Fig. 1(b)]. Assuming that each firing is independently detected with probability , the probability of consecutively missing firings and then detecting the next firing is . An IDI measured between two firings with missing firings between them is a new random variable distributed according to the sum of the individual IDI random variables. If we model the TT-IDIs with a Gaussian distribution, , where is the mean IDI and is the SD

Fig. 2. Steps taken in the calculation of the observed IDI PDF: 1) gives the probability of finding F-firings within a TT-IDI, and, hence, the probability of splitting the TT-IDI into several TF- and FF-IDIs. In the text it is demonstrated that TF- and FF-IDIs follow the same distribution, and so the term XF-IDIs will be used instead to refer to both; 2) Having calculated the and distributions, evaluating for the correct in each case gives the IDI PDF for this particular value of ; 3) Finally, the observed IDI PDF can be composed by weighting each TT- and XF-IDI PDF distribution by its contribution, i.e., by the number of IDIs contributed for each value of , and normalizing the resulting function by the overall number of IDIs, which can . be directly calculated from

of the IDIs, the detected IDIs when considering missing firings is distributed according to [13] (2) Because of the detection process, some false positive firings (F-firings), which do not belong to the actual MU firing pattern, may be included in the detected train [Fig. 1(c)]. Given that F-firings come from other MU firing patterns, and neglecting synchronization between MUs, we can assume that the two point processes are independent and, consequently, expect that F-firings appear uniformly distributed along the time dimension. Therefore, the train of erroneously detected firings can be modeled as the events within a Poisson process with intensity , where (3) is the average T-firing appearance time, and and where are related in such a way that the ratio of F-firings to T-firings equals . Under this model, when F-firings are considered in isolation of T-firings, intervals between two consecutive F-firings are exponentially distributed (4) After the detection process, the observed set of firings is the combination of T- and F-firings. The observed IDIs include [Fig. 1(c)] TT-IDIs, between two true firings; TF-IDIs, between a true and a false firing (or a false and a true firing); and FF-IDIs, between two false firings. One direct consequence of the combination of both firing processes is that the observed TT- and FF-IDIs are no longer distributed according to (2) and (4), respectively. This is a result of the fact that the union of the two point processes redefines the IDIs and consequently modifies the probability of occurrence of the different IDI values.

NAVALLAS et al.: INTER-DISCHARGE INTERVAL DISTRIBUTION OF MOTOR UNIT FIRING PATTERNS WITH DETECTION ERRORS

299

Fig. 3. (a) Joint distribution, , evaluated for the first four values. Note that as increases, the distribution mean increases and the area decreases. The former effect means that higher values are more likely in larger TT intervals than in smaller TT intervals, while the latter effect means that higher values are less likely to occur when compared to lower values; (b), (c) Two examples of the IDI PDF decomposition into its components, and the effect of the classification ; (c) . In both cases: , , and . In the graphs, the model is decomposed into the contributions of , error, : (b) when , and the independent summation of the and for all . IDI PDFs obtained by applying (13) (solid line) are compared to the histograms obtained from 100 000 IDIs simulation experiments (gray area), showing an exact fit. Note that only the two first components are relevant when is small.

In order to calculate the observed IDI PDF, the following steps must be taken (Fig. 2): 1) calculate the probability distribution of observing F-firings within a TT interval; 2) calculate the PDFs of the TT-, TF- and FF-IDIs for each value, i.e., , , and ; 3) sum the contributions of the TT-, TF- and FF-IDIs to the observed IDI PDF for all the possible values of . 1) Distribution of Number of F-Firings Within TT Intervals, : First, we obtain the probability of finding F-firings between two adjacent T-firings. The number of F-firings between two adjacent T-firings separated by is distributed as the count of the Poisson process, i.e., it follows a Poisson distribution given by (5) combines the number of F-firings per where the ratio second, , and the interval between the two T-firings during which the count is carried out, . This interval is not constant, but follows a distribution, , as in (1). Hence, the joint density for and , i.e., the joint distribution of finding F-firings within two adjacent T-firings separated by , is (6) The interpretation of this joint distribution is straightforward when depicted [Fig. 3(a)]: events with higher values tend to appear in longer TT intervals, and vice versa; also events with higher values have lower probability of occurrence (lower area when evaluating the joint PDF at a given ) than events with lower . The overall probability of a TT interval having exactly F-firings can be calculated as (7) , 2) Distributions of TT-, TF- and FF-IDIS, , and : Second, we calculate the distribution of the observed IDIs of each type (i.e., TT, TF and FF) for each value of , i.e., the joint distributions.

The observed TT-IDIs are those with no F-firings within their ). The other cases ( ) are split into combiinterval ( nations of TF- and FF-IDIs [Fig. 1(c)]. Hence, the joint PDF of the observed TT-IDIs is given by and may be directly calculated from (6). For the observed TF- and FF-IDIs, the distributions depend on the number of F-firings in each TT interval, provided there is at least one F-firing within a TT interval ( ). The F-firings are the events of a Poisson process, and so they are uniformly distributed in time. Hence, for a given , firings lie uniformly distributed in the TT interval, . The time from the leading T-firing to the first F-firing is the time of the first ordered statistic, , of a sample of size drawn from a uniform distribution defined in the interval (and subsequently scaled to cover the interval), which is distributed as a [21] (correspondingly scaled). In view of the symmetry of the problem, the same distribution applies to the TF-IDI measured between the last F-firing of the interval and the closing T-firing. For FF-IDIs, each of the IDIs are the spacings, i.e., the differences between consecutive order statistics. It can be proved that, when coming from a sample from a uniform distribution, all the spacings are equally distributed [22]

(8) where is the th spacing, , and stands for equally distributed. Hence, the FF-IDIs follow exactly the same distribution as the TF-IDIs for a given . In the following, we refer to the XF-IDIs to refer indistinctly to both kind of IDIs. For a given and (9) is introduced to compensate for the where the factor linear scaling of the time variable converting a uniform sample in the interval into the interval.

300

IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING, VOL. 23, NO. 2, MARCH 2015

Given the distributed for each tribution can be calculated for any

as in (6), the joint disas

(10) : To calculate the observed 3) Observed IDI PDF, IDI PDF, the distribution of each IDI type for each value must be weighted by its relative contribution to the complete set of IDI observations. On one hand, , as given by (7), gives the probability of occurrence of F-firings within a TT interval. On the other hand, we know that IDIs are measured in each of the cases (Fig. 2). If , only one TT-IDI is observed. If , two TF-IDIs are observed. If , two TF-IDIs and one FF-IDI are observed. Generally, if , 2 TF-IDIs and FF-IDIs, hence XF-IDIs, are observed. Hence, the observed IDI PDF model, , can be obtained by adding the joint PDFs of the TT- and XF-IDIs weighted by the number of IDIs in each case, and normalizing the whole distribution according to overall number of IDIs as dictated by the distribution

Two aspects of the model must be highlighted. First, the derivation of the IDI PDF model is independent of the distribution chosen. Second, the final IDI PDF model incorporates false positives directly, while false negatives must be previously incorporated in the model. In Fig. 3(b) and (c), two examples of the accuracy of the IDI PDF model are shown. In each case, a 100 000 IDIs sample was generated for the given parameters. The histograms obtained from the simulations match the predicted model IDI PDFs exactly, demonstrating the correctness of the model. B. Approximation of IDI PDF The exact IDI PDF model given by (13)–(16) is rather too complex for use in an intensive computational environment. However, several approximations can be made without a noticeable loss of precision. First of all, in a real EMG decomposition environment, detection errors leading to F-firings are not common, due to the fact that the MUP shapes of different MU firing patterns are different, and shape information is used in the classification process. Ratios of classification errors to total number of detected firings are usually below 10% [20]. Hence we can safely assume that the mean interval between two F-firings is much higher than the mean interval between two T-firings, or (17)

(11)

Expanding the denominator, and taking into account that is a discrete distribution (12) where is the mean number of F-firings within the TT intervals and, hence, corresponds to the ratio of F-firings to T-firings, . The final IDI PDF model including both false negatives and false positives of the detection process is

Under this assumption, the distribution in (7), giving the number of F-firings between two T-firings, will be highly (measkewed, with most of the probability falling in suring a single TT-IDI) and (measuring two TF-IDIs). Hence, the occurrence of FF-IDIs will be almost negligible and, similarly, so will the number of TF-IDIs coming from more than two F-firings within the TT interval [Fig. 3(c)]. Hence, the distribution of the XF-IDIs can be approximated by the component. This, in turn, means that, as long as F-firings occur with frequency less than 10% of the total number of IDIs, it is sufficient to retain only the first component of the F-firings contribution to the overall PDF. And, given (13), the overall distribution simplifies to

(13) where, without loss of generality, the term culated as (Appendix A)

can be cal-

(18) Under a Gaussian model for the TT-IDIs PDF, as given in (2), term equals (Appendix C) the

(14) and the terms (Appendix B)

when

can be calculated as

(19) And

is (Appendix D) (20)

(15) where (16)

where (21)

NAVALLAS et al.: INTER-DISCHARGE INTERVAL DISTRIBUTION OF MOTOR UNIT FIRING PATTERNS WITH DETECTION ERRORS

Note (Appendix D) that the calculation of should include additional factors which are not included in (20). This can be understood from a practical point of view: the value of at should be equal to the weighted contribution of the free component of given in (4). This is so because there is a negligible probability of a T-firing splitting such a small FF-IDI (if ). That is (22) Comparison of this value with the contribution of shows that avoidance of the correction leads to the exact value in (22) (Appendix D). That is, by eliminating the factor in , the value of its contribution to the overall distribution at is fitted such that it is exactly the value that we would get by considering all the contributions without the approximation. This is important because the most significant contribution of the XF-IDIs to the overall PDF occurs with excessively small values, where no TT-IDI should be present ( ). The approximation in (18)–(21) can be considered to be accurate for , independently of the value of . This is the case because, even when is low and is applied through many terms in (2), each of these terms includes its own correction term. When , , and this is sufficient for (17) to be considered valid. A final approximation involves limiting the number of summated terms in the calculation of in (2). Theoretically the number of terms is infinite, but in practice it is generally sufficient to limit the calculation to [14], hence summations in (19) and (20) may also be limited to this value.

which, after manipulating (23), leads to

(25) , effectively no truncation is performed, Note that, if and (25) is the same as (18). In the following we will still refer to , although, if a is used in the automatic EMG decomposition algorithm employed to derive the IDI data set, will be used. D. Maximum Likelihood Estimation of IDI Statistics The derived IDI PDF function, , depends on four parameters: , the IDI mean of the original MU firing pattern; , the corresponding SD; , the detection probability of the T-firings; and , the F-firings to T-firings ratio. Following [14], the IDI PDF can be used to obtain a maximum likelihood estimate of these four parameters. In an experiment of limited duration, only a finite number of firings are identified. Hence the sample of IDIs will contain a limited number of observations. Given a set of IDI observations, , we can obtain the joint probability distribution of the set of observations as

(26) For a given set of parameters, , this is the likelihood function of the parameters given the IDI observations, . The MLE of the distribution parameters, , can be obtained by maximizing the natural logarithm of the likelihood function

C. Truncated IDI PDF Automatic EMG decomposition algorithms usually prevent IDI values from falling below a certain value, . This implies that, in practice, the observed IDI distribution is truncated below . In terms of the PDF, this implies a truncation and a subsequent rescaling to maintain the unit area of the PDF. Specifically, the PDF when truncating the observed IDIs below is given by

301

(27) where the final form of tion and truncation, is

, after applying the approxima-

(28) (23) where and

is a unity step function centered at , is the cumulative distribution function (CDF) of . Given that is much smaller than the smallest TT-IDI values, the PDF area lost in truncation comes from the XF-IDIs component. Additionally, is almost constant for low values, i.e., those from 0 to the range of typical values [Fig. 3(c)], hence, it can be approximated by

(24)

Maximization of (27) and (28) can be achieved by means of optimization techniques (Fig. 4). The parameter space must be bounded, given that and must be nonnegative, and and must be between 0 and 1. The use of an optimization technique is the main reason for using the approximation of derived in Section II-B in order to alleviate the computational cost. Note that the information required to solve (27) and (28) comprises the set of observed IDIs, , and which is a predefined parameter determined by the decomposition algorithm (or zero in the case that no truncation below a minimum IDI has been performed to obtain the set of observed IDIs). To select the starting point of the optimization process, following [14], [20], the mode of the IDI histogram is considered to be a good first guess for the mean value, while 0.2 times this

302

IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING, VOL. 23, NO. 2, MARCH 2015

TABLE I RANGES AND VALUES FOR IDI PDF PARAMETERS

Fig. 4. Example of the IDI PDF estimation obtained applying the MLE apms): the histogram of the IDI set (gray bars) proach to a real IDI set ( is superimposed to the estimated IDI PDF (solid lines), and to the contribution and (dashed lines) to the overall PDF. Note that the of leftmost Gaussian curve in the estimated model (dashed line) corresponds to the estimated distribution for the true IDI PDF.

and IDI mean and SD are calculated from the set of remaining IDIs. uses the same approach as MLE but uses the IDI PDF model in [13] instead, not modeling false positives. As MLE and additionally provide , and MLE provides , identification rate (IDR) and inconsistent percent of IDIs (IP) parameters [5] were also calculated as independent and estimators for comparison purposes. Namely, and . For each estimation result, the performance was evaluated by means of the normalized error (29)

mean value is a good first guess for the standard deviation, given that physiological values of the coefficient of variation (CV) usually range from 0.1 to 0.3 [10]. For detection probability and error probability, we used starting values of 0.5 and 0.05, respectively; these values were selected by experimentation. In the current implementation, the Nelder-Mead simplex optimization method is used with a variable transformation to implement bound constraints. Optimization variables which are constrained by both a lower and an upper bound ( and ) use a sine transformation, while those constrained by only a lower bound ( and ) use a quadratic transformation [23]. III. EVALUATION OF THE MODEL A. Estimation Performance With Simulated MU Firing Patterns 1) Estimation Error: Estimation performance was tested with three simulation experiments. In all the experiments, was kept fixed at 100 ms. In each experiment, one of the other three parameters of the model, , was varied to take five different values, while the other two were randomly drawn from within a given range in each realization (Table I). For each combination of the four parameters, 10 000 trials were carried out. For each simulation trial, a synthetic MU firing pattern was generated as a Gaussian renewal process with and and with enough discharges to fill 10 s. Every individual discharge was discarded with probability to model false negatives. Additionally, false positives were modeled as extra firings, drawn at rate , with a uniform distribution within the 10 s temporal span. Estimated and were obtained by applying the MLE, EFE [20], [14] and MnT [4] algorithms. The EFE algorithm uses an iterative approach to determine a lower and a higher threshold value for the IDIs, while the MnT determines these threshold values as a function of the sample SD. All the IDIs outside the range defined by these threshold values are discarded,

is each one of the four parameters under study, . Given the small number of IDIs in 10 s simulations, the actual values of and were recalculated, after counting the false positives and false negatives of the simulated MU firing patterns, as

where

and

(30)

stands for the number of true positives, is the where number of false negatives, and is the number of false positives. Fig. 5 shows the distribution of the normalized errors for each estimated parameter as a function of , , and . Results show that MLE and EFE have good accuracy, always leading to unbiased estimators for all the estimated parameters. Only when is high does MLE present a noticeable bias towards overestimation of [Fig. 5(a)] and [Fig. 5(d)] and underestimation of [Fig. 5(d)], while EFE presents a noticeable bias in high [Fig. 5(a) and (b)] and high [Fig. 5(i) and (j)] conditions. On the other hand, and MnT present noticeable bias in certain conditions. Especially noticeable are the cases of underestimating [Fig. 5(a), (e), and (i)] and overestimating [Fig. 5(b), (f), and (j)]; MnT bias occurs in high [Fig. 5(a)], low [Fig. 5(b)], low [Fig. 5(e) and (f)], and high [Fig. 5(i) and (j)] conditions. In terms of precision, MLE has always the lower error variability. In the case of EFE, precision is comparable to that of MLE in all conditions except for high [Fig. 5(b)], low [Fig. 5(e) and (f)], and high [Fig. 5(i) and (j)] conditions, where EFE shows a larger dispersion than MLE in the estimation error. presents similar precision to that of MLE for the estimation of [Fig. 5(a), (e), and (i)] but larger dispersion for the estimation of [Fig. 5(b), (f), and (e)]. MnT shows low precision in certain extreme conditions, such as, low [Fig. 5(e) and (f)].

NAVALLAS et al.: INTER-DISCHARGE INTERVAL DISTRIBUTION OF MOTOR UNIT FIRING PATTERNS WITH DETECTION ERRORS

303

Fig. 5. Normalized error of the IDI mean, , IDI SD, , detection probability, , and false positive proportion, , for the three simulation experiments: (a)–(d) varying the IDI SD, ; (e)–(h) varying the detection probability, ; and (i)–(l) varying the false positive ratio, . For each of the simulation experiments, the MLE estimation results (square symbols) are compared to that of the EFE algorithm (filled circle symbols) for the IDI mean and SD estimation, and to the and ). In each case the results are summarized with the IDR and IP (filled circle symbols) for the and estimation (noting that median value (solid line) and the 0.15 and 0.85 percentiles (dashed lines) of the 10 000 simulations carried for each parameter combination. The gap of results in . (l) corresponds to the divergence of the normalized error when

The accuracy of MLE relative to that of IDR and IP is good: IDR and IP show high bias and parameter-dependency. Especially noticeable are the IDR error dependency on , , and [Fig. 5(c), (g), and (k)], and IP dependency on [Fig. 5(l)]. In terms of precision, MLE error dispersion is always smaller than that of IDR and IP. Given that MLE and EFE performance is always either equal to or better than that of MnT and , the remaining experiments are restricted to the MLE and EFE algorithms. 2) Estimation Reliability: For the purposes of this study, estimation has been considered to be reliable if the normalized error is below 10% in 90% of the cases for the IDI mean and in 70% of the cases for the IDI SD [14]. In order to map estimation reliability over the parameter space, the space was sampled in a 15 15 grid ( from 0.1 to 0.3; from 0.3 to 1.0) for 500 simulation trials for each combination of the two parameters. This simulation was repeated for three different values, namely: 0, 0.05 and 0.1. Fig. 6 shows the regions where EFE and MLE may be considered reliable within the space, both for the IDI mean and SD estimation. Results show that the MLE reliability region is wider than that of the EFE algorithm for the three test values of . Noteworthy is the reliability of the IDI SD estimation when is high [Fig. 6(f)], where MLE remains reliable over a considerable region of the space but EFE is not reliable at all.

B. Estimation Performance With Real MU Firing Patterns The algorithm was tested with 84 EMG signals of 10 s duration recorded from the human vastus lateralis muscle with concentric needle electrodes during constant low force isometric contractions. The study was approved by the Clinical Investigation Ethics Committee of Navarra. Informed consent was obtained from all subjects. The EMG signals were completely decomposed in the EMGLab environment [25]. Only MUFPs with high-amplitude MUPs and clearly distinguishable shape were accepted for the study. Nonstationary signals and signals without reliable and complete decomposition were discarded from the analysis, reducing the sample to 126 MU firing patterns. Normality of the IDI sample of the MU firing patterns was tested with a Lilliefors test (significance level: ), and only 88 of the 126 MU firing patterns were found to be compatible with a Gaussian IDI distribution. After stationarity and normality tests, a set of 88 complete MU firing patterns extracted from 56 EMG intramuscular signals was accepted for the study. Each of the real MU firing patterns was corrupted by simulating a detection process of the individual firings with probability and adding false firings according to a uniform distribution over the 10 s span in a proportion . In order to compare the performance of the EFE and MLE algorithms for different levels of false positive contamination, the whole corruption process

304

IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING, VOL. 23, NO. 2, MARCH 2015

algorithm does. For the EFE algorithm, the number of rejections noticeably increases with , while for the MLE algorithm the increase is slight. When comparing to previous results [14], it should be noted that in the former study real trains with uncontrolled classification error rates where employed, and so values should lie in the range of 0 to 0.1. For the EFE algorithm, former results lie within the range of rejections obtained for values between 0.05 and 0.1, which is consistent with the methods employed. For the MLE algorithm, it should be noted that the previously published study did not include modeling of classification errors [14]; hence, the fact that current results are always better than the previous ones suggests that inclusion of classification error modeling is advantageous. Fig. 7 shows outcomes of the Kolmogorov-Smirnov tests together with the SD reliability regions in the parameter space for three different values of . As increases, there is an increase in the number of cases where EFE estimation is rejected but MLE estimation is not. Many of these cases correspond to regions where EFE SD estimation is not reliable, although this relationship is not without exceptions. IV. DISCUSSION A. IDI PDF Model Advantages

Fig. 6. Reliability regions of the MLE and EFE algorithms for the IDI mean space for (a), (c), and (e) and IDI SD (b), (d), and (f) estimation in the ; (c), (d) ; (e), (f) . three different values of e: (a), (b) Note that the reliability of the estimation provided by both algorithms is reduced as the number of false positives increases. Also note that the estimation of the IDI mean is much more reliable than that of the IDI SD, where the regions are considerably reduced.

TABLE II KOLMOGOROV-SMIRNOV TEST REJECTIONS

was repeated for three different values of , namely: 0, 0.05 and 0.1, while was randomly selected for each trial from the range 0.3 to 1. For each of the trials, EFE and MLE were applied to estimate and , and the Kolmogorov-Smirnov goodness-of-fit test was applied (significance level: ) to test the null hypothesis: the complete MU firing pattern conforms to a Gaussian distribution with the estimated parameters and . The results from the Kolmogorov-Smirnov test, as summarized in Table II, indicate a rejection of the null hypothesis in 13.3% of the signals (8/60) for the MLE algorithm and 21.7% of the signals (13/60) for the EFE algorithm when . When , rejections are 15.0% (9/60) and 30.0% (18/60), respectively. When , rejections are 18.3% (11/60) and 48.3% (29/60), respectively. In total, the MLE algorithm has almost half the number of null hypothesis rejections that the EFE

The presented model has been derived to accommodate in the observed IDI PDF the possibility that some of the firings detected during MU firing pattern extraction may be false positives. The general model given by (13)–(16) allows us to incorporate false positives independently of the underlying distribution of the real IDI PDF, , as long as the T-firings can be assumed to be the outcome of a renewal process. The set of T-firings, however, may lack some firings if these are not detected in the detection process. It is in the base-model that such false negatives can be incorporated, as done in previous models [13]. In order to reduce computational cost of the IDI PDF calculations, several approximations have been proposed. These specific approximations allow a fast and precise calculation of the observed IDI PDF in normal EMG decomposition conditions ( ), under a Gaussian model for the real IDI PDF. Additionally, the truncation effect of establishing a minimum allowable IDI in the MU decomposition procedure has been incorporated into the model, and a model-based-estimation has been proposed. The MLE approach enables estimation of two statistics of the MU firing pattern distribution (the mean and SD of the IDIs) as well as two parameters of the detection process (the ratio of false negatives, and the ratio of false positives). The mean and SD IDI estimates obtained by means of MLE have lower bias and standard error than the ones obtained with other methods such as EFE [20]. The current MLE algorithm improves on a previous MLE approach [14] through inclusion of false positives in the model. The above advantages are demonstrated in terms of accuracy and precision of the estimator (Fig. 5), reliability (Fig. 6) and goodness-of-fit test results (Fig. 7). Finally, MLE estimates of the false positive and false negative ratios are obtained simultaneously and provide potentially useful information about the MU decomposition process as it occurs.

NAVALLAS et al.: INTER-DISCHARGE INTERVAL DISTRIBUTION OF MOTOR UNIT FIRING PATTERNS WITH DETECTION ERRORS

305

Fig. 7. Outcomes of Kolmogorov-Smirnov tests for MLE and EFE algorithms. Outcomes are depicted within the parameter space. The reliability limits obtained for the IDI SD estimation are the same as in Fig. 6. For each real signal, is calculated from the complete MU firing pattern statistics, and is calculated as the ratio of firings in the incomplete pattern to firings in the complete pattern. Note that the advantage of the MLE over the EFE algorithm ( symbols) is concentrated in a region with moderate detection probability and low .

B. IDI PDF Model Limitations Real data may present a certain amount of skewness. In this case, a Gamma distribution can, in certain conditions, represent more accurately the underlying IDI PDF [12]. The larger the skewness, the larger the estimation error which can be expected with a Gaussian model. However, reports on real data suggest highly symmetrical distributions [10], [26] or slight positive skewness [24], [27]. Highly skewed data reported in other studies may come from the nonstationarity of the experimental setup [9], [12]. All this evidence points towards a smoothed negative impact of the chosen Gaussian model on the estimator performance. All the tests have been made on 10 s recordings. Recordings of shorter lengths will provide less IDIs within the sample, and this will increase the estimation error variability. However, while some other algorithms discard some IDIs by thresholding [4], [20], the MLE approach takes into account all the IDIs within the sample and this is expected to be advantageous when data are scarce [14]. Regarding computation time, previous results [14] demonstrated that the EFE algorithm can be one order of magnitude faster than . Given that in terms of computational complexity, MLE is similar to , we would expect EFE to be about ten times faster than MLE. However, MLE computation time is fast: on our system it takes less than 10 ms to estimate the IDI statistics of a MUFP extracted from a 10 s recording. C. IDI PDF Model and Previous Developments The current model and its parameters are closely related to those of other approaches published in the recent literature. The model is similar to such other approaches in terms of estimation of MU firing pattern statistics, classification of the firings, or assessment of decomposition certainty. Some of the other recent models accommodate false negatives, while others focus on an attempt to accommodate false positives. The model presented here is the first model to include both types of error. In [14], the IDI PDF model including only false negatives derived in [13] is used to obtain estimates of the IDI mean and SD, and the detection probability, . The MLE model-driven estimator [14] outperforms previous algorithmic estimators available in the literature [4], [13], [20] in the absence of false positives. In [7], the same IDI PDF model is used to obtain a firing-

pattern certainty indicator, which is then applied in the process of classifying firings into decomposed MUs. In [4], a term related to the MU firing statistics is included to calculate the a priori validity of the sequence of detected IDIs. A Gaussian IDI model with no false negatives is used, and a specific term, in the form of a uniform distribution over the observed IDI interval, is included to consider false positives. Such a model, although considering false positives, is not taking into account the interference between the real MU firing pattern and the occurrence of false positives that modify the IDI distributions. The importance of the IDI PDF model presented here is that it is the first model in which the treatment of both false positives and negatives is rigorous. In terms of the distribution parameters that can be estimated with the current model, and provide useful information that other approaches try to obtain in different manner. In [6], firingtime information is used to reduce the possibility of mixed clusters, by calculating the ratio of the number of inconsistencies to the number of IDIs. This parameter provides the same information as the estimate . In [5], several statistical parameters are extracted from the IDI histogram to classify firing-patterns as valid, incomplete or contaminated. In essence, two of these parameters, IDR and IP, estimate the ratio of false negatives and false positives, respectively, and are thus estimates of and . Hence, and can be seen to provide an excellent way to classify firing-patterns into valid, incomplete and contaminated [5] and help determine further decomposition actions. In the model presented here, both and estimates are model-driven and estimated simultaneously to and , and this results in more consistent estimation, as demonstrated by the Kolmogorov-Smirnov tests on real data. V. CONCLUSION An IDI PDF model including accommodation for false positives and false negatives has been analytically derived. The model can be applied within the process of MU firing pattern parameter estimation by means of maximum likelihood estimation. Although simplified equations are obtained on the basis of the assumption that MU firing pattern IDIs follow a Gaussian distribution, the model can be extended to other IDI PDFs as long as the MU firing pattern can be safely considered to be a renewal process.

306

IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING, VOL. 23, NO. 2, MARCH 2015

APPENDIX C

APPENDIX A CALCULATION OF

CALCULATION OF

Following (6) and replacing the conditional distribution according to (5) evaluated at

FOR THE

GAUSSIAN MODEL

Given (14), if we assume (2) as the PDF for the TT-IDIs, we obtain

(31) (36) APPENDIX B CALCULATION OF Following (10) and substituting the corresponding distributions as in (9) and (5)

Each term can be further simplified by noting that the exponential term applied to each Gaussian distribution produces two effects: scaling down the function and shifting the mean. Specifically

(37) where and rewritten as

are defined as in (21). Hence (36) can be

(38) (32)

APPENDIX D CALCULATION OF

Note that, although the integral should be carried over all the possible values of , given that the support of the Beta distribution for is limited to , can only take values in the interval. The last form of the equation can be rewritten as (33) The term

Evaluating (15) at

FOR THE

GAUSSIAN MODEL

leads to (39)

When calculating the integral, we note that for each Gaussian term in , , the integral corresponds to the complementary cumulative distribution function, hence we obtain an error function (erfc) term. This leads to

can be expanded as (34)

(40) For reasons explained in the main text, we eliminate the terms in the previous expression, obtaining (41)

Hence (33) can be expanded as

Tweaking the expression in this way, its contribution to the is overall IDI PDF when evaluated at (35) (42) which can be finally rewritten as (15) by defining (16).

given that the summation converges to 1.

NAVALLAS et al.: INTER-DISCHARGE INTERVAL DISTRIBUTION OF MOTOR UNIT FIRING PATTERNS WITH DETECTION ERRORS

REFERENCES [1] D. W. Stashuk, “Decomposition and quantitative analysis of clinical electromyographic signals,” Med. Eng. Phys., vol. 21, no. 6, pp. 389–404, 1999. [2] D. Stashuk and Y. Qu, “Adaptive motor unit action potential clustering using shape and temporal information,” Med. Biol. Eng. Comput., vol. 34, no. 1, pp. 41–49, 1996. [3] H. Parsaei, F. Nezhad, D. Stashuk, and A. Hamilton-Wright, “Validating motor unit firing patterns extracted by EMG signal decomposition,” Med. Biol. Eng. Comput., vol. 49, no. 6, pp. 649–658, 2010. [4] K. C. McGill and H. R. Marated, “Rigorous a posteriori assesment of accuracy in EMG decomposition,” IEEE Trans. Biomed. Eng., vol. 19, no. 1, pp. 54–63, Feb. 2011. [5] H. Parsaei and D. W. Stashuk, “A method for detecting and editing MUPTS contaminated by false classification errors during EMG signal decomposition,” in Proc. IEEE 33th IEEE Eng. Med. Biol. Soc. Conf., 2011, pp. 4394–4397. [6] H. R. Marateb, S. Muceli, K. C. McGill, R. Merletti, and D. Farina, “Robust decomposition of single-channel intramuscular EMG signals at low force levels,” J. Neural. Eng., vol. 8, no. 6, p. 066015, Dec. 2011. [7] H. Parsaei and D. W. Stashuk, “SVM-based validation of motor unit potential trains extracted by EMG signal decomposition,” IEEE Trans. Biomed. Eng., vol. 51, no. 1, pp. 183–191, Jan. 2012. [8] H. Parsaei and D. W. Stashuk, “EMG signal decomposition using motor unit potential train validity,” IEEE Trans. Neur. Sys. Rehab. Eng., vol. 21, no. 2, pp. 265–274, Mar. 2013. [9] C. J. De Luca, “Physiology and mathematics of myoelectric signals,” IEEE Trans. Biomed. Eng., vol. 26, no. 6, pp. 313–325, Jun. 1979. [10] H. P. Clamann, “Statistical analysis of motor unit firing patterns in a human skeletal muscle,” Biophys. J., vol. 9, no. 10, pp. 1233–1251, Oct. 1969. [11] W. S. Masland, D. Sheldon, and C. D. Hershey, “Stochastic properties of individual motor unit interspike intervals,” Am. J. Physiol., vol. 217, no. 5, pp. 1384–1388, Nov. 1969. [12] K. B. Englehart and P. A. Parker, “Single motor unit myoelectric signal analysis with nonstationary data,” IEEE Trans. Biomed. Eng., vol. 41, no. 2, pp. 168–180, Feb. 1994. [13] K. C. McGill, “A method for quantitating the clinical electromyogram,” Ph.D., Stanford Univ., Stanford, CA, USA. [14] J. Navallas, A. Malanda, and J. Rodriguez-Falces, “Maximum likelihood estimation of motor unit firing pattern statistics,” IEEE Trans. Neural. Sys. Rehab. Eng., vol. 33, no. 3, pp. 460–469, May 2014. [15] K. C. McGill, K. L. Cummins, and L. J. Dorfman, “Automatic decomposition of the clinical electromyogram,” IEEE Trans. Biomed. Eng., vol. 32, no. 7, pp. 470–477, Jul. 1985. [16] D. Stashuk, “EMG signal decomposition: How can it be accomplished and used?,” J. Electromyogr. Kinesiol., vol. 11, no. 3, pp. 151–173, Jun. 2001. [17] J. R. Florestal, P. A. Mathieu, and A. Malanda, “Automated decomposition of intramuscular electromyographic signals,” IEEE Trans. Biomed. Eng., vol. 53, no. 5, pp. 832–839, May 2006. [18] J. Florestal, P. Mathieu, and K. Mcgill, “Automatic decomposition of multichannel intramuscular EMG signals,” J. Electromyogr. Kinesiol., vol. 19, no. 1, pp. 1–9, 2009. [19] H. Parsaei, D. W. Stashuk, S. Rasheed, C. Farkas, and A. Hamilton-Wright, “Intramuscular EMG signal decomposition,” Crit. Rev. Biomed. Eng., vol. 38, no. 5, pp. 435–465, 2010. [20] D. W. Stashuk and Y. Qu, “Robust method for estimating motor unit firing-pattern statistics,” Med. Biol. Eng. Comput., vol. 34, no. 1, pp. 50–57, 1996. [21] A. Stuart and J. K. Ord, Kendall's Advanced Theory of Statistics. Volume 1: Distribution Theory, 6th ed. London, U.K.: Hodder Education, 1994.

307

[22] R. Pyke, “Spacings,” J. Roy. Stat. Soc. B, vol. 27, no. 3, pp. 395–449, 1965. [23] J. D'Errico, Fminsearchbnd Documentation 2005 [Online]. Available: http://www.mathworks.com/matlabcentral/fileexchange/8277-fminsearchbnd-fminsearchcon [24] R. S. Person and L. P. Kudina, “Discharge frequency and discharge pattern of human motor units during voluntary contraction of muscle,” Electroencephalogr. Clin. Neurophysiol., vol. 32, no. 5, pp. 471–483, May 1972. [25] K. C. McGill, Z. C. Lateva, and H. R. Marateb, “EMGLAB: An interactive EMG decomposition program,” J. Neurosci. Methods, vol. 149, no. 2, pp. 121–133, 2005. [26] F. Buchthal, P. Pinelli, and P. Rosenfalck, “Action potential parameters in normal human muscle and their physiological determinants,” Acta. Physiol. Scand., vol. 32, pp. 219–229, 1954. [27] C. J. De Luca and W. J. Forrest, “Some properties of motor unit action potential trains recorded during constant force isometric contractions in man,” Kybernetik, vol. 12, no. 3, pp. 160–168, Mar. 1973.

Javier Navallas was born in Pamplona in 1976. He graduated in 2002, and he received the Ph.D. degree in Telecommunication Engineering from the Public University of Navarra, Pamplona, Spain, in 2008. He has also worked as a Software Engineer. He is presently an Associate Professor in the Electrical and Electronics Engineering Department, University of Navarra. His research interests are modeling of the neuromuscular systems and neural information processing.

Javier Rodriguez-Falces received the Ph.D. degree in electromyography from the Public University of Navarra, Navarra, Spain, in 2007. He is a Professor in the Department of Electrical and Electronic Engineering, Public University of Navarra. He teaches courses in the Master of Biomedical Engineering and in the Bachelor Degree of Telecommunication. His research interests include biomedical signal processing, quantitative analysis of electromyographic signals, modeling, and analysis of muscle electrical and mechanical responses and neuromuscular adaptations to exercise.

Armando Malanda was born in Madrid, Spain, in 1967. He graduated in telecommunication engineering from Madrid Polytechnic University, in 1992. He received the Ph.D. degree from Carlos III University, Madrid, in 1999. In 1992 he joined the School of Telecommunication and Industrial Engineering of the Public University of Navarra, Navarra, Spain. In 2003 he became an Associate Professor in the Electrical and Electronics Engineering Department of the same university. During this period he has been teaching several subjects related to digital signal processing, image processing and biomedical engineering. His areas of interest include the analysis, modeling and simulation of bioelectrical signals, particularly EEG and EMG.