Exact inter-discharge interval distribution of motor unit firing patterns with gamma model

Inter-discharge interval distribution modeling of the motor unit firing pattern plays an important role in electromyographic decomposition and the statistical analysis of firing patterns. When modeling firing patterns obtained from automatic procedures, false positives and false negatives can be taken into account to enhance performance in estimating firing pattern statistics. Available models of this type, however, are only approximate and use Gaussian distributions, which are not strictly suitable for modeling renewal point processes. In this paper, the theory of point processes is used to derive an exact solution to the distribution when a gamma distribution is used to model the physiological firing pattern. Besides being exact, the solution provides a way to model the skewness of the inter-discharge distribution, and this may make it possible to obtain a better fit with available experimental data. In order to demonstrate potential applications of the model, we use it to obtain a maximum likelihood estimator of firing pattern statistics. Our tests found this estimator to be reliable over a wide range of firing conditions, whether dealing with real or simulated firing patterns, the proposed solution had better agreement than other models. Graphical Abstract Model of the MU firing pattern generation and detection: fT,1(τ), IDI PDF of the physiological firing pattern; fT(τ), IDI PDF after modeling undetected firings (false negatives); fS(τ), IDI PDF after modeling classification errors (false positives)


Introduction
Analysis of the motor unit (MU) firing pattern provides invaluable information for EMG analysis [1][2][3][4][5] and the automation and evaluation of EMG decomposition [6][7][8][9][10][11]. The firing pattern under certain physiological conditions can be modeled as a renewal point process, and this approach has been demonstrated in physiological studies to be adequate [1][2][3][4]7]. When a firing pattern is modeled as a renewal point process, the inter-discharge intervals (IDIs) are independent and equally distributed following a certain probability density function (PDF).
Two kinds of error may arise when the MU firing pattern is obtained from automatic decomposition methods [12][13][14][15][16]: false negatives, which are firings not detected in the decomposition process; false positives, which are firings not belonging to the MU but included in its firing pattern, i.e., classification errors. In the literature, there are models of IDI PDF that take into account false negatives [17,18] and false positives [19]; all of these models assume a Gaussian model for the physiological IDI PDF.
Instead of a Gaussian distribution, a more suitable distribution for modeling the physiological IDI PDF is the gamma distribution [4]. This is so because the gamma distribution has nonnegative support, whereas the Gaussian distribution extends into negative values. Thus, models based on a Gaussian distribution can generate calculated IDIs that are negative, while measured IDIs can never be negative [20]. Additionally, the gamma distribution allows for some degree of skewness in the IDI PDF, which may enable models based on the gamma distribution to better reflect reality since physiological evidence indicates that the IDI PDF shows low-to-moderate skewness [2,5,21].
More importantly, the only published model that accommodates false positives [19] provides only an approximation of the IDI PDF. In the current paper, we show that, by applying statistical analysis to model the superposition of renewal point processes [22], it is possible to derive an exact solution of the IDI PDF that incorporates not just false positives but also false negatives.
The aim of the present study is to derive an exact IDI PDF model that takes into account false negatives and false positives, and in which a gamma distribution is used to model physiological IDIs. This paper presents the mathematical derivation of the model, and evaluates its usefulness in the context of estimation of MU firing pattern statistics.

