Alterations of Motor Unit Characteristics Associated With Muscle Fatigue

This study aims to characterize motor unit (MU) features associated with muscle fatigue, using high-density surface electromyography (HD-sEMG). The same MUs recruited before/after, and during muscle fatigue were identified for analysis. The surface location of the innervation zones (IZs) of the MUs was identified from the HD-sEMG bipolar motor unit action potential (MUAP) map. The depth of the MU was also identified from the decay pattern of the MUAP along the muscle fiber transverse direction. Both the surface IZ location and the MU depth information were utilized to ensure the same MU was examined during the contraction before/after muscle fatigue. The MUAP similarity, defined as the correlation coefficient between MUAP morphology, was adopted to reveal the alterations in MU characteristics under the condition of fatigue. The biomarkers of the same MUs were compared before/after fatigue (task 1) at 5%, 10%, and 15% maximal voluntary contraction (MVC) and in the process of continuous fatigue (task 2) at 20% MVC. Our results indicate that the MUAP morphology similarity of the same MUs was 0.91 ± 0.06 (task 1) and 0.93 ± 0.04 (task 2). The results showed that MUAP morphology maintained good stability before/after, and during muscle fatigue. The findings of this study may advance our understanding of the mechanism of MU neuromuscular fatigue.


I. INTRODUCTION
M USCLE fatigue is complex due to the various physiological phenomena which contribute to it.Access to physiological data within the muscle or the nervous system could reveal time-dependent changes indicative of a muscle fatigue process [1], [2].Muscle fatigue can be classified into two categories: central fatigue and peripheral fatigue.Central fatigue refers to a decline in nerve drive, which manifests as impaired activation of motor neurons; on the other hand, peripheral fatigue refers to impairment in cross-bridge function and the end of the neuromuscular junction, in the form of a decrease in muscle fiber conduction velocity (MFCV) and a decrease in contractile speeds.Multiple factors may contribute to muscle fatigue, including decreased motor neuron output (If the output of some motor neurons is decreased, the output of other motor neurons may be increased, which will lead to other motor neurons fatigue), increased motor neuron synaptic inhibition, and motor neuron intrinsic fitness; these factors make these neurons less responsive to synaptic excitation during sustained activity [3], [4].
In the past few decades, surface electromyography (sEMG) has been widely used to study the biochemical and physiological alterations of muscle fatigue due to its noninvasive nature and convenience of use.Some researchers have studied the sEMG signal amplitude and frequency alterations during muscle fatigue [5], [6], [7].Some findings [8], [9], [10] suggested that the mean frequency (MF) of sEMG is particularly sensitive to local metabolic changes.As we can see, sEMG is a promising technique to study muscle fatigue.
With the process of muscle fatigue, alterations may occur at the MU level, including MU firing frequency, muscle fiber conduction velocity (MFCV) and MU location [11], [12].By taking advantage of advancements in sEMG decomposition techniques [13], [14], the abundant spatial-temporal information provided by high-density sEMG (HD-sEMG) enables us to study muscle activation down to the MU level.Luca [1] reported that when muscle fatigue occurs, the muscle increases motor unit (MU) recruitment and MU firing frequency to overcome the loss of muscle force, which usually results in an increase in EMG amplitude.Previous findings demonstrated that MU synchronization increased during muscle fatigue [15], [16], and the same outcome was observed for the normalized mutual information [17].
Although HD-sEMG has obvious advantages in the study of muscle fatigue and enables the noninvasive study of muscle fatigue down to the MU level, some questions remain unanswered.For example, how the characteristics of a specific MUAP change during the process of muscle fatigue and whether it changes synchronously with the macro parameters of the HD-sEMG signal.This study aims to answer these questions to better understand the neuromuscular control strategy and the neurophysiological mechanism of muscle fatigue.

