EMG Probability Density Function: A New Way to Look at EMG Signal Filling From Single Motor Unit Potential to Full Interference Pattern

An analytical derivation of the EMG signal’s amplitude probability density function (EMG PDF) is presented and used to study how an EMG signal builds-up, or fills, as the degree of muscle contraction increases. The EMG PDF is found to change from a semi-degenerate distribution to a Laplacian-like distribution and finally to a Gaussian-like distribution. We present a measure, the EMG filling factor, to quantify the degree to which an EMG signal has been built-up. This factor is calculated from the ratio of two non-central moments of the rectified EMG signal. The curve of the EMG filling factor as a function of the mean rectified amplitude shows a progressive and mostly linear increase during early recruitment, and saturation is observed when the EMG signal distribution becomes approximately Gaussian. Having presented the analytical tools used to derive the EMG PDF, we demonstrate the usefulness of the EMG filling factor and curve in studies with both simulated signals and real signals obtained from the tibialis anterior muscle of 10 subjects. Both simulated and real EMG filling curves start within the 0.2 to 0.35 range and rapidly rise towards 0.5 (Laplacian) before stabilizing at around 0.637 (Gaussian). Filling curves for the real signals consistently followed this pattern (100% repeatability within trials in 100% of the subjects). The theory of EMG signal filling derived in this work provides (a) an analytically consistent derivation of the EMG PDF as a function of motor unit potentials and motor unit firing patterns; (b) an explanation of the change in the EMG PDF according to degree of muscle contraction; and (c) a way (the EMG filling factor) to quantify the degree to which an EMG signal has been built-up.

no complete analysis, in terms of the EMG signal's amplitude probability density function (EMG PDF), of the way in which the EMG signal progressively builds-up and the baseline is filled up with motor unit potentials (MUPs) as muscle activity increases.
Knowledge of the EMG signal filling process is of great value in prosthesis control, where analysis of the EMG signal is extensively applied [1], [2], [3] in order to determine the intended degree of muscle activation [4], [5]. One of the main limitations for robust control is related to the stochastic behaviour of the signals [6], and reliable modelling of the EMG signal as a random process is found to be be useful in these applications [7], [8], [9], [10], [11]. Another important application of EMG recruitment analysis deals with motor unit (MU) firing pattern extraction [1] and its use in the investigation of neural drive strategies [5].
Analysis of EMG filling as the level of muscle activation increases is also widely employed in clinical practice in interference pattern analysis [12], [13]. However, this is usually attempted by EMG waveform analysis, by looking for qualitative or quantitative descriptors of the degree to which an EMG signal has been filled [14], such as, turnsamplitude analysis [15], or number of short segments and activity analysis [16], [17].
Previous studies on EMG amplitude distribution have demonstrated that EMG PDF consistently shows shape variation according to contraction level [18], [19], [20], [21]. It has been reported that the EMG PDF lies between a Laplacian and a Gaussian distribution at low contraction levels [18], [19], [22] and tends towards a Gaussian distribution at higher activation levels [8], [20], [21]. On the other hand, when the EMG interference pattern is completely formed, it is assumed that the signal is equivalent to a Gaussian process, in accordance with the central limit theorem [23], [24].
There are several analyses of the EMG signal in the time domain that have been used to model the pattern of EMG recruitment. There are detailed models describing the convolutional theory for EMG build-up [25], [26], models describing the MU firing patterns in steady isometric conditions as a renewal point process [27], and models describing the recruitment and firing rate of the MUs as a function of muscle activation [28], [29]. Although in general this type of modelling assumes randomness, the modelled EMG signals are usually described as being quasi-deterministic [24], while more comprehensive This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ descriptions should be grounded on filtered point process theory [30], [31].
As motor unit recruitment increases, individual MUP contributions can no longer be identified, and the so-called EMG interference pattern arises [13]. In this scenario, the EMG signal can be modelled as a Gaussian random process with a band-limited power spectral density [23]. However, there is evidence that the EMG signal is not a Gaussian random process when the level of contraction is low [9], [11], [18], [19], [20], [21] and few MUs have been recruited [24], that is, before the interference pattern is completely formed [15]. To our knowledge, there is no published derivation of the complete EMG PDF from minimum to maximum contraction or description of why, in terms of EMG signal generation, the distribution changes its characteristics as force increases.
In this work, firstly, we provide a complete analytical derivation of the EMG PDF valid for the whole range of muscle activation, based on MUP waveforms and the convolutional model of the EMG. Secondly, we define a new measure to quantify EMG filling (the EMG filling factor) which is a ratio of the first two non-central moments of the rectified EMG signal. Finally, we illustrate, with the help of experiments with simulated and real signals, the usefulness of the proposed factor; explain why, as the activation level of a muscle increases, the EMG distribution changes; and discuss how the EMG filling factor can be used to track level of muscle activation.