Firing pattern and false negatives
The MU firing pattern can be interpreted as a point process whose events take place in instants {t k } [1]. An MU firing pattern of short duration (< 10 s), under constant contraction conditions [3], and for which any possible MU synchronization is neglected, can be modeled as a renewal process [1]. In such a model, the physiological IDIs, calculated as {τ k } = {t k+1 − t k }, are independent and identically distributed following the same PDF, f T ,1 (τ ) (Fig. 1a).
A shifted gamma distribution, sometimes referred to as the Pearson type III distribution, of the IDI PDF is proposed in [4,21] as a model that enables incorporation of a degree of skewness in the IDI PDF (Fig. 2a). The analytical expression is as follows: where α is the location parameter, β is the scale parameter, and ρ is the shape parameter. These three parameters allow independent control of the physiological IDI mean μ, standard deviation σ , and skewness γ . The exact relationships between the parameters and the IDI PDF moments are [4] as follows: To incorporate false negatives, the resulting point process is modeled as a thinned version [20] of the physiological MU firing pattern (Fig. 1b). In this way, individual firings can be regarded as being independently detected with probability p; hence, the probability that a firing go undetected is 1 − p. When a firing is not detected, the result will be a new IDI that is the summation of consecutive IDIs including undetected firings. The probability of not detecting n−1 consecutive firings but detecting the nth equals p (1 − p) n−1 . Hence, the IDI PDF including false negatives can be expressed [4,17,18] as follows: where f T ,n (τ ) is the PDF of the IDIs obtained by summation of n consecutive IDIs. For the gamma model, the summation of n gammadistributed independent variables with equal β also follows a gamma distribution [23] that is calculated as follows: where x i = τ − α i and y = τ − n i=1 α i . Given that all the IDIs are independent and identically distributed Note that lost firings lead to side lobes of the PDF with means in nμ, and increasing standard deviation √ nσ ; (c) IDI PDF of the superposed process f S (τ ), obtained by joining the thinned process f T (τ ) and the Poisson process f E (τ ) modeling false positives according to (1), the summation of n IDIs is distributed as follows: Hence, Eq. 3 can be developed to obtain the IDI PDF of the thinned gamma process (Fig. 2a) as follows: Note that, according to Eq. 2, the mean IDI for each of the f T ,n (τ ) PDFs is as follows: while the mean IDI of the MU firing pattern process after incorporating false negatives is as follows:

Superposition of renewal processes
The possibility of false positives can be incorporated into the model by regarding the occurrence of false positives as a superposed point process [19]. The work of Cox and Smith [22] allows exact calculation of the PDF intervals between successive events of a superposition of renewal processes. Given a renewal point process with intervals between successive events distributed according to f (τ ), it is straightforward that where F (τ ) is the cumulative distribution function and F C (τ ) is the complementary cumulative distribution function or survival function.
The delay function is defined as the PDF of the intervals measured from a fixed time to the immediately preceding event [22]. The equilibrium delay function, g(τ ), is the delay function a long time after the beginning of the process, and can be obtained [22] as follows: where μ is the mean IDI of the process. This equation can be reversed and combined with Eq. 9 to obtain [22] as follows: We can define the complementary cumulative delay, G C (τ ), as follows: From this definition, we can elicit two relationships needed in later calculations as follows: The equilibrium delay functions are necessary because the superposition of N independent renewal point processes leads to a new renewal point process [22] where the complementary cummulative delay is as follows:

Firing pattern and false positives
If we neglect MU synchronization, firings from different MUs are independent. Hence, if a false positive is a random firing from a MU firing pattern other than that under analysis, the false positive error process can be modeled as a Poisson point process [19] with the following: where λ is the intensity of the false positive point process and can be calculated as follows: where μ E is the mean IDI between false positives and e is the ratio of false positive firings to true positive firings. The equality of F C E (τ ) and G C E (τ ) is due to the memorylessness of the Poisson point process [24].
When only two processes defined by PDFs f T (τ ) and f E (τ ) are superposed, reversing Eq. 14 and applying Eq. 15 leads to the IDI PDF of the superposed process (Fig. 1c) as follows: where the mean IDI for the superposed process is as follows: Deriving Eq. 20 gives the following: and substituting Eqs. 13 and 14 into Eq. 22 leads to the following: Incorporating Eqs. 16, 17, and 18 into Eq. 23 enables calculation of the IDI PDF of any MU firing pattern, independently of its model, when superposed to a Poisson process modeling false positives. The resulting equation is as follows: The survival function of the thinned process can be calculated as (Appendix A.1) follows: where (Appendix A.2) The complementary cumulative delay of the thinned process can be calculated as (Appendix A.3) follows: where (Appendix A.4) Finally, Eqs. 25 and 27 make it possible to develop Eq. 24 so that we can obtain the exact PDF of the MU firing pattern with a gamma model taking into account both false negatives and false positives ( Fig. 2c) as follows:

