## Abstract

Over the past four decades, various methods have been implemented to measure synchronization of motor-unit firings. In this work, we provide evidence that prior reports of the existence of universal common inputs to all motoneurons and the presence of long-term synchronization are misleading, because they did not use sufficiently rigorous statistical tests to detect synchronization. We developed a statistically based method (SigMax) for computing synchronization and tested it with data from 17,736 motor-unit pairs containing 1,035,225 firing instances from the first dorsal interosseous and vastus lateralis muscles—a data set one order of magnitude greater than that reported in previous studies. Only firing data, obtained from surface electromyographic signal decomposition with >95% accuracy, were used in the study. The data were not subjectively selected in any manner. Because of the size of our data set and the statistical rigor inherent to SigMax, we have confidence that the synchronization values that we calculated provide an improved estimate of physiologically driven synchronization. Compared with three other commonly used techniques, ours revealed three types of discrepancies that result from failing to use sufficient statistical tests necessary to detect synchronization. *1*) On average, the z-score method falsely detected synchronization at 16 separate latencies in each motor-unit pair. *2*) The cumulative sum method missed one out of every four synchronization identifications found by SigMax. *3*) The common input assumption method identified synchronization from 100% of motor-unit pairs studied. SigMax revealed that only 50% of motor-unit pairs actually manifested synchronization.

- motor-unit firings
- synchronization
- common input
- synchronization methods

over the past four decades, measurements of synchronization of motor-unit firing instances have been used to infer the existence of common presynaptic inputs to motoneurons. Typically, these measurements have been made from cross-correlation histograms, calculated between the firing instances of paired motor-unit action-potential trains (MUAPTs), as described by Perkel et al. (1967b), a work that has provided the foundation of all subsequent synchronization measurement techniques. By detecting peaks in the cross-correlation histogram, McIsaac and Fuglevand (2007), Nordstrom et al. (1992), and Sears and Stagg (1976) identified firing instances, separated by a fixed latency, that occurred more often than would be expected by chance. They reported that these peaks indicated the presence of synchronization that resulted from common presynaptic inputs shared by the motoneurons, but it had been demonstrated previously by Moore et al. (1970) and Perkel et al. (1967a, b) that cross-correlation peaks could also result from moderately nonstationary firing trains, as well as from refractoriness inherent to each neuron. Therefore, the detection of peaks in the cross-correlation histogram alone was known to be insufficient proof that motor-unit firing instances are synchronized occasionally as a result of common inputs. For this reason, preliminary statistical tests are essential to mitigate the influence of moderate nonstationarities and refractoriness before assessing the degree of synchronization between motor units. However, virtually all synchronization-detection methods have not considered these tests and instead, have applied assumptions or approximations that remain to be proven by empirical data. Indeed, Perkel et al. (1967b) cautioned that even the most basic statistical assumptions can result in false conclusions from the cross-correlation data, such as false detections of synchronization.

The cumulative sum (cusum) detection method, or cusum technique, is one approach commonly used to detect synchronization. This method was originally applied to peristimulus time histograms by Ellaway (1978) to detect changes in motor-unit mean firing rates in response to applied stimuli. Later, studies by Adams et al. (1989), Connell et al. (1986), and Keen and Fuglevand (2004) used the cusum method to identify deviations in the mean value of bin amplitudes from the cross-correlation histogram. These works proposed, but did not prove, that a change in the mean value of the histogram, beyond a preset threshold, was indicative of synchronization. To our knowledge, no one has documented a correlation between changes in the mean value of the cross-correlation histogram and the relatively high-density regions in the histogram that could result from synchronized motor-unit firing instances. Hence, the robustness of the cusum detection method against erroneous detections of synchronization is yet to be established.

After identifying the location of a synchronization peak, many studies rely on the z-score synchronization detection method to compute the statistical significance of the synchronization peak (Sears and Stagg 1976; Wiegner and Wierzbicka 1987). Use of the z-score has been justified by the assumption that bin amplitudes in the cross-correlation histogram of firing instances can be approximated by a normal distribution. Sears and Stagg (1976) and Wiegner and Wierzbicka (1987) advanced this notion by proposing that the Poisson statistics of neuron-firing instances, reported by Cox and Smith (1954), could be approximated with normal statistics. Although normal statistics may describe the firing instances of some neurons, the actual statistics that describe motoneuron firing instances remain disputed (Clamann 1969; De Luca and Forrest 1973; Lippold et al. 1960; Person and Kudina 1972). Nonetheless, studies by Keen et al. (2012), McIsaac and Fuglevand (2007), and Nordstrom et al. (1990) have not considered empirical reports of non-normal motoneuron-firing statistics and instead, applied the z-score to detect synchronization.

Perkel et al. (1967b) illustrated the inadequacy of the z-score method, when they noted that the variance of bin amplitudes in the cross-correlation histogram is greater than that predicted by normal statistics. Specifically, refractoriness inherent to the motoneuron-firing process induces dependence between adjacent bins in the histogram. As a result, bins with relatively high numbers of occurrences are more likely to be followed by bins with relatively low numbers of occurrences. They cautioned that failure to account for the subsequently large variance of bin amplitudes “can lead to false attributions of dependence to cells that are, in fact, firing independently.” In spite of this warning, normal approximations of the cross-correlation histogram bin amplitudes have been used frequently to implement the z-score synchronization detection method.

Even when a statistically significant synchronization peak is not detected in the cross-correlation histogram, some studies nevertheless report the degree of synchronization within a fixed, 11-ms time duration, centered at 0 ms latency (Dartnall et al. 2008; Keen and Fuglevand 2004; Semmler and Nordstrom 1995). This method is based on the notion that common inputs produce synchronized firing instances among all pairs of motoneurons. We refer to this practice as the common input assumption synchronization detection method. It raises three concerns.

*1*) It is a subjective approach. The method assumes that synchronization occurs within ±5.5 ms, centered at a 0-ms latency, even though studies by Datta and Stephens (1990), De Luca et al. (1993), Kirkwood et al. (1982), Schmied et al. (1993), and Semmler and Nordstrom (1995) have demonstrated that synchronization peaks in the cross-correlation histogram exist over latencies spanning ±20 ms, with peak widths ranging from 6 to 40 ms. In fact, differences in motoneuron conduction velocities and innervation locations alone could easily produce synchronization latencies as great as 12 ms between motor units from the first dorsal interosseous (FDI) muscle, as shown by Andreassen and Arendt-Nielsen (1987), Dengler et al. (1988), and Saitou et al. (2000).