A. EMG Convolutional Model
A widely accepted approach to modelling the EMG signal during static isometric contractions is to use a convolutional model [25]. In this derivation, each active MU contributes to the EMG signal its own distinct MUP waveform p i (t), at each of the time instants that the MU is firing, according to its MU firing pattern f i (t). The MU firing pattern can be expressed as a train of impulses [23] where K i is the number of firings of the ith MU within the time interval under analysis. In static isometric contractions, inter discharge intervals τ ik = t i,k+1 − t ik can be modelled by a renewal point process [27], [33] with τ ik ∼ N(µ τ i , σ τ i ).
The mean inter discharge interval µ τ i is the inverse of the MU firing frequency. The MUP train of the ith MU can be expressed as [23] x where * is the convolution operator. Given the random nature of the firing instants, (2) can be interpreted as the response to the impulses of a marked point process [30], and so (2) can itself be regarded as a filtered point process [31].
The EMG signal can be calculated as the summated contribution of the MUP trains from the M active MUs as [25] x While this model provides a simple and accurate description of the EMG signal in the time domain [23], [34], it does not help us to understand the characteristics of EMG PDF. In some approximations [23], the EMG signal has been modelled as a Gaussian noise process subject to the application of the central limit theorem when the number of active MUs is large enough; while this approach may reflect reality for high-contraction recordings, it does not explain the EMG PDF under low and moderate levels of muscle contraction [7], [9], [11], [22].