Model-based maximum likelihood estimation
To test the usefulness of the IDI PDF model provided in Eqs. 29 and 30, a maximum likelihood estimator (MLE) of the model parameters was implemented. In essence, the MLE of the model parameters (μ,σ ,γ ,p,ê) can be obtained for a given set of IDIs, {τ k } K k=1 , by applying optimization techniques [18,19] to maximize the following: The only limit that needs to be specified before being able to perform a numerical optimization of Eq. 31 is the truncation of the infinite summation in Eq. 29. Given that the nth term of the summation is weighted by p (1 − p) n−1 , the truncation error is related both to N (the upper limit of the truncated summation) and p. Even with a detection probability as low as p = 0.2, the weight of the 30th component, K 30 (τ ), is just 3 · 10 −4 , which, being 646 times smaller than the weight of the first component, is negligible in practical terms. Thus, by truncating the infinite summation after the 30th component (N = 30), precise numerical calculation of Eq. 29 is possible; smaller values of N may be chosen if it is deemed appropriate to trade precision for calculation speed.
In the current implementation, optimization was performed using Matlab's patternsearch algorithm (MATLAB, The MathWorks, Inc., Natick, MA, USA) running on an Intel i7-2640M 2.8 GHz PC. Computation times for the MLE were always under 100 ms.
Parameter values chosen for the initial point from which optimization proceeds were similar to those used in other, previous estimation approaches (Table 1).μ is initialized as the sample mode of the IDIs [18, 19,25]; the most frequent value is expected to be close to the peak of the distribution, which is expected to be close to the real mean (Fig. 2c).σ is initialized as 0.2 times the sample mode of the IDIs [18, 19,25], which leads to a coefficient of variation in the middle of the physiological range (from 0.1 to 0.3 [1]).γ is initialized to 0.2, in anticipation of a nonsymmetrical PDF with moderate skewness.p is initialized to 0.5 [18], a value in the lower range of the expected detection probability for automatically detected firing patterns [26,27].ê is initialized to 0.05 [18,19], which is a value in the middle of the range of the expected false positive rates for automatically detected firing patterns [26].
While in Eqs. 29 and 30, the solution is given in terms of μ, α, β, ρ, p, and e, a single set of parameters, such as (μ, σ, γ, p, e) or (α, β, ρ, p, e), can be used as the parameter-space for MLE optimization. If optimizing in the (μ, σ, γ, p, e)-space (hereafter, moment-space), the values of the IDI PDF parameters can be obtained by applying Eq. 2. Conversely, if optimizing in the (α, β, ρ, p, e)-space (hereafter, parameter-space), the IDI PDF moments can be obtained by applying as follows:

Log-likelihood curves
To evaluate any differences between the two optimizationspace approaches, we carried a case-study simulation of  [1,2]); b (see [5]) the sensitivity of log-likelihood as a function of the model parameters. In this experiment, a total of 100 MU firing patterns were simulated with μ = 100, σ = 10, γ = 0.5, p = 0.7, and e = 0.05. For each simulation trial, a synthetic MU firing pattern was generated as a gamma renewal process with μ, σ , and γ and with enough discharges to fill 10 s. To model false negatives, each individual discharge had a probability 1 − p of being discarded. Additionally, false positives were modeled as extra firings, drawn at a rate λ = ep/μ, with a uniform distribution over the 10s temporal span. Log-likelihood curves for each trial were calculated for a range of values for each parameter above and below the parameter's real value while keeping the remaining parameters constant at their real value. In this way, the effect of IDI sample variability on log-likelihood maxima for each parameter was evaluated.