*2*) The common input assumption synchronization detection method does not apply the mathematically rigorous and empirically tested results of Perkel et al. (1967b), warning that false conjectures of neuron connectivity can arise from computations that lack a sufficient statistical test for dependent firing behavior.

*3*) The approach occasionally produces negative values of synchronization (Nordstrom et al. 1992), the meaning of which is unclear.

To mitigate these shortcomings, we set out to design an improved approach that is not subject to the drawbacks indicated above. Our statistically based synchronization detection method, which we will refer to as the SigMax detection method, is based on a previous synchronization method developed by De Luca et al. (1993).

Importantly, SigMax does not rely on the assumption that synchronization exists amongst all pairs of motoneurons, and it does not depend on the underlying distribution of the firing instances from each motor unit. The SigMax method is comprised of three tests.

*1*) Test for statistically significant nonstationarities. The limitation of the analysis of synchronization to pairs of stationary MUAPTs is necessary to ensure that synchronization measurements are not biased by moderately nonstationary firing behavior.

*2*) Test for statistically significant dependent firing instances. The dependence test is applied to pairs of stationary MUAPTs. The test is robust to false detections of synchronization that could result from motoneuron refractoriness.

*3*) Test for the most statistically significant peak in the cross-correlation data. Only pairs of stationary MUAPTs with dependent firing instances are tested. The detected peak provides the latency, peak width, and magnitude of synchronization between the MUAPTs.

Synchronization was measured using the SigMax detection method and was compared with synchronization calculated by the z-score detection method, the cusum detection method, and the common input assumption detection method. Our analysis demonstrated that these previously used synchronization methods are subject to additional detections, missed observations, and disparate estimation of synchronization between MUAPTs.

## METHODS

### Experimental Design and Protocol

The experimental design and protocol implemented in this study are described in the accompanying report by Kline and De Luca (2014) and will be summarized here. Six healthy subjects, four men and two women, ages ranging from 21 to 23 yr, all with no known history of neuromuscular disorders, volunteered for the study. Before participating, all subjects read, indicated they understood, and signed a consent form, approved by the Institutional Review Board at Boston University. All experiments were performed on the FDI muscle of the hand and the vastus lateralis (VL) muscle of the lower limb. Isometric force was measured during index-finger abduction and leg extension via load cells. Target trajectories and visual feedback of the isometric contraction force were displayed for the subject on a computer monitor.

The surface electromyographic (sEMG) signals were recorded with a five-pin sensor, described previously in De Luca et al. (2006). The surface sensor was placed on the skin over the center of the muscle belly. Signals from the four pairs of electrodes in the sensor were differentially amplified and filtered with a bandwidth of 20–450 Hz. The signals were sampled at 20 kHz and stored in computer memory for offline data analysis. Before recording data, we measured the maximal voluntary contraction (MVC) force by three brief maximal contractions, each with a duration of 3 s, separated by a rest period of 3 min. The MVC of greatest value was chosen to normalize the force level of all following contractions for later comparison across subjects. Subjects proceeded to track a series of target trapezoidal trajectories displayed on the computer screen with the output of the force sensor. For the FDI muscle, trajectories increased at a rate of 10% MVC/s; were sustained at 5, 10, 15, 20, 25, or 30% MVC for 35 s; and were then decreased back to zero at 10% MVC/s. For the VL muscle, trajectories again increased at a rate of 10% MVC/s; were sustained at 20, 25, 30, 35, 40, or 50% MVC for 35 s; and were then decreased back to zero at 10% MVC/s. The recorded force output was high-pass filtered from 0 to 450 Hz, digitized at 20 kHz, and stored in computer memory for offline data analysis. At least 5 min of rest was allotted between contractions.

### EMG Signal Decomposition and Error Reduction

The sEMG signals from four channels of the decomposition EMG (dEMG) sensor were decomposed into their constituent MUAPTs, using the dEMG algorithms described by De Luca et al. (2006), substantially improved in Nawab et al. (2010), and independently verified with three different methods, including direct visual comparison, by Hu et al. (2013a, b, c, 2014). The output of the algorithm provided the firing instances of all MUAPTs obtained from the decomposition. Each firing instance, as measured by the algorithm, was defined by the time of the greatest absolute value of the action potential.

The occasional errors made by our sEMG signal decomposition algorithm were measured using the decompose-synthesize-decompose-compare validation method, described by Nawab et al. (2010). These errors were then mitigated using the error-reduction technique, described in the accompanying report by Kline and De Luca (2014). In brief, we obtained multiple independent decomposition estimates, each from the sEMG signal after adding Gaussian white noise, equal in root mean square to the baseline noise of the sEMG signal. These estimates were then applied to our error-reduction algorithm to produce a new, more probable estimate of the MUAPTs. We implemented the error-reduction procedure using 30 decomposition estimates for each contraction. Only MUAPTs obtained from decomposition with accuracy >95% were retained for further analysis.

### SigMax Statistical Synchronization Computations

All computations were performed on a 25-s epoch of sEMG signal, recorded during a constant force isometric contraction. Our SigMax detection method used three statistical tests to detect and measure synchronization.

#### Test for statistically significant nonstationarities.

We implemented the now widely used Kwiatkowski, Philips, Schmidt, and Shin (KPSS) test to detect statistically significant, nonstationary MUAPTs (Kwiatkowski et al. 1992). Details of the test are provided in *step 1* of appendix. Because moderate nonstationarities could alter the detection and estimation of dependent firing instances from the cross-correlation data (Moore et al. 1970), it was necessary to perform our synchronization calculations only on stationary MUAPTs.

#### Test for statistically significant dependent firing instances.

Dependence was assessed from the cumulative distribution of the cross-correlations, measured using the recurrence time analysis, described previously by Perkel et al. (1967b), and depicted in Fig. 1. For each pair of MUAPTs, a reference and alternate MAUPT were assigned, with the reference MUAPT having fewer firing instances. Forward and backward recurrence times were measured as the latency between each firing instance of the reference MUAPT and forward and backward firing instance of the alternate MUAPT, denoted by *t*_{f} and *t*_{b}, respectively. Different orders of recurrence times were computed, such that the *i*^{th} order recurrence times, *t*_{fi} and *t*_{bi}, were measured between the reference firing instance and the *i*^{th} forward and *i*^{th} backward firing instance of the alternate MUAPT, respectively, as depicted in Fig. 1. The value of *i* ranged from one to five, such that no more than five forward and backward recurrence times were recorded for each reference MUAPT firing instance.