B. Derivation of the EMG Probability Density Function
We now present the theory behind our derivation of the EMG PDF valid both when contraction levels, and consequently activity in the EMG signal, is low [19], and when contraction levels are high and the EMG interference pattern is formed [23]. The objective is to find a new way to describe, in terms of the signal's amplitude distribution, the EMG signal and the way that that signal builds up or fills out. In this way, the resulting new tool for EMG analysis will encapsulate the relationship between the shape of the EMG PDF and the activation level of muscle.
Given a MUP, p i (t), the sampling process to obtain the MUP can be regarded as an extraction of the amplitude values of the MUP at random instants. Since the sampling process is not synchronized with the MUP's firing instants, the distribution of the sampling times over the course of the support corresponding to the duration of the MUP is uniform. At each particular time instant, the MUP sample amplitude will be given by p i (t). Hence, the sampling process can be regarded as a random variable transformation of the uniform time distribution through the p i (t) function.
As an example, a very simple model for the MUP waveform can be formed by two triangular functions ( Fig. 1(a)) as in where (·) is the triangular function, a i is MUP amplitude and d i is MUP duration. In order to obtain the sampled-MUP PDF, we have to take into account that sampling of the MUP occurs with equal probability for every time instant within the MUP duration. Random variable transformation of a uniform PDF in the temporal domain ( Fig. 1(a)) leads to a uniform PDF model for the sampled-MUP PDF (Appendix A) where (·) is the rectangular function. Another simple model for the MUP waveform would be one period of the sine function ( Fig. 1(b)) as in  and in this case, the random variable transformation ( Fig. 1(b)) leads to an arcsine distribution for the sampled-MUP PDF The sine MUP model provides a more accurate description of the MUP waveform and a better approximation to the corresponding sampled-MUP PDF of a realistic signal ( Fig. 2(a)). However, the exact MUP PDF would always require a precise transformation on the exact counterpart MUP waveform.
In a MUP train segment of duration τ i s between two MU discharges, and neglecting any source of noise, the MUP itself occupies d i seconds while the remaining τ i − d i seconds of the segment are empty, i.e., have zero amplitude. Hence the MUP train will have a mixture PDF formed by the MUP PDF modelling the amplitude distribution within the MUP duration and a Dirac delta distribution modelling the zero signal value during the interval between MUPs (Fig. 2(b)).
Both contributions are mixed in a proportion determined by the ratio η i which is the probability of a sample of pertaining to the MUP, namely where δ(·) is the Dirac delta distribution, φ i (x) is the MUP train PDF of the ith MU, and η i = d i /µ τ i . The latter can be obtained by using the elementary renewal theorem [32] that, in our context, states that the limit number of firings per unit time is 1/µ τ i . Hence, η i can be obtained as the MUP duration d i divided by the expectation of the inter discharge interval τ i of the MU firing pattern driving the MUP train. Note that an increased firing rate in a single MUP train will simultaneously lower the Dirac delta coefficient (1−η i ) and increase the MUP PDF coefficient (η i ), modulating its relative contribution to the MUP train PDF in (8).
When several MUP trains add up to form the EMG, and neglecting any form of MU synchronization [24], the corresponding operation in the PDF domain is a convolution: as we add independent random variables, the resulting PDFs are convolved. Given two MUP trains, the resulting PDF is where the mixing coefficients have a direct interpretation as the probabilities of sampling a superimposition of the two MUPs, η 1 η 2 ; sampling only the first MUP, η 1 (1 − η 2 ); sampling only the second MUP, (1 − η 1 )η 2 ; or sampling an empty segment In general, for a set of M active MUs, the resulting EMG PDF (Fig. 2(c)) will be formed by the convolution of all the PDFs of all the active MUP trains, namely where the K operator is used to indicate the iterated convolution of the φ i (x) distributions. In this expression, each of the contributing MUPs has its own MUP train PDF, φ i (x), obtained for a given η i from its MUP PDF, θ i (x), which, in turn, depends on its particular MUP shape, p i (x).
This expression can be developed to read which can be viewed as a mixture distribution where the terms involving η values are the mixing coefficients of the Dirac delta, MUP PDFs and convolved MUP PDFs distributions. It is important to note that changes in the firing rates are accommodated in the mixing coefficients of the distributions. When summing several MUP trains to form the EMG signal, increased firing rates will turn into increased η values, causing the mixing coefficients of higher order iterated convolutions to increase. This reflects the fact that the amount of MUP superpositions increases as firing rates increase.
Another important observation about the iterated convolution terms in (11) is that, while it is correct to assert that the support of the PDF after the convolution is the summation of the supports of the convolved PDFs, the more terms that are included in the convolution, the more the resulting distribution tends to get smoothed out and tends to concentrate around its central values (the distribution tends to acquire a Gaussianlike shape).This effect is commonly known as EMG amplitude cancellation.
For the sake of simplicity, if we assume an averageη = η i for all the active MUs, and an average MUP PDF,θ (x), the resulting EMG PDF can be expressed as In this simplified expression, the mixing coefficients of the distributions are directly interpretable as the probability of sampling an empty section of the EMG signal Note that, in this simplified scenario, these probabilities follow exactly a binomial distribution B(M,η).