Simulated MU firing patterns
Estimation performance was tested with four simulation experiments. In all of these experiments, μ was kept fixed at 100 ms. In each run of each experiment, one of the other four parameters of the model, (σ, γ, p, e) was varied to take one of six different values, while the other three parameters were kept fixed at the simulation point values (Table 1). For each combination of the four parameters, 3 series of 1,000 trials were carried out. Firing patterns with false positives and negatives were simulated by means of the same procedure as that described in the previous section. Estimations of parameters were obtained by independently applying MLE with a Gaussian model [19] and a gamma model. Estimation with the gamma model was performed twice: once in the moment-space and once in the parameter-space. Each model was tested with one of the three independent simulation series. For each estimation result, the normalized error was calculated. Given the small number of IDIs in 10-s simulations, the actual values of p and e (that is, the values based on occurrence in the simulated patterns) were recalculated [19], after counting the false positives and false negatives of the simulated MU firing patterns, as T P /(T P + F N) and F P /T P respectively, where T P stands for the number of true positives, F N is the number of false negatives, and F P is the number of false positives. The sample distributions of the estimates for each parameter combination were tested for bias by means of a one-sample t test (significance level: α = 0.01). The variances of the distributions were also tested with a one-sided F test for equal variance (significance level: α = 0.01), which indicates whether the variance of one distribution is significantly larger than that of another. This latter test was performed for the three MLE methods in pairs and in both senses (A bigger than B and B bigger than A).
The results were then arranged to determine if the variance of a given method is significantly larger than that of just one of the other methods or both of them.

Real MU firing patterns
The algorithm was tested with 84 EMG signals of a 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 Research Ethics Committee of Navarra. Informed consent was obtained from all subjects. The EMG signals were completely decomposed in the EMGLab environment [28]. The signal to interference ratio (SIR) [29] of the decomposed signals was calculated as a measure of the completeness of the decomposition. Only MU firing patterns with high-amplitude MUPs with a clearly distinguishable shape were accepted for the study. The decomposability index (DI) [15] of the decomposed MUPs was calculated as a measure of the MUP distinguishability. Nonstationary signals and signals without reliable and complete decomposition were discarded from the analysis, reducing the sample to 103 MU firing patterns.
Normality of the IDI sample of the MU firing patterns was tested with a Lilliefors test (significance level: α = 0.005); only 72 of the 103 MU firing patterns were found to be compatible with a Gaussian IDI distribution. The mean, standard deviation, and skewness of sample IDI statistics were calculated for each of the MU firing patterns.
Each of the 103 real MU firing patterns was corrupted by simulating a detection process of the individual firings with probability p, and adding false firings according to a uniform distribution over the 10-s span in a proportion e. The corruption process of each firing pattern was repeated independently five times for each combination of ten different values of p between 0.5 and 1, and ten different values of e between 0 and 0.5. In this way, we obtained 51500 corrupted MU firing patterns simulated to cover 100 different (p, e) combinations. Each of these trials was subjected to MLE estimation with a Gaussian model and to MLE estimation with a gamma model in the parameterspace in order to obtain estimates ofμ andσ , and, with the gamma model, additionallyγ . In view of results obtained from the experiments described in the previous section (see Section 3.2), the moment-space version of gamma-based MLE method was not included in the analysis. The resulting parameter estimates were tested, with the Kolmogorov-Smirnov goodness-of-fit test (significance level: α = 0.05), to determine whether they were in agreement with the null hypothesis that the complete MU firing pattern conformed to a Gaussian distribution with the estimated parameterŝ μ andσ , or to gamma distribution with the estimated parametersμ,σ , andγ .
The percentage of Kolmogorov-Smirnov tests in which the Gaussian and gamma model estimates were not rejected were calculated for each of the 100 (p, e) combinations. These percentages were calculated independently within the two MU firing pattern sets defined by the results of the Lilliefors test. To provide a comparison of the gamma and Gaussian models, we calculated the difference in the percentage of Kolmogorov-Smirnov rejections between both models.