A. The Experimental Setup
Five subjects (3 female, 2 male, age 28-36 years, weight 50-65 kg) without neuromuscular disorders were recruited.All subjects gave signed consent approved by the Institutional Review Board of University of Houston, Houston, TX, USA.
Surface EMG signals were recorded from the biceps brachii during isometric voluntary contractions.Force and sEMG signals were recorded simultaneously.The subjects were seated comfortably on a height-adjustable chair and held the handle of the force sensor using their right hands.The elbow was flexed at approximately 135 • and bent naturally and remained in this position throughout the experiment, including rest.The forearm was arranged to lie naturally and was not supinated, or pronated.The forearm was almost parallel to the ground.The experiment included four stages in order: (1) subjects were asked to hold the handle using their right hand for three maximum contractions separated by 2 min rest, and the maximum of these three contractions was considered the maximum voluntary contraction (MVC); (2) subjects were instructed to hold the handle at three levels of 5%, 10% and 15% MVCs.The duration of each contraction was 15 seconds and between them, a full recovery time of 2 min was given; (3) subjects continued to keep isometric contractions at 20% MVC until the force dropped by 50%, even with verbal encouragement; (4) The previous experimental tasks were followed by 5%, 10% and 15% MVC isometric contractions after muscle fatigue, and no break was given between each contraction to record the sEMG signal.During each contraction, the target contraction level, i.e., the force target was displayed on a computer monitor as a red line, and the subjects' force trace was shown on the monitor as a white line.Subjects were asked to match their force trace (white line) to the target red line as precisely as they could.Force signals were recorded at a sampling frequency of 2000 Hz.

B. HD-sEMG Recordings
After skin preparation, a high-density surface EMG electrode array (8 by 8 sensors) with an electrode diameter of 4.5 mm and an interelectrode distance of 8.5 mm in two directions was placed on the muscle belly of the biceps brachii.Surface EMG was recorded using a REFA system at a sampling frequency of 2048 Hz (Enschede, TMSI, Netherlands).The reference electrode was placed on the lateral epicondyle of the humerus, and the ground electrode was attached to the wrist of the contralateral arm with a fully soaked Velcro wrist band (TMSi, Enschede, the Netherlands).

C. HD-sEMG Signal Processing
The stable sEMG signals of 10 seconds are first filtered by a bandpass filter (5-500 Hz) and a notched filter to eliminate the 60 Hz frequency.Then, the K-means convolution kernel compensation (KmCKC) algorithm was employed to decompose the composite sEMG signal into constituent pulse trains [13], [18].In brief, the K-means clustering method is used to classify the peaks detected from the global activity index into three groups.The group with the largest number of firing instants was used to construct the initial innervation pulse train (IPT) by using a linear minimum mean square error method (LMMSE) [19].Then, a modified multistep iterative convolution kernel compensation (CKC) method was adopted to update the estimated IPT.The number of iterations is usually set to 30-40.Finally, spike-triggered averaging (STA) algorithms were adopted to calculate the motor unit action potential (MUAP) map.That is, the MUAP is obtained by averaging the waveforms of a certain length of time before and after each spike of IPT.The bipolar MUAP mapping was obtained by subtracting the consecutive monopolar MUAP from its proximal neighbor in the muscle fiber direction.
According to the expression of LMMSE in literature [19], a statistically linear estimate with HD-sEMG signals can be written: where ŝ j is an estimation of IPT.
] is an expectation.Usually, it is necessary to extend the initial sEMG signals X with L − 1 delayed repetitions to improve the knows-to-unknowns ratio [19].The expression X for can be written as: If all the firing times of IPT are known in advance, the estimated equation for C X θ can be written as [18]: where φ j represents a set of firing times for j − th MU, and car d(φ n ) represents the cardinality of φ j .For the j − th MU, K-means clustering algorithm is employed to classify the firing times, and the one with the most firing times is reserved.All the firing times in this class are calculated by Equation 3, resulting in C X θ .
Examples of a one-channel sEMG signal and an IPT are shown in Figures 1(B) and 1(C), respectively.