C. EMG Filling Factor
Having derived the EMG PDF, the next objective is to quantify EMG PDF changes corresponding to changes in muscle activation level. We will present the calculation of the filling factor and show the exact derivation for the MUP, MUP train, and EMG PDFs presented in section II-B.
A simple way to quantify the change in EMG PDF shape as an EMG signal is being filled with the contributions from newly recruited MUs is to calculate a ratio of the first two non-central moments of the rectified signal. We will refer to this ratio as the EMG filling factor This parameter is the inverse of the square of the form factor as introduced in [35] to interpret MUNIX results and more recently used in [36] to derive a modification of the MUNIX procedure.
To begin, in order to calculate R for a single MUP train, we need to calculate it for the MUP. The folded distribution, which results from rectifying the MUP waveform, can be obtained by applying where u(·) is the Heaviside step function. Developing on the example, introduced in section II-B, of a simplified triangular MUP waveform ( Fig. 1(a)) that has a uniform distribution (14) gives Hence, the EMG filling factor of a triangular MUP is 3/4 (Table I).
In the case of the sine MUP ( Fig. 1(b)) with an arcsine distribution as in (5), the folded distribution is a half-arcsine distribution (Appendix B) given by and calculation of the non-central moments (Appendix B) gives µ ′ 1 = 2a i /π and µ ′ 2 = a 2 i /2. Hence, the EMG filling factor of a sine MUP is 8/π 2 (Table I).
For the MUP train PDF obtained in (8), the folded distribution corresponding to the rectified signal, φ F (x), can be calculated as because the Dirac delta distribution is concentrated at x = 0, and hence it does not change on folding.
Calculation of the moments of the folded distribution of the MUP train is straightforward as it can be demonstrated that given that the nth central moment of a Dirac delta distribution is equal to 0. In other words, the non-central moments of any MUP PDF are scaled by factor η when considering the counterpart MUP train PDF. Hence, for the triangular MUP, the MUP train filling factor becomes 3η/4. With a value of η around 0.2 (MUP duration of 20 ms and mean inter discharge interval of 100 ms), the MUP train filling factor is calculated to be 0.15 (1st row in Fig. 4(a)). Using the same values for the sine MUP model, the MUP train filling factor of 8η/π 2 is calculated to be slightly over 0. 16. Calculation of the moments of the folded distribution of the EMG stem from (11), after folding each of the linearly combined distributions. To be precise, where the last equation stems from the general form of a convolution of L MUP PDFs with indices given in the set {I 1 , . . . , I L }. It is important to note that the folding operation must be applied after solving the MUP superimposition driven by each iterated convolution operation.
Since each summated term in (11) refers to an independent contribution to the mixture, the overall moments will follow the same linear combination of the iterated convolution terms. Hence, taking into account (18), (19), and (20), the moments of the folded distribution of the EMG are Due to the intrinsically complex nature of the analytical derivation of the moments of the folded distributions obtained after the iterated convolutions, a descriptive approach will be adopted to interpret the evolution of the EMG filling factor as M, the number of active MUs in (11) and (21), increases.
At low levels of muscle contraction, a Laplacian-like distribution of the EMG amplitude arises [18], [19], [22]. This shape reflects the fact that in low force recordings it is typical to find large sections of the EMG signal with little or no contribution from MUPs [13].
If we assume that the EMG PDF follows a zero-mean Laplacian distribution L(0, b), the folded counterpart would be an exponential distribution Exp(1/b). The non-central moments are µ ′ 1 = b and µ ′ 2 = 2b 2 , resulting in a filling factor of 1/2 (Table I), (4th row in Fig. 4(a)).
When, as a consequence of an increasing number of recruited MUs, the EMG signal builds up and constitutes a full interference pattern [13], the form of the EMG PDF tends towards a Gaussian distribution, as theoretically called for by the central limit theorem [23].