Log-likelihood curves
The log-likelihood curves are depicted in Fig. 3. Curves fall into three groups: first, curves with clear maxima within the estimation range, such curves occurred for parameters such as μ, α, β, and ρ; second, curves with shallower variation within the estimation range, these curves were found for parameters such as σ , p, and e. Finally, curves that were almost flat, such as curves for γ . Directly related to shallowness in log-likelihood curves variation is the dispersion of the solutions (points in Fig. 3) around the real value, i.e., the standard deviation of the normalized estimation error (values in the legends in Fig. 3). This is the explanation of the higher spread of the results for some of the parameters: in decreasing order of standard deviation, Although these results relate to estimation in a single point of the optimization-space, they suggest that the MLE optimization procedure performed better in the parameterspace than in the moment-space.

Simulated MU firing patterns
Results of the experiments testing estimation performance are depicted in Fig. 4, with the estimated parameters arranged in rows and the varying parameters arranged in columns.
In terms ofμ (Fig. 4a-d), andσ (Fig. 4e-h), the Gaussian model presents bias in 45/48 cases, and tended to underestimate these parameters more severely with increasing σ , γ , and e, and decreasing p. The gamma model in the moment-space also presented bias, in 32/48 cases, although this bias is not so noticeable in the percentile plots.  Regardingγ ( Fig. 4i-l), note that the Gaussian model does not provide an estimate of this parameter. The gamma model in the moment-space was biased in 24/24 cases, tending to overestimate the parameter. The gamma model in the parameter-space was biased in 19/24 cases, although this bias is not so noticeable in the percentile plots. The moment-space model had larger variance in 18/24 tests, while the parameter-space model had larger variance in only 3/24 tests. The two gamma models provided values ofγ within ± 100% of the real value.
In the case ofp ( Fig. 4m-p), the gamma model in moment-space showed lower bias than the other models (22/24, 6/24, and 16/24 biased cases for the Gaussian, gamma in moment-space, and gamma in parameter-space, respectively). Although the Gaussian model tended to underestimate the parameter, for high values of e, it overestimated. In terms of variance, the Gaussian model had larger variance in 8/48 tests, the gamma model in the moment-space in 13/48 tests, and the gamma model in the parameter-space 14/48 tests. The three algorithms provided values ofp within ± 5% of the real value. Forê (Fig. 4q-t), all the models were essentially biased (21/24, 15/24, and 17/24 biased cases for the Gaussian, gamma in moment-space, and gamma in parameter-space, respectively). In terms of variance, the gamma models were clearly superior to the Gaussian model (39/48, 6/48, and 6/48 cases of larger variance for the Gaussian, gamma in moment-space, and gamma in parameter-space, respectively). The three algorithms provided values ofê within ± 50% of the real value, with the best results being obtained in low σ and high p conditions.