D. Identical MU Identification
MUs with the same three-dimensional spatial location were considered the same.The three-dimensional spatial location of the MU can be determined in two steps: (1) characterization of the surface MU innervation zone location from the bipolar MUAP map and (2) determination of the MU depth from the monopolar MUAP map.Since the MUAP initiates from the neuromuscular junction (indicated as IZ) and propagates in two opposite directions toward the fiber endings, the surface location of MU IZ can be visually identified from the bipolar MUAP mapping by observing either the reversal of signal polarity or minimum amplitude [20], [21].In addition, according to previous studies [22], [23], the depth of MU recruited at a low contraction level can be estimated according to the full width half maximum (FWHM), defined as the muscle-fiber-transverse-direction distance over the skin surface where the absolute value of the MUAP negative peak is higher than 50% of the maximal absolute value.Therefore, by combining the surface MU IZ location with the identified depth information, the same MU recruited at different contraction levels can be identified.Figure 1 (A) shows the flow chart for the identification of the same MU.

E. Data Analysis
The biomarkers of the same MU were compared before / after fatigue (task 1) at 5%, 10%, and 15% maximal voluntary contraction (MVC) and in the process of continuous fatigue (task 2) at 20% MVC.
After the same MU before/after fatigue is determined, the similarity of MUAP between the two MUs is studied.The MUAP similarity calculation steps are: (1) The STA algorithm is used to calculate the average waveform of a certain length of time before and after each spike of IPT as a MUAP waveform; (2) Determine the surface MU IZ location (shown in Figure 1D); (3) MU depth was estimated by FWHM [22], [23].As shown in Figure 1E, take the waveform on the electrodes in the column (perpendicular to the direction of the muscle fibers) where the position of the MU IZ is located, and estimate the MU depth; (4) Confirm the same MUs extracted before/after fatigue at the same IZ location and depth; (5) The MUAP similarity is calculated by matching the maximum of two MUAP waveforms.That is, first, the waveforms are aligned based on the maximum, and then the MUAP similarity is calculated.
Note that the " * " shown in Figure 1D is the location of the IZ, and select the unipolar waveform on the channel closest to the IZ as the MUAP to calculate the similarity.In Figure 1D, " * "is in the middle of two channels, either of which can be a MUAP waveform (in the red box in Figure 1E).
The similarity of two MUAPs was evaluated by the correlation coefficient (CC) between the IZ channel (closest to the identified surface IZ location) identified from monopolar MUAP maps of two MUs, defined as follows: where C O V represents the covariance matrix.MUAP1 and MUAP2 represent two MUAPs with the same length of 20 ms.
If MUAP1 and MUAP2 are the same, the CC would achieve its maximum value of 1; if MUAP1 and MUAP2 are opposite, the CC would achieve its minimum value of −1.

A. MU Characteristics Before / After Muscle Fatigue
Figure 2 shows the amplitude and frequency of the 19 MUs from 5 subjects before / after fatigue.Figure 3 shows examples of the same action potential across before/after fatigue.Figure 4 shows the CCs of 19 MUs before / after fatigue.As we can see from this figure, the MUAP similarity is high, even between those before / after fatigue.

B. MU Characteristics During Continuous Muscle Fatigue
Figure 5 shows the time-dependent curves of the mean frequency (MF) and mean power frequency (MPF) in the 20% MVC fatigue test.The slopes of MF and MPF were negative, and MF and MPF decreased with the extension of time.
Figure 6 lists the amplitude and frequency of 21 MUs from 5 subjects at 10-sec windows during continuous fatigue.Even under the same force, the amplitude and frequency of the same MU changed irregularly during continuous fatigue.In other words, the frequency and amplitude of the same MU can increase or decrease with the development of fatigue.This result demonstrates the same pattern as Figure 2.
Figure 7 presents the CC of 21 MUs during continuous fatigue.As we can see from this figure, the MUAP during continuous fatigue remains almost constant, presenting the same pattern as that shown in Figure 4.
Figure 9 shows part of the MU firing time.Figure 10 is an example of a change in the recruiting region.The MUAP amplitudes of each MU on the electrode are extracted, and then all the MUAP amplitudes on the same electrode are superimposed.Finally, convert to a range of 0-1, where 0 corresponds to the minimum and 1 to the maximum.