D. Simulation Experiments
Simulation experiments were designed with a view to analysing, as degree of muscle contraction is increased from low to maximum, how the shape of the EMG PDF evolves and how the EMG filling factor changes.
A combination of state of the art models was employed in the simulation experiments. Muscle architecture with MU territory sizing and placement was modelled as in [37], ensuring a uniform overlapping of MU territories and MU fiber density over the muscle cross-section. MU fiber number and MU territory area were modelled as an exponential function of the motor unit index [28], [37]. Complete individual muscle fiber innervation was modelled as in [38] and [37], with MU fractions [39] modelled as in [40] with a uniform distribution of the motor end-plate means of each fraction around the muscle mid-length and a uniform distribution around the mean of the motor end-plates within each fraction. Muscle fiber conduction velocities were modelled with a Gaussian distribution for each MU with the mean determined as an exponential function of the MU index [40] and a fixed 0.1 coefficient of variance.
Recruitment and firing characteristics were modelled as in [29], with the average firing rate for each MU provided as a function of the muscle activation, measured as a percentage of the maximum voluntary contraction (MVC) level. MU firing patterns were modelled as renewal point processes with a truncated Gaussian distribution for the inter discharge intervals [33] with the mean obtained as the inverse of the average firing rate and a fixed 0.15 coefficient of variance. In order to simulate a bipolar surface EMG recording of each MUP, recording of the individual muscle fiber's activity was modelled as in [41], and MUPs were obtained by summating the contributions from the fibers belonging to each MU [42]. The final simulated EMG recording was calculated as the summated convolution of the MUP firing patterns and their corresponding MUPs [25].
With this simulation model, with randomized MU territory placement and randomized MU fiber innervation, we obtained 1000 different muscle realizations. For each of the muscle realizations, 121 EMG signals of 1 s duration sampled at 2 kHz at different MVC levels were simulated. In order to obtain the complete recruitment pattern as force increases, the 121 MVC levels were selected following the exponential distribution of the recruitment thresholds of the 120 MUs and adding an extra level at 100% MVC. For each of the 121 simulated signals, µ ′ 1 and µ ′ 2 were estimated as the non-central sample moments m ′ 1 and m ′ 2 respectively. The non-central sample moments were calculated according to where x[n] is the sampled EMG signal and N = 2000 is the number of samples in each recording. The EMG filling factor R was then calculated as The 121 5th, 50th, and 95th percentiles of m ′ 1 , m ′ 2 , and the EMG filling factor R were calculated at each MVC value for the 1000 simulated muscle realizations available.

E. Real Signal Experiments
Ten voluntary subjects (5 female and 5 male) aged between 20 and 28 years (mean ± SD: 23 ± 2 years) participated in this study. Informed consent was obtained from all subjects. None of them reported any current or recent neuromuscular disorder. The experiments were conducted following the guidelines of the Declaration of Helsinki and were approved by the Ethics Committee Board of the Public University of Navarra, Spain (PI-023/19 approved on 11/11/2019).
Experiments involved gradually increasing the isometric dorsiflexion force of the tibialis anterior muscle. Each subject sat on an adjustable chair in a slightly reclined position with the right foot strapped to a footplate of an ankle ergometer. The plate was inclined at an angle of 45 • relative to the floor and the seat was adjusted so that ankle and knee joint angles were at 90 • and 120 • , respectively.
Surface EMG signals were recorded from the tibialis anterior muscle, using self-adhesive circular surface electrodes (Ag/AgCl, 10 mm diameter, Kendall Meditrace 100). In each experiment, two pair of electrodes were placed in differential (bipolar) configuration (see EMG proximal and distal electrodes in Fig. 1 in [43]), positioned lengthwise over the belly of the muscle. The reference electrodes were located over the tendon of the tibialis anterior, and the ground electrode was placed over the tibia. Before electrode placement, the skin was adequately prepared (shaving, light abrasion with sandpaper, and cleaning with rubbing alcohol). Surface EMG signals were amplified (bandwidth from 10 Hz to 1 kHz) and sampled at 5 kHz using an analog-to-digital conversion system (MP150, BIOPAC, Goleta, CA). The differential EMG signal recorded was obtained by subtracting the proximal and distal EMG signals.
In each experimental session, subjects performed three isometric dorsiflexion ramp contractions of 60 s duration, with a resting interval of 10 min between contractions. In order to obtain signals at different activation levels, the electromyographist gave qualitative indications to the subjects to modulate their muscle contraction. Special care was taken to obtain valid recordings at low and very low contraction levels.
The signal acquired from each ramp contraction, was cut into segments of 1 s duration. Non-central sample moments were calculated for each segment as (22) and (23), and the EMG filling factor was calculated as (24). Figure 3(a) is a scatter plot of the EMG filling factor as a function of the mean absolute amplitude value of all the simulated EMG signals, R(m ′ 1 ), where each dot corresponds to a 1 s EMG signal. There is a clear pattern, and this is confirmed by a plot of the 5th and 95th percentile limits of the observations shown in Figure 3(b). Starting from EMG filling values of around 0.1 to 0.2, when the mean rectified amplitude is low, the R values increase towards 0.5, the reference value of a Laplacian distribution, and then begin to saturate with an asymptote at 2/π, the reference value of a Gaussian distribution. After reaching the saturation region, the only observable change is a gradual increase in the rectified mean amplitude of the EMG signals.