According to basic statistics theory, if two point processes are independent, then they will be uncorrelated. For the specific case of firing instances and recurrence times from pairs of neurons, McFadden (1962) proved that

*if the stationary firing instances from two neurons occur independently, then the recurrence times will be uniformly distributed (uncorrelated)*.

This is a crucial property of recurrence times, fundamental to the study of synchronization. It enables the detection of synchronization that results from dependent firing instances between motor units. With the use of the contrapositive of the proof by McFadden (1962)

*if the recurrence times are not uniformly distributed (correlated), then the stationary firing instances from two neurons do not occur independently*.

We tested pairs of stationary MUAPTs for significantly correlated firing instances, indicative of dependence. Correlation was evaluated by computing the goodness of fit between the empirical cumulative distribution of recurrence times and the predicted uniform cumulative distribution. The specific details of the procedure used for the goodness-of-fit test are provided in *step 2* of appendix and are summarized here. Recurrence times were divided into different orders at intervals of the mean interpulse interval (IPI) of the alternate MUAPT, as described by the equation
(1)

where *i* ranged from −4 to 4. Recurrence times of MUAPTs with independent firing instances were uniformly distributed within each interval of the mean IPI of the alternate MUAPT. Figure 2*A* shows an example of the empirical cumulative distribution of recurrence times superimposed over the predicted uniform cumulative distribution function. We tested the null hypothesis that the empirical recurrence time data were uniformly distributed using the goodness-of-fit method with the Cramér-von Mises test statistic (Cramér 1928; von Mises 1931). Because of its sensitivity to deviations in the mean rather than the variance of the data from uniformity (Stephens 1974), the Cramér-von Mises test was robust to false detections of synchronization that could result from motoneuron refractoriness. The null hypothesis was rejected at the 0.05 significance level, corresponding to a Cramér-von Mises statistic >0.461, as indicated by Stephens (1970). Recurrence time distributions that deviated significantly from uniformity indicated a statistically significant correlation between the MUAPTs. Significantly correlated firing instances detected between stationary MUAPTs indicated the pair of motor units manifested dependent firing instances.

#### Test for the most statistically significant peak in the cross-correlation histogram.

This test was performed to identify and quantify the amount of synchronization between stationary MUAPTs with dependent firing instances. Specific details of this test are provided in *step 3* of appendix and are summarized here. Each of the nine recurrence time intervals from each pair of MUAPTs was tested separately for synchronization by detecting clusters of recurrence times with a density that exceeded what would be expected due to chance. These clusters, or peaks, could occur at different latencies and last for different durations, or peak widths. Our approach iteratively tested all possible latencies and peak widths of the recurrence time data to identify the most statistically significant occurrence of synchronization. Specifically, the peak widths ranged from 1 ms to one-half of the mean IPI of the alternate MUAPT. For each peak width, we detected the latency that produced the greatest number of recurrence-time occurrences (*k*) and computed the statistical significance of the detection. The peak width that produced the number of occurrences with the greatest statistical significance, beyond the 0.05 significance level, indicated a detection of synchronization. An example of synchronization detected from the empirical cumulative distribution of recurrence times is shown in Fig. 2*B*. For each detection, the synchronization peak width (*W*) and latency (*L*) were recorded for further analysis. The magnitude of synchronization was measured using the synchronization index
(2)

where *k*_{max} was the maximum number of recurrence times detected within the peak of width *W*, and *k̄* was the average number of recurrence times expected by chance within the peak of width *W*, computed as
(3)

where *n* is the number of firing instances of the reference MUAPT, and *m* is equal to the mean IPI of the alternate MUAPT divided by the peak width *W* (see *Eq. 7* in appendix). The synchronization index provided the percentage of firing instances between two MUAPTs that occurred in excess of chance (De Luca et al. 1993).

### Synchronization Measured By Other Methods

Our SigMax detection method addresses the required statistical considerations for identifying the latency, width, and magnitude of synchronization with the greatest statistical significance from each pair of stationary MUAPTs with dependent firing instances. We compared synchronization results obtained from SigMax with those obtained by three other methods published previously that do not incorporate statistical tests necessary to measure synchronization adequately. This comparison was performed on the same set of stationary MUAPTs evaluated using the SigMax method. Synchronization was quantified using the synchronization index in *Eq. 2* for all detection methods tested.

#### Z-Score synchronization detection method.

We followed the practice of assuming a normal distribution of bin amplitudes in the cross-correlation histogram to compute the z-score significance threshold. The baseline mean amplitude of the histogram was measured over a ±200-ms region of the histogram, excluding the ±20-ms region centered ∼0 ms latency, as prescribed by previous methodologies used by Keen and Fuglevand (2004); Schmied et al. (1993); and Sears and Stagg (1976). We then identified peaks in the histogram that exceeded the 0.05 normal significance threshold, corresponding to a z-score > 1.96. The number of significant peaks detected and the corresponding synchronization index, peak width, and latency were compared with the synchronization statistics measured for the same pairs of MUAPTs using our SigMax detection method.