Real MU firing patterns
Histograms of the SIR of the decomposed EMG signals, and DI of the MUs are depicted in Fig. 5. The SIR is always greater than 27 (5th percentile is 38.4 and median value is 72.7) indicating the completeness of the decomposition. The DI is always greater than 5 (5th percentile is 7.5 and median value is 21.3) indicating the distinguishability of the MUPs from which the MU firing patterns are extracted [15].
Histograms of the physiological IDI sample statistics of the 103 real MU firing patterns analyzed are depicted in Fig. 6. The observed values are in accordance with physiological studies reporting IDI means in the (30, 160) ms interval [1,2], IDI standard deviations in the (5, 25) ms interval [1,2], IDI coefficients of variance in the (0.1, 0.33) interval [1], and IDI skewness in the (− 0.5, 2) interval [5]. The Lilliefors test rejected a hypothesis of normality in 44 of the 103 complete MU firing patterns; in the remaining 59 patterns, normality was not rejected.
Results of Kolmogorov-Smirnov tests of goodness-of-fit of the Gaussian and gamma models' estimates are depicted in Fig. 7. The Gaussian model's estimates fell into two Rejection of the hypothesis that there is no bias is marked with a small symbol under each of the distribution's boxplots (circle, triangle, and square for Gaussian, gamma in moment-space, and gamma in parameter-space, respectively). Rejection of the hypothesis that the variance is not greater is marked with a small symbol (same ones as before) over each of the distribution's boxplots, the symbol being doubled if the variance is significantly larger than that of the other two methods. The omission in the first set of results in t is only apparent and corresponds to the divergence of the normalized error when e = 0 distinct regions (Fig. 7a, d): first, when e < 0.3, the lower the e and the higher the p, the better the estimate; second, when e > 0.3, the estimate quality was almost independent of p and degradation of estimation was mainly a result of an increase in e. In the case of the gamma model (Fig. 7b,  e), there was no differentiation of behavior as a function of e, and it was generally the case that the lower the e and the higher the p the better the estimate. These observations are corroborated by the figures for the differences in rejection percentages (Fig. 7c, f), which indicate that the main advantage of the gamma model over the Gaussian model in this respect appeared when e > 0.3 and this advantage increased as e increased. For both models, there was a general loss in estimate quality when dealing with MU firing patterns whose normality had been rejected by the Lilliefors test. Quantitatively, normality-rejected firing patterns resulted in, on average, 15-20% more rejected estimates than did nonnormality-rejected firing patterns. With regard to the differences in estimate rejection percentages between models, the advantage of the gamma model over the Gaussian model was 5-10% greater when dealing with normality-rejected firing patterns than when dealing with nonnormalityrejected ones (Fig. 7c, d).