III. RESULTS
Hence, the simulation analysis indicates that the EMG filling factor, R, follows a characteristic curve when plotted as a function of m ′ 1 . The overall form of this EMG filling curve is described by the 50th percentile ( Fig. 3(b)). There is a monotonical increase of the EMG filling factor as the contraction level increases. The curve has two clear sections:  Figure 5. In all of the graphs, the dashed horizontal lines indicate the EMG filling factors of a pure Laplacian (1/2) and a pure Gaussian (2/π) distribution (Table I).
it is almost linear for the lowest contraction values, and there is saturation of the filling factor for the upper half of the curve; between these sections is a transition region.
Results with real signals depicted in Figure 3(c) are in good agreement with the simulation results. EMG filling factor values start around 0.3 for the lowest contraction levels and rise rapidly towards 0.5. As contraction levels increase, the curve flattens and stabilizes and the filling factor approaches saturation at around 2/π. This form for the EMG filling curve was seen in all of the three ramp-contraction trials (100% repeatability) with all of the ten subjects. Although a big effort was made to record signals at low contraction levels, most of the low contraction signals lie in the Laplacian region, and no recordings were obtained with an EMG filling factor below 0.29.
In order to analyze how the EMG signal, EMG PDF, and filling factor change with level of muscle contraction, we selected one of the simulation realizations and looked at six activation values corresponding to 1, 7, 15, 30, 60, and 120 active MUs. The signal, PDF, and filling factor are represented in Figure 4, in which each row corresponds to an activation value.
In terms of the EMG signal ( Fig. 4(a)), the effect of increasing the voluntary contraction level is twofold: the signal progressively fills up with MUPs, and the amplitude of the signal progressively increases as more MUs are recruited. At the lowest contraction levels, when few MUs are recruited, the EMG signal has wide empty regions with low or zero signal level. As recruitment advances, the EMG signal quickly becomes filled, forming the so-called interference pattern. At the highest contraction levels, the interference pattern is seen to be completely formed, signal complexity increases (the number of turns increases) and there is an increase in the EMG amplitude. Figure 4(b) illustrates how EMG PDF changes as muscle contraction increases. At the lowest contraction levels, when there are too few MUPs to fill the EMG signal, the EMG PDF is a semi-degenerate distribution (1st and 2nd rows). This part of the PDF is represented primarily by a Dirac delta distribution. The contribution of the Dirac delta distribution to the EMG PDF becomes smaller as the number of MUPs in the signal increases. As recruitment increases, but before the EMG signal is completely filled (while some segments of low amplitude activity are still present in the signal), the EMG PDF becomes more and more Laplacian-like (3rd and 4th rows). When the EMG becomes filled and the interference pattern is fully formed, the EMG PDF approximates more and more to a Gaussian distribution (5th and 6th rows). The convergence to a Gaussian PDF in this high activity context is in accordance with the central limit theorem [23], [35].
With regard to the curve for the EMG filling factor (Fig. 4(c)), as the degree of muscle contraction increases and the PDF evolves, the value of the EMG filling factor increases. EMG filling factor values start at around 0.15 and rise, initially quickly, to values in accordance with a Laplacian distribution (0.5). This rise in the EMG filling factor continues until saturating at values approximating to a maximum of 2/π, corresponding to a Gaussian distribution.
Three real recordings of one-second duration are presented in Figure 5 together with corresponding reference distributions (Laplacian and Gaussian). The three signals were selected to illustrate different degrees of EMG filling. The first case ( Fig. 5(a)) is an EMG signal with low activity; there is apparently just one contributing MU, and the EMG PDF is almost flat on its support but with a sharp peak in its central part. The EMG filling factor is 0.29, which is in accordance with the fact that the distribution is sharper than the reference Laplacian distribution. However, the observed peak is not as narrow as expected for a Dirac delta distribution. The broadness of the peak is due to additive electrical noise in the recording, and such noise may explain why we were unable to record EMG signals with EMG filling factor values much below 0.3. The second case (Fig. 5(b)) concerns a recording at a higher level of muscle contraction. There are few segments of the EMG signal with low activity. By visual inspection, the corresponding EMG PDF is very close to the reference Laplacian distribution; the EMG filling factor is 0.47, which confirms that the EMG signal lies in the Laplacian region of the filling curve. The third case (Fig. 5(c)) shows a signal that looks like a full interference pattern. The EMG PDF closely resembles the reference Gaussian distribution, as indicated by the EMG filling factor of 0.63. Finally, the results for real signals in Figure 5 are in good agreement with the simulation results in Figure 4. showing how the PDF changes from a semi-degenerate PDF in which the Dirac delta contribution is still noticeable to an exponential-like distribution (dashed lines), corresponding to a Laplacian distribution of EMG signals, and finally to a form akin to a half-normal distribution (dotted lines), corresponding to a Gaussian distribution of EMG signals. Note the changes of scale on both axes. (c) Position of the EMG filling factor on the EMG filling curve (circles), illustrating the monotonic increase in EMG filling factor as the level of contraction increases, and showing differentiated sections of the curve: an almost linear section for the low levels of contraction, a transition section crossing the Laplacian filling factor (dashed lines at 1/2) and a saturation section tending towards the Gaussian filling factor (dotted lines at 2/π) at high levels of contraction.