We also tested the underlying assumption used to justify the appropriateness of the z-score for detecting synchronization. As stated in introduction, the z-score method is based on the assumption that motor-unit firing instances are distributed normally. Previously used distribution tests relied on χ^{2} or Kolmogorov-Smirnov (KS) test statistics to identify the deviations of empirical data from the expected distribution (Clamann 1969; De Luca and Forrest 1973). Whereas these methods were the most practical for the time of their use, later work by Stephens (1974) showed that the χ^{2} and KS tests have lower statistical power than other goodness-of-fit test statistics. For example, in the specific case of normal hypothesis testing, D'Agostino et al. (1990) showed that these tests are prone to missed detections of deviations from normality. In contrast, the D'Agostino-Pearson omnibus test (D'Agostino and Pearson 1973) overcomes these drawbacks, not only because it is more sensitive to non-normal statistics, but also, because it quantifies the nature of deviation from normality using skewness and kurtosis statistics. Therefore, we implemented the D'Agostino-Pearson omnibus test to determine if the motor-unit IPIs deviated significantly from a normal distribution. The null hypothesis of normality was rejected at the 0.05 significance level.

#### Cusum synchronization detection method.

This method detects a synchronization peak using the cusum of bin amplitudes from the cross-correlation histogram. We applied the cusum detection method by first measuring the baseline mean of the cross-correlation histogram within ±200 ms, excluding the ±20 ms centered ∼0 ms latency. The cusum data were then computed as the running sum of the difference between the baseline mean and the amplitude of each bin in the cross-correlation histogram [for additional details, see Ellaway (1978)]. The cusum data were normalized by the difference between the maximum and minimum values. Locations where the normalized cusum data crossed the 0.1 and 0.9 thresholds indicated the positive and negative boundaries of a synchronization peak, respectively (Keen and Fuglevand 2004; Schmied et al. 1993). The synchronization peak was calculated between the boundaries of the cusum data and was compared with the synchronization peak obtained from our SigMax detection method.

#### Common input assumption synchronization detection method.

This method computes the magnitude of synchronization from a fixed region of the cross-correlation histogram, regardless of the statistical significance computed for the detected amount of synchronization. According to this approach, synchronization was measured from a fixed, 11-ms region of the cross-correlation histogram, centered at 0 ms latency (Keen et al. 2012; McIsaac and Fuglevand 2007; Semmler and Nordstrom 1995). Results from the common input assumption synchronization detection method were compared with synchronization results obtained using our SigMax detection method.

## RESULTS

### SigMax Synchronization Results

All MUAPTs were obtained using 30 iterations of the error-reduction algorithm performed on 144 recorded sEMG signals; a subset of 36 of these signals was used in the accompanying report by Kline and De Luca (2014) to test thoroughly the error-reduction process. In total, 2,287 MUAPTs were obtained with accuracies >95%; 894 were from 72 FDI contractions, ranging from 5% to 30% MVC; and 1,393 were from 72 VL contractions, ranging from 20% to 50% MVC. In total, statistically significant nonstationarities were detected in 100 or 11.2% of MUAPTs from FDI contractions and 187 or 13.4% of MUAPTs from the VL contractions. The remaining stationary data consisted of 794 MUAPTs with 333,633 firing instances from the FDI and 1,206 MUAPTs with 701,592 firing instances from the VL.

Of the stationary data, the IPIs from 98.6% of MUAPTs from the FDI and 99.8% of MUAPTs from the VL deviated significantly from a normal distribution. Examples of non-normal IPI histograms from several MUAPTs are shown in Fig. 3. Normally distributed IPIs would manifest skewness equal to zero and kurtosis equal to three (D'Agostino et al. 1990), but according to the D'Agostino-Pearson omnibus test, the IPIs of all six MUAPTs in Fig. 3 manifested skewness and kurtosis values that deviated significantly from normality (*P* < 0.0001). On average, IPIs from all MUAPTs tested were positively skewed and had relatively greater occurrences in the distribution tails than expected by normal statistics (Table 1).

In total, 6,453 and 11,283 pairs of stationary MUAPTs were tested for synchronization from FDI and VL contractions, respectively. Figure 4 summarizes the results. Our SigMax detection method found statistically significant synchronization in 42.0% of FDI MUAPT pairs and 54.8% of VL MUAPT pairs. The average synchronization index was 19.8 and 16.9, the synchronization peak width averaged 25.8 ms and 18.5 ms, and the synchronization latency averaged −0.1 ms and −0.3 ms in FDI and VL data, respectively.

#### Effects of decomposition errors on calculations of synchronization.

To illustrate the importance of mitigating decomposition errors before measuring synchronization, we compared results from our SigMax method applied to pairs of MUAPTs, obtained with and without the error-reduction algorithm. Synchronization detections were compared between the same pairs of MUAPTs in both muscles. According to the results shown in Fig. 5, unmitigated decomposition errors resulted in 11.2% or 1,268 additional detections of synchronization and 22.1% or 2,498 missed detections of synchronization. For those MUAPTs between which synchronization was correctly identified, on average, decomposition errors produced a synchronization index that differed by −21.5%, a peak width that differed by −25.1%, and a peak latency that differed by more than ±100% (Fig. 5).

#### Comparison of the z-score synchronization detection method with SigMax.

Figure 6 presents an example of z-score synchronization detections from the cross-correlation histogram of recurrence times for one pair of stationary MUAPTs. The average bin amplitude of the histogram and the 1.96 z-score statistical significance threshold are shown. In total, 12 peaks manifested an average value above the normal significance threshold, indicating 12 separate detections of synchronization from the z-score method. For the same data, the SigMax detection method identified only one statistically significant synchronization peak. When applied to all of the stationary MUAPTs, the z-score synchronization detection method identified 104,662 synchronization peaks from 6,453 pairs of FDI MUAPTs and 177,894 synchronization peaks from 11,283 pairs of VL MUAPTs (Fig. 7). On average, the z-score detection method found 16 synchronization peaks from each pair of MUAPTs. This compares with a single synchronization detection found by SigMax in only 50% of the paired MUAPTs studied.

We evaluated the capabilities of the z-score detection method to estimate synchronization in the most centrally located peak in each cross-correlation histogram. The average values and 95% confidence intervals of the synchronization index, peak width, and latency (shown in Fig. 4) were calculated for the same set of stationary MUAPTs, evaluated using the SigMax method. Relative to the results from the SigMax method, on average, the z-score method produced a synchronization index that differed by −71.9%, a synchronization peak width that differed by −85.8%, and a synchronization latency that differed beyond ±100% (shown in Fig. 7).

#### Comparison of the cusum synchronization detection method with SigMax.

Two examples of the cusum synchronization detection method are shown in Fig. 8. Figure 8*A* provides an example of recurrence time data from the first recurrence interval of MUAPT *pair #1*. The corresponding normalized cusum is plotted in Fig. 8*B*. The 0.1 and 0.9 confidence intervals, typically used in the cusum detection method (Keen and Fuglevand 2004; Schmied et al. 1993), are illustrated. The cusum data crossed the 0.9 confidence interval near 19 ms. However, the 0.1 confidence interval was not crossed within the first recurrence interval shown. The same data analyzed using our SigMax detection method revealed a statistically significant peak centered near 10 ms that spanned 14 ms in width (Fig. 8*C*). Because the cusum method was unable to identify a synchronization peak within the first recurrence interval, MUAPT *pair #1* was marked as a missed detection.

The cross-correlation histogram measured from MUAPT *pair #2* is shown in Fig. 8*D*, with the corresponding normalized cusum displayed in Fig. 8*E*. A synchronization peak was detected by the cusum method between −9 and 4 ms. Relative to the 22.7 ms-wide synchronization peak found by our SigMax method (Fig. 8*F*), the cusum detection method produced a synchronization peak width that differed by −9.7 ms or −42.7%.

When we applied the cusum method to the same set of stationary MUAPTs evaluated using the SigMax method, synchronization peaks were identified from 99.8% of the paired MUAPTs tested. The peak widths of the cusum detections spanned ±200 ms of the cross-correlation histogram tested. However, to evaluate better the differences between the cusum detection method and our SigMax method, we limited the analysis of cusum peaks to only those detected within the first interval of recurrence times. For these data, the cusum method detected synchronization in 43.3% of paired MUAPTs from the FDI and 36.9% of paired MUAPTs from the VL (see Fig. 4). Relative to synchronization detected by the SigMax method, the cusum method produced 6.2% or 549 additional detections of synchronization and 27.5% or 2,447 missed detections of synchronization. When comparing the detections made by both methods, the greatest discrepancies in the synchronization statistics were observed for the peak widths, 95% of which differed within the range of −67.6% to 101% (see Fig. 7).

#### Comparison of the common input assumption synchronization detection method with SigMax.

Two sample synchronization peaks are shown in Fig. 9. For one pair of MUAPTs, the common input assumption method assumed that the synchronization peak was located at 0 ms latency with a peak width of 11 ms, eliciting a synchronization index of −8.14 (Fig. 9*A*). For the same pair of MUAPTs, our SigMax detection method indicated that no statistically significant synchronization peak was present.

The cross-correlation histogram measured from the recurrence times of a second pair of MUAPTs is shown in Fig. 9*B*. The common input assumption detection method assumed that the synchronization peak was located at 0 ms latency with a peak width of 11 ms. This assumed detection yielded a synchronization index of 0.536. However, our SigMax detection method found that the most statistically significant peak in the data centered near 10.4 ms latency with a peak width of 6.78 ms and synchronization index of 5.88.

The common input assumption method detected synchronization peaks in 100% of the pairs of stationary MUAPTs, whereas the SigMax method detected synchronization in only 50% of stationary MUAPT pairs. Compared with the SigMax method, on average, the synchronization index differed by −35.0%, the synchronization peak width differed by −38.8%, and the synchronization latency differed by −100% (see Fig. 7). However, most peculiarly, the common input assumption method produced negative values of synchronization in 2,339 or 13.2% of paired MUAPTs (see Fig. 4).

### Long-term Synchronization

The percentage of paired MUAPTs, detected with synchronization by our SigMax method at each recurrence interval, is shown in Fig. 10, *A* and *C*, for FDI and VL data, respectively. The first recurrence interval manifested the greatest percentage of pairs of MUAPTs with synchronization. Relatively higher-order recurrence intervals had progressively fewer synchronized MUAPT pairs.

Perkel et al. (1967b) observed that peaks in relatively lower orders of recurrence times often produced harmonic peaks in relatively higher orders of recurrence times. Their analysis demonstrated that these harmonic peaks are merely statistical artifacts of recurrence time analysis. Such occurrences do not necessarily indicate unique occurrences of synchronization. Figure 10, *A* and *C*, illustrates the percentage of pairs of MUAPTs with synchronization that did not result from harmonics. All recurrence time intervals beyond the first interval manifested <4% of paired MUAPTs with synchronization. Because our SigMax detection method used a minimum significance threshold of 0.05, as many as 5% of MUAPT pairs could manifest significant peaks as a consequence of random variance in the recurrence time data that were not indicative of synchronization. Therefore, all recurrence intervals with <5% of synchronized MUAPT pairs were considered statistically insignificant detections. The remaining statistically significant detections of synchronization occurred within the first recurrence time interval. The final distributions of the latencies of these synchronization detections are shown in Fig. 10, *B* and *D*. Overall, 95% of the synchronization latency data ranged from −7.4 to 7.2 ms in the FDI and from −6.0 to 5.5 ms in the VL.

## DISCUSSION

This study revealed two physiological findings.

*1*) Our analysis did not support the previously reported assumption that common inputs cause synchronization amongst all motoneurons in a pool.

*2*) Long-term synchronization should not be interpreted as representing a physiological event.