Mathematical model of the MU firing pattern
The IDI PDF model that we developed in Section 2.1 and that assumes a gamma distribution for physiological IDIs has two important advantages over previous models: firstly, it is based on an exact analytical calculation, and this enables calculation of the PDF with unprecedented precision; secondly, it introduces a skewed distribution to model an IDI distribution. Several important features stem from each of these advantages.
The exactness of the model allows it to be used without restricting the values of its parameters. The validity of the earlier model [19] was restricted to conditions in which e < 0.2, that is, to circumstances in which false positives represented less than 17% of the overall IDIs. This limitation was due to the loss of precision for higher values of e. The existence of this shortcoming was corroborated in the present study; in experiments with both simulated and real signals, estimates with the Gaussian model worsened with increasing e (Fig. 4d, h, p, and t), and the goodness-of-fit plummeted as e increased beyond 0.3 (Fig. 7a, d). While this limitation may not be troublesome with some decomposed MU firing patterns of high quality [26,27], it can be a severe problem with recordings taken under conditions of medium-to-high force muscle contraction, where it is not always possible to completely identify all the discharges due to superposition, and algorithms generally make more classification errors [10,11]. As demonstrated in our experimental results, the current model, being exact, can be employed even in high e environments: Histograms of the physiological IDI sample statistics of the 103 real MU firing patterns that were analyzed: a IDI mean, μ; b IDI standard deviation, σ ; c IDI skewness, γ ; and d IDI coefficient of variance, σ/μ. The observed values are in accordance with physiological studies, which report IDI means in the (30, 160) ms interval [1,2], IDI standard deviations in the (5, 25) ms interval [1,2], IDI coefficients of variance in the (0.1, 0.33) interval [1], and IDI skewness in the (−0.5, 2) interval [5]  estimates showed neither significantly more bias nor greater variance when e increased (Fig. 4d, h, p, and t), and the goodness-of-fit remained dependent on both p and e, even for e > 0.3 (Fig. 7b, e).
The use of a gamma model instead of a Gaussian model for the physiological MU firing pattern has two main benefits. On the one hand, the gamma distribution is strictly nonnegative and hence can be legitimately used as an inter-event distribution of a renewal point process in Eq. 1 [20]. The Gaussian distribution, because it allows interevent values to be negative, can not be used unless suitably modified. Although the Gaussian distribution is commonly used because of the simplicity of the equations [4,[17][18][19], it is never formally correct.
In addition, use of a gamma distribution allows for a certain degree of skewness that may better reflect what is known about real MU firing patterns than a symmetric (zero skewness) distribution [4]. Although experimental evidence is not conclusive, different researchers have reported a small-to-moderate degree of skewness in the IDI distribution [2,5,21]. Whether a product of firing dynamics [5] or of unsteadiness in recording conditions [3,4], variability needs to be accommodated, and accepting the possibility of a degree of skewness in the IDI distribution can help in this respect.

Application to MU firing pattern estimation
In the set of real signals employed in the experiments reported here, a moderate degree of skewness was observed to be present (Fig. 6c) in MU firing patterns obtained in recording conditions that would traditionally be regarded as providing firing pattern stationarity [25]. Furthermore, over 40% of the MU firing patterns (44/103) were not compatible with the normality assumption in the Lilliefors test, and the Kolmogorov-Smirnov goodness-of-fit test consistently indicated that gamma model estimates had better fit than Gaussian model estimates (Fig. 7c, f).
With optimization procedures, special care must be taken to select the most appropriate parameter-space to evaluate (29). In our simulation experiments using the gamma model, results suggest that the parameter-driven MLE performed better than the moment-driven MLE in terms of the bias but worse in terms of the variance (Fig. 4). Summarizing, the gamma model in the moment-space was biased in 77/120 cases while the gamma model in the parameterspace was biased in 67/120 cases; significantly larger variance occurred in 81/216 cases and in 107/216 cases, respectively.
From the results obtained in the simulation experiments, we can ensure that, for the same amount of data available (length of the MU firing pattern), IDI PDF mean and standard deviation estimates obtained from the gamma model are less prone to bias than the estimates obtained from the Gaussian model at the cost of an increased variance. Hence, the relative amount of data to obtain reliable estimates of the gamma model parameters when applied to EMG decomposition should be of no relevance, given that it is has been shown to be superior to the Gaussian model-based estimates in 10 s recordings, which are typical in this kind of experiments. However, the stationarity requirement can be more problematic, given that IDI PDF asymmetry can arise from the nonstationarity of the MU firing pattern [4,30]. Hence, careful assessment of MU firing pattern stationarity is always needed when trying to apply a stationary model to MU firing pattern parameter estimation [4].

Further applicability of the model
Applicability of the IDI PDF model presented here is not limited to MU firing statistics estimation: the model can be useful in many statistically rigorous calculations concerning MU firing patterns, for example, in the development of EMG decomposition algorithms and in evaluation of EMG decomposition. Many automatic algorithms for EMG decomposition incorporate IDI PDF models [11,[13][14][15]28], and IDI PDF statistics have been used to validate the MU firing patterns extracted by automatic EMG decomposition methods [27] and to derive a rigorous a posteriori calculation of EMG decomposition accuracy [7].
Besides being applied into new frameworks, new IDI PDF models can be easily implemented from Eq. 24 to accommodate other distributions for the physiological MU firing pattern, such as the Weibull distribution [30]. The resulting models will then incorporate modeling of both false negatives and false positives. In the same way, if a simpler, truncated-Gaussian model that overcomes the negative support of the distribution is developed, a Gaussian IDI can be implemented, reducing the physiological model parameters to μ and σ .
Although a Poisson distribution is a general and effective approach to model the false positive error process PDF [7,19], it may not fulfill the criteria used during EMG signal decomposition. Hence, the false positive error process PDF could also be refined to fit the particular restrictions of a specific EMG decomposition algorithm, e.g., in [19], a correction to include a minimum allowable IDI value has been incorporated in terms of a truncation of the distribution.

Conclusions
The presented IDI PDF model based on a gamma distribution is exact, is strictly nonnegative, and allows the introduction of a controlled degree of skewness into the physiological IDI distribution. Our test experiments demonstrate the feasibility of deriving an accurate MLE of physiological and detection parameters, and indicate that such a MLE can provide better results than previous models can.