A. EMG PDF
The analytical model of EMG PDF developed in this work embraces the deterministic EMG convolutional model [25] and the description of the EMG as a Gaussian random process [23] and fills the gap between these two models.
An EMG signal is usually described as quasideterministic [24] when the contraction level is low. At the other end of the recruitment range, the assumption that the EMG signal is a band-limited Gaussian-distributed stochastic process with zero mean [23], [24] is only valid when the interference pattern is fully formed: when there are sufficient MUs recruited to justify application of the central limit theorem [23], [35] to the summation of individual MUP trains. In the current study, we provide an analytical description of EMG PDF that covers the complete recruitment range.
Our results are in agreement with previous descriptions of transitions within the EMG PDF in low to mid force conditions [18], [19], [20]. Specifically, our results indicate that the statistical model encompasses the transitions of the distribution from semi-degenerate, when few MUs are recruited, to Laplacian-like, when the level of contraction is still low but more MUs contribute to the EMG [9], [11], [18], [19], [22], and then to Gaussian-like, when the EMG interference pattern has been formed [8], [18], [19], [23].
When the level of muscle contraction is low-to-moderate, that is, too low for the interference pattern to be fully formed, stationarity and ergodicity do not apply. However, this does not imply that EMG PDF models in this activation range of the distribution are useless to research or clinical practice [19], [20]. On the contrary, they provide an alternative way to track force and muscle activation [7], [8], [9], [10], [11], [21]. The reality of such signals is what may validate the use of non-Gaussian signal processing methods for EMG analysis [18].
Our analytical formulation does not currently take MU synchronization into account. To further consider MU synchronization, we refer the reader back to the coefficients in (9). As described in section II-B, when the two MU firing patterns are independent, the probability of sampling a superimposition of the two MUPs is η 1 η 2 . However, if two MUs are synchronized, close firings that cause superimposition of MUPs will be more likely than under complete independence [24], and so the probability of sampling a superimposition will be greater than η 1 η 2 . If the MU firing rates of each MU do not change, the greater probability of superimposition will also affect the EMG signal in terms of both the empty signal contribution (1−η 1 )(1−η 2 ) and the two individual MUP contributions η 1 (1 − η 2 ) and (1−η 1 )η 2 , causing them to decay. To take MU synchronization into account, the generalization of the EMG distribution in (10) and (11) needs to incorporate, for each set of synchronized MUs, the above changes in the coefficients. Although we do not expect accommodation of a moderate degree of synchronization to have a big impact on the final EMG PDF when activation level is high, further research would be needed in order to determine the impact of synchronization on the EMG PDF in low and moderate activation scenarios.