### The Common Input Assumption

We found synchronous and dependent firing behavior in only 50% of the motor-unit pairs studied. Our result differs from that predicted by the notion that common inputs cause synchronization amongst all pairs of motoneurons in the pool of a given muscle. Yet, our finding is consistent with the original work of Sears and Stagg (1976), who concluded that the common input notion “does not mean that each contributing motoneuron necessarily has common presynaptic inputs with every other motoneuron of the pool.” We first reported our doubts concerning the common input assumption in De Luca et al. (1993). More recent studies by Dartnall et al. (2008), Hockensmith et al. (2005), and Keen and Fuglevand (2004) have also reported that only a fraction of motor-unit pairs manifests statistically significant synchronization. However, these later studies did not consider the statistical significance of their results and continued to support the notion of common inputs among all pairs of motoneurons.

Consider the underlying requirements to prove the existence of the common input assumption. Dependence can only be proven between stationary MUAPTs using methods capable of detecting significantly correlated firing instances that are robust to false detections from motoneuron refractoriness. According to basic statistical theory

*if the firing instances of two motor units are independent, then the stationary firing instances are uncorrelated*.

Then, according to the contrapositive

*if the stationary firing instances are correlated, then the firing instances of two motor units are not independent*.

Dependence of the firing instances between motor units can only be proven by first establishing that the firing instances manifest a statistically significant correlation. Although this principle is fundamental to measuring synchronization, it has not been taken into account in the majority of previous studies. Keen et al. (2012), McIsaac and Fuglevand (2007), and Semmler and Nordstrom et al. (1995) proposed that the magnitude of synchronization quantified the strength of common inputs received by motoneurons. Their approach assumed dependence of all motoneurons on common inputs without proving that the firing instances between motor units were correlated. Yet, the inverse of the conditional relationship above indicates that

*if the stationary firing instances are uncorrelated, then the firing instances of two motor units may or may not be dependent*.

Therefore, measurements of synchronization from paired motor units with uncorrelated firing instances do not prove that motoneurons receive common inputs.

Even if Kirkwood and Sears (1978), Nordstrom et al. (1992), and Sears and Stagg (1976) observed significantly correlated firing instances among some pairs of motor units, their observations still would not prove that common inputs cause synchronization. Consider the following simple example. An observer records a person walking on the beach every ∼12 h. The observer also records that the tide of the ocean is highest every time the person walks on the beach. After repeatedly making similar observations over several days, the observer notices a correlation between the person walking on the beach and the high tide. Although the correlation in the data indicates some degree of dependence between the person walking and the high tide, it does not prove that the person causes the high tide. Similarly, correlated firing instances indicative of synchronization between stationary MUAPTs do not prove that synchronized firing instances are caused by common inputs to the motoneurons.