IV. DISCUSSION
In this paper, by using the KmCKC algorithm, the sEMG signal was first decomposed into constituent pulse trains from which monopolar and bipolar MUAPs can be constructed.The same MU was recruited at different contraction levels before / after fatigue and at different time windows during continuous fatigue.This MU was extracted by jointly considering the surface MU IZ location identified from the bipolar MUAP map and its depth information identified from the monopolar MUAP map.The MUAP similarity was studied and compared to reveal the characteristics of the same MU in the state of muscle fatigue.
An important contribution of this paper is to verify that the morphology of the MU maintains good stability before, during and after fatigue.This phenomenon clearly shows that the same MU under different levels of muscle contractions presents high similarity.This conclusion is consistent with previous findings that a MUAP can be used to track the same MU within and across test sessions [24], [25], [26], [27].A possible explanation for this phenomenon is that MUAP is mainly related to individual factors such as muscle structure, Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.length and shape, and these factors remain stable in isometric contractions, even in fatigue.Therefore, MUAP may be a reliable diagnostic tool to use to assess the status of the neuromuscular system [28].
In Figure 2, there were three types of changes in the firing frequency of the same MU before / after muscle fatigue: (1) Most MUs decrease in firing frequency and increase in amplitude.The decrease in firing frequency is due to fatigue.The increase in amplitude is probably due to the post-fatigue recovery period of MU.
(2) The firing frequency and amplitude of a few MUs decrease.The decrease in firing frequency is due to fatigue.The decrease in amplitude is probably due to the fact that MU has not yet entered the recovery period after fatigue.
(3) The firing frequency and amplitude of a few MUs increase.This reason is probably due to MU non-fatigue.That is, muscle fatigue and MU fatigue process are not consistent.
In Figure 5, MF and MPF decreased during muscle fatigue.At the same time, it can be seen in Figure 6 that the frequency changes of the same MU (21 MUs in total) can be divided into two categories as follows: (1) The firing frequencies of MU (15 MUs in total) decrease with muscle fatigue progression, which is divided into three cases here: ① The firing frequency decreases continuously, such as MU4, MU11, MU12, MU15, MU20, and MU21.It shows that these MUs are synchronous with the process of muscle fatigue; ② The firing frequency increases, then decreases, such as MU1, MU2, MU3, MU6, MU7, and MU13.The recruitments of these MUs are shown in Figure 9, which shows that they are recruited later.This was the reason for temporally selective MU recruitments, that is, some MUs are recruited early in the contractions, and some MUs are recruited later; ③ the firing frequency is stable at first, then increases, and finally decreases, such as MU5, MU16, and MU18.This indicates that these MUs are not fatigued in the early stage but begin to be fatigued in the late stage.That is, the process of muscle fatigue and individual MU fatigue are not completely synchronized.
(2) The firing frequency of MU (6 MUs in total) does not decrease as muscle fatigue progressed, such as MU8, MU9, MU10, MU14, MU17, and MU19.Many scholars have studied the effect of muscle fatigue on MUs; for example, [29] and [30] mentioned that after fatigue, an MU is either enhanced by increasing the central frequency to drive or the other MUs are activated.However, these conclusions are based on the whole muscle and do not study the changes in the same MU in different states.Because this paper only focuses on the change in amplitude and frequency of the same MU, the fatigue state of a single MU during muscle fatigue is probably different from the whole muscle fatigue state.In other words, whole muscle fatigue, but for an individual MU, is not necessarily synchronous with the whole muscle fatigue state.There are two possible reasons for non-synchronization: ① MU recruitment regional rotation [31], In this case, rotations that occur within the pool of MUs that receive co-input will be triggered by intrinsic changes in MU excitability [32], without violating the size principle (see [33], P 1730).As shown in Figure 10, the change of MU recruitment region can be seen.It indicates   that the number of MU recruitment and the specific MUs are changing, which is likely to affect the firing frequency of MU.  ② because an MU is either enhanced by increasing the firing frequency to drive or the other MUs are activated.An increase in the firing frequency is an increase in the neural drive in response to muscle fatigue.Therefore, the firing frequency of individual MU increase can be observed during muscle fatigue.Therefore, the firing frequency of MU does not decrease with the progress of muscle fatigue and does not synchronize with the whole muscle fatigue.
In this paper, the surface IZ location and FWHM of MUs are used to determine the three-dimensional spatial location of MUs, and MUs are then separated according to their threedimensional spatial locations.In practice, the amplitude and waveform similarity of MUAPs are also taken into account to distinguish MUs that overlap in three-dimensional space.The distribution and structure of muscle fibers is complex, a large variation exists in the distribution of fibers within the territories of biceps brachii MUs [34], [35], [36], and different MUs may also have different structures.Those MU variations will lead to variations of their MUAP waveforms, therefore MUAP waveforms which convey MU-specific features, were employed in this study to distinguish MUs, jointly with 3D location information of MUs.It is rare if not impossible that two MUs overlap exactly in the three-dimensional space, and also share exactly same size and structure, then our method may fail in distinguishing them.
The current study presents the first effort to investigate the characteristic alterations of the same MU before / after fatigue, as well as during continuous fatigue, based on HD-sEMG recordings.Despite the promising results offered by this approach, it does have some limitations.The first is that the efficacy of this approach highly depends on the decomposition results.As we can see from Figures 2 and 6, some MUs cannot be decomposed from some contraction trials.In the current study, the contraction force was focused on low contraction levels of 10%-20% MVCs.At higher contraction levels, the difficulty of detecting the same MU in different contraction trials may increase due to signal superimposition and cancellation.
FWHM is essentially a continuous variation of amplitude value of MUAP at the electrode.Although MU depth and MU size are related to MUAP, MU size has much less effect on FWHM compared to MU depth, which has a high correlation with FWHM [22], [23], [37], [38].As such, FWHM has been used in MU depth estimation in previous studies [39], [40].A question was raised in [41], how similar should the surface representation be for paired action potentials to be classified as belonging to the same motor unit?This study hypothesized that two MUs with the same spatial location and high similarity of MUAP waveform are more likely to be the same MU.To test this hypothesis, the same MU before and after fatigue was compared in this study, and the CC value of the same MU before and after fatigue was about 0.91%.There are two reasons leading to the high CC value of the same MU before and after fatigue: 1) In these experiments, the muscle morphology can be considered stable under the conditions of constant force and isometric contractions, so the shape of motion unit action potential is also stable [24], [25], [41]; 2) when obtaining MUAP in this paper, STA technology is used, which weakens the distortion of the waveform.It is concluded that the waveform similarity of the same MU is generally higher than that of different MU.However, two MUs with high waveform similarity are not necessarily the same MU.If the waveform similarity is high, the possibility of the same MU is high.
One limitation of this study is, there are other factors which may affect UMAP but have not been fully taken into account, such as conduction velocity, MU size, recruitment principle, distribution of muscle fibers etc.A future study is needed to address this limitation.