B. EMG Filling
The EMG filling factor and EMG filling curve presented in this work provide a new way to look at EMG signals recorded under low to moderate levels of contraction.
The EMG filling factor can be regarded as the squared inverse of the EMG-waveform's form factor. In [35], it was shown that the form factor of EMG waveforms stabilizes as the EMG PDF approaches a Gaussian distribution. That is, when an EMG signal is completely filled and the interference pattern is fully formed, the form factor ceases to provide further information on EMG recruitment [35]. However, if the EMG is not completely filled, the form factor changes in relation to the degree of filling.
This change in the form factor, although not a subject of investigation in [35], is key to the work we report here. The form factor is stable and uninformative in the saturation zone of the EMG filling curve (Fig. 3), but variable and hence informative in the curve's linear region, which corresponds to EMG signals at low and moderate levels of muscle activation.
In essence, what we report here is the relation between the degree of EMG filling and the EMG PDF and then a way to quantify the degree of EMG filling through the EMG filling factor (Fig. 4).
Researchers have demonstrated that the form factor can be used in the analysis of surface EMG signals from patients with the neurogenic condition amyotrophic lateral sclerosis in order to detect loss of MU activity [36]. We suggest that such a loss of MU activity can be interpreted as an alteration in the EMG filling curve, and that such an alteration could be described, quantified and tracked in terms of a suitable parameterization of the curve.

C. Limitations and Strengths of the Model
While the analytical EMG PDF is exact under the assumptions made in its derivation, the complete calculation for a given set of active MUs would require the exact MUPs and firing rates. Obtaining, processing and interpreting this data would be a complex process, which highlights the need for a simpler way, such as the filling factor, to encapsulate useful physiological information from the EMG PDF.
The EMG filling factor provides a good indicator of recruitment in low force conditions. However, it saturates relatively quickly as force is increased and may prove less effective in tracking EMG changes after the interference pattern is fully formed; this has been stated to be the case for the form factor [35]. Nevertheless, any pathology that prevents development of a full interference pattern in EMG should still affect the EMG filling factor [36].
An important practical aspect to consider is how the filling curve will be affected by electrode set-up (e.g. surface or intramuscular EMG, electrode derivation) and filtering settings. Both electrodes and filtering will only affect the waveform of constituent MUPs by changing the number of phases (positive and negative lobes of the MUP) or modifying their shapes. However, all MUP phases will always have a sine-like waveform, and hence the sampled-MUP PDFs would be expected to have platykurtic distributions and an EMG filling factor that is high (above 0.637). Recall that low EMG filling factor values are achieved by the emptiness of the signal; i.e., it is the η coefficient when forming the MUP train PDF in (8) that lowers the EMG filling factor to the 0.2 to 0.35 range. Additionally, the iterated convolutions in (11) representing MUP superimpositions will always tend towards Gaussian distributions, as dictated by the central limit theorem. All this suggests, then, that different electrode set-ups and filter settings would not affect the validity of the theoretical results, although they might affect the exact shape of the EMG filling curve. Further studies are needed to exactly determine how the EMG filling curve is affected by the recording method (e.g. surface or intramuscular EMG, electrode derivation, and filter settings) and also to validate the tool with different target muscles.

V. CONCLUSION
The statistical model of EMG filling derived in this work (a) provides an analytically consistent derivation of the EMG PDF as a function of MUPs and MU firing patterns; (b) explains the relationship between level of muscle contraction and EMG PDF; and (c) provides a way (the EMG filling factor) to quantify EMG PDF shape in terms of where it lies on the EMG filling curve and hence a way to estimate the degree of muscle activity.
where, in order to interpret the rectangular function, it has to be taken into account that it is valued as 1 in the (−1/2, 1/2) support of its argument; hence, the support as a function of x can be calculated by solving Hence, the sampled-MUP PDF equals The folded arcsine distribution or half-arcsine distribution can then be defined as f |X | (x; a) = HAS(a) = 2 π √ a 2 − x 2 x − a/2 a .
The noncentral moments of the half-arcsine distribution can be calculated as When n = 1 and solving the corresponding integral, we get and for n = 2 we obtain