### Long-Term Synchronization

Long-term synchronization with latencies beyond ±20 ms has been reported by Datta and Stephens (1990), De Luca et al. (1993), Kirkwood et al. (1982), Schmied et al. (1993), and Semmler and Nordstrom (1995). With the use of the SigMax detection method, we found that these occurrences are not likely to be a physiological event but are an artifact of false detections produced by two factors inherent to previous synchronization detection methods.

The first factor is the use of an insufficient and relatively low significance threshold to detect synchronization peaks. This effect may be seen in Fig. 6 for the z-score method, where an overwhelming number of synchronous peaks were detected at long latencies. Similar detections of long-term synchronization were obtained in 18% of the motor-unit pairs using a binomial significance threshold, as done by De Luca et al. (1993). A comparison of their results with those presented in this study suggests that their binomial significance threshold underestimated the actual statistical significance of synchronization detections.

The second factor is the presence of harmonics in the cross-correlation histogram (Fig. 10). Specifically, when synchronization is observed amongst first-order recurrence times, there is a greater likelihood that synchronization will also be observed at higher orders of recurrence times. The presence of these harmonics has been well documented in the work of Perkel et al. (1967b) but has not been considered by most studies of synchronization [for examples, see Bremner et al. (1991); Kirkwood et al. (1982); Semmler et al. (2002)]. The error can be avoided by calculating synchronization, exclusively from first-order recurrence times, using an objectively derived and statistically reasoned method to detect statistically significant occurrences of synchronization from stationary MUAPTs with dependent firing instances.

### Quality of Data

In addition to using a more robust and statistically more rigorous approach to measuring synchronization, we had a vetted data set representing physiologically meaningful information relevant to manifestations of synchronization during isometric contractions.

*1)* The data had a known, high accuracy level. The decomposition algorithm used in this study has been validated extensively by De Luca et al. (2006), De Luca and Contessa (2012), and Nawab et al. (2010). Independent verification, using three different methods by Hu et al. (2013a, b, c, 2014), has confirmed that our dEMG algorithms can yield firing instances of MUAPTs having an average accuracy of 95%. Furthermore, we applied the error-reduction algorithm described in the accompanying report by Kline and De Luca (2014) to reduce decomposition errors and improve the estimate of the identification and the location of the firing instances. As may be seen in Fig. 5, our analysis showed that unmitigated decomposition errors result in false identifications, missed detections, and incorrect estimation of synchronization. Other studies of synchronization have reported no measure of the decomposition accuracy. Instead, to address potential decomposition errors, firing instances located amongst superpositions of other action potentials were discarded (Nordstrom et al. 1990, 1992; Semmler and Nordstrom 1995). It is in these locations that the action potentials of different motor units occur at or near the same time, indicating that the firing instances of different motor units are most likely to be synchronized. This tailoring practice would likely contribute to the lesser synchronization levels reported previously.

*2)* The data set was assembled objectively. MUAPTs were not subjectively culled or altered in any manner. Instead, we performed our synchronization analysis on the entire set of MUAPTs obtained from sEMG decomposition within the accuracy bounds we report.

*3)* Combined with our technology for recording and decomposing sEMG signals, our experiments were designed to study synchronization in natural voluntary contractions where the subjects were instructed to maintain a constant force. This paradigm is similar to that used by Contessa et al. (2009), De Luca et al. (1993), and Defreitas et al. (2013). Due to technical limitations, other studies were constrained to study MUAPTs during artificial contractions, where the subjects were typically instructed to manipulate their contractions so as to maintain fixed motor-unit firing rates ∼8–10 pulses/s (Dartnall et al. 2011; Datta and Stephens 1990; Keen and Fuglevand 2004; Nordstrom et al. 1992; Semmler et al. 1997). Force levels were rarely reported, but from the description of the tests, they seemed to be <5% MVC. In our experiments, data were obtained from more natural contractions whose force levels ranged from 5 to 50% MVC.

### Evolution of the SigMax Detection Method

Our approach for measuring synchronized firing instances between MUAPTs considered the consequences of incorrect statistical assumptions cautioned by Perkel et al. (1967b). We addressed their warnings by designing a synchronization method that detected the most statistically significant incidence of synchronization between stationary MUAPTs with dependent firing instances. The tests for stationarity and dependence are fundamental precursory measures for calculating the degree of synchronization between motor units. For the past four decades, all of the methods that have been developed to measure synchronization have not applied these considerations. Some proposed that the motor-unit firing instances had a Gaussian distribution but never tested their assumptions, others established values of the latency and peak width of synchronization without testing statistical significance, and almost all methods did not test for statistical dependence between MUAPTs but were used nonetheless to substantiate the notion that motoneurons depend on common inputs.

The SigMax detection method overcomes the shortcomings of other approaches. It is an extension of the method we described in De Luca et al. (1993). Like the earlier report, it uses binomial statistics to detect synchronization. Our earlier work, although more statistically rigorous than other methods, had a drawback. The method underestimated the probability of detecting regions of statistically significant density in the cross-correlation histogram. To overcome the shortcoming of our earlier work, instead of computing the 95% binomial significance threshold—that each bin has a given number of occurrences—we computed the 95% binomial significance threshold—that any bin has a given number of occurrences. Mathematically, this difference is illustrated by comparing the single binomial probability equation provided in De Luca et al. (1993) with the more detailed series of equations required for our approach, as detailed in *step 3* of appendix. Probabilistically, the difference becomes visible when considering a simple experiment of rolling dice. The probability of rolling three *number 3s* from a six-sided die after *n* rolls is lower than the probability of rolling three of any kind from a six-sided die after *n* rolls. As a result, the 95% binomial significance threshold of rolling three *number 3s* is relatively lower than the 95% binomial significance threshold of rolling three of any kind. This accounts for the lower significance threshold used by De Luca et al. (1993). In contrast, our significance threshold is relatively higher and therefore, less susceptible to the false detections. Specifically, the following improvements were made.

*1)* Our SigMax method does not require any a priori assumptions of the underlying motor-unit firing statistics. This contrasts with the z-score synchronization detection method that depends on a normal distribution of bin amplitudes in the cross-correlation histogram—an assumption that results in erroneous detections of synchronization. Instead, SigMax identifies the most statistically significant detection of synchronization, regardless of the statistics that describe motor-unit firing instances.