V. CONCLUSION
In the present study, the muscle fatigue mechanism was investigated at the MU level based on HD-sEMG recordings.By taking advantage of the sEMG decomposition technique, the same MU was identified in different contraction trials, and the similarity in MUAP morphology was calculated.The current results demonstrate that MUAP morphology maintained good stability which greatly enhances our understanding pertaining to the mechanism of muscle fatigue.Behaviors of MU frequency and amplitude vary largely among MUs, reflecting possible compensation mechanisms among MUs in the process of muscle fatigue.

Fig. 1 .
Fig. 1.Algorithm block diagram and an example of identifying the same MU.(A) The algorithm flow of this paper.(B) An example of a one-channel sEMG signal.(C) IPT, which is decomposed by the KmCKC method.(D) The bipolar waveform of the MUAP, with the arrow indicating the direction of the muscle fibers and the " * " indicating the surface position of the MU.A waveform on a row of electrodes in the direction perpendicular to muscle fibers and is used to calculate the arc length.(E) MUAPs from a row of electrodes perpendicular to the direction of the muscle fibers.

Fig. 5 .
Fig. 5. MF and MPF of the fatigue test in 5 subjects.The red line corresponds to the slope of the MPF.The blue line corresponds to the slope of the MF.

Fig. 6 .
Fig. 6.Comparison of amplitude (top) and frequency (bottom) of the same MU during fatigue (20%MVC).T1-T5 correspond to specific time windows of the fatigue trial, and each time window is 10s long.T1 and T5 are the first 10 seconds and the last 10 seconds of the signal, and T1, T2, T3, T4, T5 have equal time intervals.

Fig. 9 .
Fig. 9. Part of the MU firing time, the time period is T1 in Figure 6.