*2)* We implemented a test to identify MUAPTs that contained statistically significant nonstationarities—a fundamental, preliminary measure that must precede any synchronization analysis. Only stationary MUAPTs were used for the analysis of synchronization. Data from Moore et al. (1970) demonstrated that moderate nonstationarities hinder the detection and estimation of dependent firing behavior between MUAPTs.

*3)* Pairs of stationary MUAPTs were tested directly for significantly dependent firing instances. Our goodness-of-fit test evaluated by the Cramér-von Mises test statistic ensured that all measurements of synchronization made between stationary MUAPTs were proven to have dependent firing instances. All other reported approaches used to measure synchronization did not apply tests for dependent motor-unit firing behavior.

*4)* We measured synchronization from the cumulative distribution function of cross-correlation data, rather than from the cross-correlation histogram of the data. In so doing, we avoided the false detections of synchronization that could occur from relatively high fluctuations of bin amplitudes as a result of the arbitrary nature of histogram binning. The use of the cumulative distribution function also allowed us to avoid false detections of synchronization from relatively high fluctuations in bin amplitudes that typically result from motoneuron refractoriness (Perkel et al. 1967b).

### Consequence of the Z-Score Synchronization Detection Method

The z-score detection method relies on the assumption that the distribution of the IPIs is Gaussian. Contrary to previous reports by Andreassen and Rosenfalck (1980), Buchthal et al. (1954), and Clamann (1969), the normality of the IPI distribution has been questioned before by De Luca and Forrest (1973), Lippold et al. (1960), and Person and Kudina (1972). The overwhelmingly conclusive results from the motor-unit IPI normality test performed in this study demonstrate that IPIs of MUAPTs in EMG signals, obtained from natural (unconstrained) isometric contractions, are not normally distributed. We are confident in this observation because of the following: *1*) the size of the data from individual contractions was considerably greater (more than one order of magnitude) than that used in previous studies; *2*) we used the D'Agostino-Pearson omnibus test for normality, which has been widely accepted as a superior test than those available to the earlier reports; and *3*) the number of observations used for each normality test was tested previously by D'Agostino et al. (1990) and found to be sufficient to reveal statistics underlying motor-unit firing instances. Some previous reports of the z-score method (Keen et al. 2012; Nordstrom et al. 1992; Sears and Stagg 1976) may have used data sets with a greater number of firing instances, but the D'Agostino-Pearson omnibus test (D'Agostino and Pearson 1973) indicates that such large data sets are not required. In fact, the presumed need for large data sets may have cajoled investigators to collect data under the artificial (constrained) contraction paradigm discussed previously.

### Consequences of the Cusum Synchronization Detection Method

Relative to our SigMax method, the cusum method produced two types of discrepancies when estimating the synchronization peak: *1*) approximately one in every four synchronization peak detections is missed (Fig. 7, *A* and *B*); and *2*) correctly identified peaks manifest peak-width discrepancies that range beyond ±50% (Fig. 7, *C–H*). Although synchronization detections from the cusum method are subject to relatively large estimation discrepancies, the average synchronization values are relatively similar to those measured using SigMax. This result indicates that comparisons of average values may not necessarily reveal the discrepancies produced by a given method. Overall, the discrepancies reported in Fig. 7 demonstrate that changes in the cusum from the mean value of the cross-correlation histogram do not guarantee adequate detections of synchronization.

### Consequence of the Common Input Assumption Synchronization Detection Method

The common input assumption synchronization detection method established de facto that the latency of synchronized firing instances occurs at 0 ms, spanning 11 ms in peak width. Our analysis has demonstrated that this assumption is factually incorrect and resulted in the detection of negative values of synchronization from 13.2% of pairs of FDI and VL MUAPTs. Nordstrom et al. (1992) also reported negative synchronization. The physiological interpretation of negative synchronization has never been explained and remains a confounding detail that poses a physiological conundrum.

Overall, our analysis motivates the need for robust statistical methods when measuring synchronization. Failure to apply necessary fundamental tests, such as those for stationarity or dependence, leads to the discrepancies between synchronization that we observed and synchronization values reported in previous studies. It is only from empirical data obtained by factually substantiated statistical measures that physiological mechanisms of motor-unit control can be revealed properly. As such, the validity of the notion that common inputs cause synchronization of motor-unit firing instances awaits experimental verification.

## GRANTS

Support for this work was provided, in part, by the National Institute of Child Health and Human Development (HD05011/HD/NICHD) and National Institute of Neurological Disorders and Stroke (NS077526-01/NS/NINDS) and by funds from Delsys.

## DISCLOSURES

C. J. De Luca is the president and CEO of Delsys, the company that developed the sEMG decomposition technology.

## AUTHOR CONTRIBUTIONS

Author contributions: C.J.D.L. and J.C.K. conception and design of research; J.C.K. performed experiments; J.C.K. analyzed data; C.J.D.L. and J.C.K. interpreted results of experiments; J.C.K. prepared figures; C.J.D.L. and J.C.K. drafted manuscript; C.J.D.L. and J.C.K. edited and revised manuscript; C.J.D.L. and J.C.K. approved final version of manuscript.

## ACKNOWLEDGMENTS

We are grateful to Dr. R. A. D'Agostino for guidance in development of the statistical procedures used for this analysis. We thank the subjects who painstakingly participated in the experiments.

## Appendix: SIGMAX SYNCHRONIZATION DETECTION METHOD

#### Step 1: Test for Statistically Significant Nonstationarities

Methods that were used previously to detect nonstationary MUAPTs have relied on qualitative analyses, such as observing scatter plots of successive IPIs (Clamann 1969; Masland et al. 1969; Person and Kudina 1972). Although innovative at their time of use, these techniques limit the detection of nonstationarities to specific orders of IPIs. To overcome this shortcoming, we implemented the now widely used KPSS test to detect statistically significant, nonstationary MUAPTs (Kwiatkowski et al. 1992). According to DeJong et al. (1992) and Diebold and Rudebusch (1991), the KPSS test has greater statistical power for detecting nonstationary data than the standard unit root tests presented by Dickey and Fuller (1979). For each train of IPIs that we tested, the null hypothesis was that the motor-unit firing instances separated by different lag times were stationary. We implemented the test with T^{1/2} different lag times, where T is the sample size of IPIs for each MUAPT, as recommended by Andrews (1991) and Kwiatkowski et al. (1992). The null hypothesis was rejected at the 0.05 significance level, corresponding to a KPSS test statistic of 0.463. All MUAPTs that produced a KPSS test statistic >0.463 contained statistically significant, nonstationary firing instances. Only stationary MUAPTs were tested for synchronization.

The IPIs from two example MUAPTs are shown in Fig. 11. Specifically, in Fig. 11*A*, not only was no visible trend apparent in the data, but also, the KPSS test statistic was 0.0643, well below the 0.463 detection threshold for nonstationary MUAPTs. In Fig. 11*B*, the IPIs of a moderately nonstationary MUAPT are shown. Although only a slight positive trend was visible in the IPI data as a function of time, the KPSS test statistic was 1.19, indicating that the MUAPT contained statistically significant, nonstationary firing instances.

#### Step 2: Test for Statistically Significant, Dependent Firing Instances

We assessed dependence of the firing instances between stationary MUAPTs by identifying statistically significant correlations in the recurrence time data using a goodness-of-fit test between the empirical cumulative distribution function of recurrence times and the predicted uniform cumulative distribution function. Further details of this and other goodness-of-fit tests are provided in D'Agostino and Stephens (1986). Specifically, we tested the null hypothesis that if the firing instances of two motor units are independent, then their recurrence times will be uniformly distributed (McFadden 1962). Mathematically, this was represented by the hypothesis

*H*_{0}: *t*_{1}, *t*_{2}, …, *t*_{n} comes from *U*(*t*)

where *n* recurrence times, *t*_{i}, come from a uniform distribution. By detecting deviations of the empirical cumulative distribution function of recurrence times, *F*(*t*), from the predicted uniform cumulative distribution function, *U*(*t*), we could test the null hypothesis and identify motor units with dependent firing instances. We evaluated the goodness of fit between the empirical and uniform cumulative distributions using the Cramér-von Mises test statistic, *W*_{n}^{2}, first presented by Cramér (1928) and von Mises (1931) and then later improved upon by Smirnov (1936, 1937), as
(4)

We chose the Cramér-von Mises test statistic, *W*_{n}^{2}, as opposed to the χ^{2} or KS statistics, because the Cramér-von Mises test statistic tends to have greater statistical power (Stephens 1974). To compute the Cramér-von Mises goodness-of-fit test, we arranged the recurrence times, *t*_{i}, in ascending order, such that

*t*_{1} ≤ *t*_{2} ≤ … ≤ *t*_{n}

We then performed a probability integral transformation (Fisher 1930) of the recurrence time data to obtain the sequence of values *u*_{i}, such that

*u*_{i} = *U*(*t*_{i}); for *i* = 1…*n*

These data were used to compute the Cramér-von Mises test statistic specific for the uniform distribution goodness-of-fit test derived by Pearson and Stephens (1962) and Stephens and Magg (1968) as (5)

Following the statistical analysis of Stephens (1970), we adjusted *W*_{n}^{2} to produce a modified Cramér-von Mises test statistic, *W**, as
(6)

With the use of Table 1 in Stephens (1970), the modified Cramér-von Mises test statistic, *W**, was compared with the derived percentage points of *W** to determine the statistical significance of the goodness of fit. We rejected the null hypothesis, *H*_{0}, at the 0.05 significance level, demonstrated by the conditional

Rejection of the null hypothesis, *H*_{0}, established the existence of statistically significant correlations, indicating dependent firing instances between the stationary MUAPTs.

#### Step 3: Test for the Most Statistically Significant Peak in the Cross-Correlation Histogram

We identified the most statistically significant occurrences of synchronization between stationary MUAPTs with dependent firing instances. Figure 12 illustrates a pictorial demonstration of our algorithmic approach. For each empirical cumulative distribution of recurrence times, we first selected a temporal window of width *W*_{j}, such that
(7)

where *m*_{j} ranged from two to the closest integer, greater than or equal to the mean IPI of the alternate MUAPT (). The window width *W*_{j} was chosen, such that the recurrence time data over the given interval could be sectioned entirely into *m*_{j} equal sections of width *W*_{j}. For each window width, the number of *k* recurrence times within the window was counted at different latencies, *L* (Fig. 12). The latency, *L*, was increased incrementally by 0.1 ms. After computing the amplitude *k* within the window of width *W*_{j} at all possible latencies, we identified the latency, *L*, that contained the maximum amplitude of recurrence times, denoted *k*_{max}.

We computed the probability of finding at least *k*_{max} occurrences within the interval of width *W*_{j}. The probability that any of *m*_{j} equally spaced sections, *bin*_{1,2,…m}, of width *W*_{j} contains at least *k*_{max} occurrences is equal to the union of the probability that each section contains at least *k*_{max} occurrences shown by the equation
(8)

With the use of the inclusion-exclusion principle of probability, we can decompose *Eq. 8* as follows
(9)

Because all sections, *bin*_{1,2,…m}, are of equal width *W*_{j}, and the probability of incrementing each section is equal under the null hypothesis (McFadden 1962), we can simplify *Eq. 9* as follows
(10)

which can be reduced further to (11)

where *n* is the number of firing instances from the reference MUAPT or total number of possible occurrences within a section, or *bin*, of width *W*_{j}. The probability that *i* bins have *k* occurrences can be found using Bayes conditional probability rule as follows
(12)

where (13)

and (14)

where the probability *p*_{i} is defined by
(15)

In such a manner, we calculated the probability *P*_{j} of finding at least *k*_{max} occurrences within any window of width *W*_{j}.

We repeated the above procedure using windows with progressively greater widths defined by *Eq. 7*. For each window of width *W*_{j}, we found the latency *L* that provided the maximum number of occurrences *k*_{max}. Subsequently, we computed the probability *P*_{j} of finding *k*_{max} occurrences within each window of width *W*_{j}. After repeating the procedure for all possible window widths, we identified the *j*^{th} window that provided the minimum probability *P*_{min} from all measured probabilities *P*_{j}. If the minimum probability, *P*_{min}, were <0.05 significance threshold, then that pair of MUAPTs was determined to have statistically significant synchronization within the given interval of recurrence times tested.

Notice that *Eqs. 12*–*14* are well suited for detecting a statistically significant number of synchronized occurrences at relatively high significance levels, such as the 0.05 level used in this study. These equations approach the boundary of their efficacy when evaluating the statistical significance for a measured number of recurrence times relatively close to the number of recurrence times expected by chance (see *Eq. 3*). However, under these conditions, the similarity between the number of recurrence times measured and the number of recurrence times expected by chance would indicate that the firing instances from the two motor units do not occur with statistically significant dependence. Therefore, this boundary condition has virtually no effect on the detections of statistically significant synchronization reported in this study.

- Copyright © 2014 the American Physiological Society