Muscle fatigue analysis during dynamic contractions based on biomechanical features and Permutation Entropy

1.   Introduction

Muscle fatigue is a reduction in the maximum voluntary force induced by exercise ^[1]. This topic has been addressed by different exercise physiology researchers, as it has a great impact on athlete performance, and is an important risk factor for the occurrence of overload injuries that, in some cases, can be irreversible ^[2,3]. It has been shown that adequate physical conditioning can help avoid fatigue ^[4]. Such conditioning involves modifying or avoiding exercises or tasks that cause fatigue, in turn preventing muscle injuries. Muscle fatigue causes numerous injuries in workers who perform mechanical or strenuous tasks, generally associated with overload efforts. For example, in Colombia, 38.04% of the population declares to be working in conditions of fatigue most of the time ^[5], and according to occupational disease data, between 2009 and 2012, an average 88.25% of the reported accidents were associated with musculoskeletal injuries. Surface electromyography (sEMG) has been used for muscle fatigue analysis and many researchers have studied the changes during prolonged exercise. In this regard, there are some important topics to address, such as contraction type, fatigue states and sEMG features.

The study of muscle fatigue has been addressed mainly for isometric contractions ^[6,7,8,9], where the length of the muscle does not change during exercise. However, dynamic contractions (where the length of the muscle changes during exercise), represent a more realistic view of many daily activities as well as sports (i.e. walk, carry objects, heavy lifting, run, cycle, among others), but these contractions have not been explored in depth. Dynamic contractions are more challenging; first, the muscle movement adds low frequency noise and artifacts to the sEMG, as well as hindering its acquisition and processing. Second, the activation of sEMG is not continuous, which is why it is necessary in many cases to segment the signal. Finally, the performance of sEMG features decrease in this signal.

Most of the work in this field approaches muscle fatigue as a two state process: Fatigue and Non-Fatigue ^[10,11]. However, some authors propose to address this problem in three states (Non-Fatigue, Transition-to-Fatigue, and Fatigue), because to introduce the transient-to-fatigue state before the onset of fatigue allows the development of different prediction and/or detection algorithms of muscle fatigue. In ^[8] they used this approach to detect the Transition-to-Fatigue state, and based on it made predictions regarding how long it takes the bicep muscles to reach the Fatigue state during isometric contractions. Other studies have used sEMG and deep learning to detect muscle fatigue states ^[12]. Comparisons between classifiers has also been conducted, showing that Auto-Regressive Models has the best performance for detecting the three states of muscle fatigue, while Support Vector Machines, are more useful for Fatigue detection ^[13].

Many studies have reported on the performance of different indexes in time domain, spectral domain or time-frequency domain ^[9]. The spectral domain indexes have been specially used in studies working with dynamic contractions ^[14,15,16]; however, of late, time-frequency indexes have become a trend which seems to be more effective ^[17,18]. Generally many studies relate several indexes like Mean Frequency (FM), Median Frequency (Fmed), Dimitrov's spectral index (Flnsmk), the Energy of the Detail Coefficients (EDC$_k$) of Wavelet transform, and even some time domain indexes like Zero Crossing (ZC), Root Mean Square (RMS), Waveform Lenght (WL), and so forth ^{[2,13,14,19,20,21]}. Non-linear parameters, which are not as commonly used, can also be implemented in the analysis of muscle fatigue ^[9]. The use of Non-linear paremeters was proposed by Nieminen and Takala (1996) ^[9] who suggested that sEMG signals could be better modeled as outputs of a nonlinear dynamic system, which opens up more possibilities in the assessment of muscle fatigue and new techniques like entropy analysis, fractal analysis and recurrence quantification analysis. Permutation Entropy (PE) is a recently proposed non-linear parameter that has been proved in many biomedical applications ^[11,22]; however, there are limited studies assessing its performance in sEMG signals during conditions of muscle fatigue. PE is interesting to consider in this application, because unlike other non-linear features, PE has low computational complexity and is almost parameter free.

The aim of this paper is to asses the performance of PE in distinguishing between three muscle states (Non-Fatigue, Transition-to-Fatigue and Fatigue) during dynamic contractions using sEMG and biomechanical signals.

2.   Methods

This study was divided into 5 stages, described in Figure 1. First, a muscle fatigue induction protocol was established in order to acquire the data. Second, the signals were segmented to obtain each muscle contraction. Third, the signals were characterized, including time and frequency features and PE. Then, clustering based on biomechanical features was used to label the signals of each contraction into three classes: Non-Fatigue, Transition-to-Fatigue, and Fatigue. Finally, we assessed the feature performance in distinguishing among the three muscle states.

Figure 1.  Methodological structure.

DownLoad: Full-Size Img PowerPoint

2.1.   Fatigue protocol and data acquisition

Nine healthy young adults (5 women and 4 men) aged 21 +/- 3, who do not frequently perform physical activity participated in this study. All subjects with injuries in upper limbs and cardiovascular or neuromuscular pathologies were excluded. This study was approved by the Ethics Committee for Scientific Research of the Instituto Tecnologico Metropolitano of Medellin. All the participants signed an informed consent which explained the study and the risks of their participation.

This study focuses on muscle fatigue induced in the deltoid muscle, as the shoulder is considered a critical joint in musculoskeletal injuries, therefore the study of the deltoid muscle group has been useful in implementing algorithms for muscle fatigue detection during dynamic contractions ^[14,23]. For that reason, we propose a lateral dumbbell lifting exercise that produces a right or obtuse angle between the arm and the trunk. During this exercise, the back is held straight without any lateral movement of the trunk. In this protocol, the participants were seated during the exercise and posture was visually controlled by an instructor.

The load during the exercise, was determined by the Maximum Repetition Value (MRV) of each participant and was computed by the Brzycki expression ^[24] as:

where $w$ is the weight of the dumbbell and $C_n$ is the number of contractions. Each individual warmed-up and was given an adjustable dumbbell (21.15 +/- 2.25 lb for men, and 10.2 +/- 1.2 lb for women). At that point, they were asked to perform as many lateral dumbbell lifts as they could with each arm (Figure 2)

Figure 2.  Threshold based on the 25th percentile.

DownLoad: Full-Size Img PowerPoint

Once the MRV was determined, the participants performed a first round of the exercise, at 60% of the MRV with each arm as long as they could or until they started to express fatigue. After a 30-minute break, they repeated the routine with 50% of the MRV. A total of 4 tests per participant were conducted. We used a 25-beats-per-minute metronome to control the cadence of the contractions, and participants were asked to perform as many contractions as possible with the beat of the metronome.

PowerLab 16/35 (AD Instruments Inc.) equipment with a sampling frequency of 2KHz and a Notch filter for deleting the line interference was used for acquiring sEMG signals. Subsequently, a butterworth pass band filter of second order between 30-300 Hz was applied, based on the methodology proposed in ^[25] in order to improve burst detection. PowerLab was connected to the participant with Ag/AgCl surface electrodes of a 20 mm diameter disc, placed on the half zone of the lateral deltoid with a 20 mm inter-electrode distance and with a reference electrode in the back at the height of the C7. The study area was previously cleaned with alcohol and was not shaved.

Each test was recorded with Basler high-speed cameras to capture the biomechanical signals. Reflective markers were placed on the ribs, forearm, and Acromion of the participants (Figure 2). Kinovea software was used to track the markers and obtain their position signals.

2.2.   Contraction segmentation

sEMG segmentation is an important task in the analysis of dynamic contractions as the muscle activation is not continuous, unlike isometric contractions where a simple moving window is enough to analyze fatigue progression. However, in dynamic contractions it is necessary to obtain the segments of muscle activation to avoid rest segments, where muscles fibers are not activated (sEMG does not carry any information), enabling us to analyze the fatigue progression based only on the activation segments. Therefore, it was necessary to implement a preprocessing stage to segment the sEMG and biomechanical signals (angle, angular speed, and angular acceleration).

2.2.1.   sEMG signals

The threshold-based method proposed in ^[25] was used for contraction (burst) segmentation. According to their method, initially, a passband filter between 30 and 300 Hz was applied. Subsequently, we used the Teager-Kaiser (TKEO) energy coefficient ($\psi$) defined as 2.2:

where $X$ is a sEMG signal of length $N$, defined as $X = \{x_1, x_2, \cdots, x_N\}$. Afterwards, the signals were rectified and filtered with a 50-Hz low-pass filter. Finally, burst detection was carried out based on the Approximate Generalized Likelihood Ratio (AGLR) method ^[26,27,28]. AGLR is based on statistical testing of the null hypotheses $H_0$ (sEMG signal is not a burst) and the alternate hypothesis $H_1$ (sEMG signal is a burst), describing the statistical properties of a series of sEMG samples and uses a thresholding approach defined as ^[25]:

where $p_1$ and $p_0$, represent the probability that sEMG segments correspond to a contraction (burst) or not, respectively.

Algorithm 1, describes the pseudocode for the segmentation of a sEMG signal. Parameters $\mu_0$, $\sigma_0$, $\mu_1$ and $\sigma_1$, refer to the mean and deviation of univariate Gaussian Probability Density Function, which were obtained after manually segmenting 6 bursts and 6 baseline signals of each subject. Parameters $k$ (window length) and Th (Threshold) of the AGLR method where set to 300 and 15 respectively after a heuristic search (results of this task are not shown in this publication, as it is not relevant to the main objective).

After segmentation, we removed the border of the sEMG segments in order to ensure that all segments had the same length (1000 points, equivalent to 0.5 seconds of signal) and to avoid border effects related to the beginning and the end of the contraction.

2.2.2.   Biomechanical signals

The reflective markers were placed on each participant forming a triangle between marker 1 on the forearm ($M1$), marker 2 on the acromion ($M2$) and marker 3 on the ribs ($M3$), (Figure 3), in order to obtain the angle between the arm and the trunk during exercise.

Figure 3.  Angle calculation during exercise. M1: ulnar fossa. M2: acromion. M3: trunk.

DownLoad: Full-Size Img PowerPoint

Kinovea software was used to track the positions of the three markers. Kinovea returns the position of the marker as Cartesian coordinates (based on pixels) in every frame of the video (500 frames per second). Then, we computed the distance $d1$ (distance between forearm and ribs), $d2$ (distance between forearm and acromion) and $d3$ (distance between acromion and ribs) in pixels using the Euclidean norm. Distances $d2$ and $d3$ are constant, representing the anthropometric measurements of any subject. However, these distances presented small changes between frames, due to tracking errors and involuntary subject movement, for this reason, we set $d2$ and $d3$ as the average distance during all frames.

Being $\theta_n$ the angle signal formed between the trunk and the forearm, it was computed using law of cosines as:

where $d1_n$ is the distance between forearm and ribs, in the moment $n$.

Biomechanical signals angular velocity ($\omega_n$) and angular acceleration ($\alpha_n$) were computed as the first and second derivative of $\theta_n$. Then, the signals were filtered using third order Butterworth low-pass filters at cut-off frequencies of 100 Hz (for $\theta_n$ and $\omega_n$) and 30 Hz (for $\alpha_n$).

Algorithm 2 shows the process for computing biomechanical signals $\theta_n$, $\omega_n$ and $\alpha_n$. Inputs M1, M2 and M3 are the matrix of dimensions $N \times 2$, that represent the Cartesian position at $N$ time of the markers according to the figure 3.

Finally, every contraction in the signals ($\theta$, $\omega$ and $\alpha$) was segmented. For this purpose, we employed a segmentation approach based on a threshold for the $\theta$ signal, and used the same times to construct a remaining signal, as follows:

when $\mathrm{Mask}_n = 1$ refers to a contraction, in this work we used the $25^{th}$percentile of $\theta$ as the threshold ($Th$) (Figure 4).

Figure 4.  Threshold based on the 25th percentile.

DownLoad: Full-Size Img PowerPoint

2.3.   Fatigue progression labeling based on clustering and biomechanical features

Assuming that biomechanical changes are due to muscle fatigue progression, we propose to obtain features from $\theta$, $\omega$ and $\alpha$ as described in Table 1.

Table 1.  Computed biomechanical features.

Feature	Abbreviation	Formula
Maximum angle	Max$_\theta$	$\max \theta_n$
Mean angle	Mean$_\theta$	$\frac{1}{N}\sum_{n=1}^N \theta_n$
Integral of angle	Int$_\theta$	$\sum_{n=1}^N \theta_n$
Maximum velocity	Max$_\omega$	$\max \omega_n$
Minimum velocity	Min$_\omega$	$\min \omega_n$
Mean velocity	Mean$_\omega$	$\frac{1}{N}\sum_{n=1}^N \omega_n$
Integral of velocity	Int$_\omega$	$\sum_{n=1}^N \omega_n$
Mean speed	Mean$_{|\omega|}$	$\frac{1}{N}\sum_{n=1}^N|\omega_|$
Maximum acceleration	Max$_\alpha$	$\max \alpha_n$
Minimum acceleration	Min$_\alpha$	$\min \alpha_n$
Mean acceleration	Mean$_\alpha$	$\frac{1}{N}\sum_{n=1}^N \alpha_n$
Integral of acceleration	Int$_\alpha$	$\sum_{n=1}^N \alpha_n$

---

 | Show Table

DownLoad: CSV

Each of the biomechanical features obtained was normalized intra-subject, as some biomechanical variables such as Max$_\theta$ and Max$_\omega$ present great variability between subjects. Next, the features were correlated with the contraction number in order to identify biomechanical features related to muscle fatigue progression. In this study, values with an absolute correlation greater or equal to 70% were chosen as relevant biomechanical features and were used for clustering.

We employed an unsupervised learning algorithm in order to label the contraction into three groups, which has a direct relation with changes in the biomechanical variables of the subject, for practical effects we refer to these different stages as: Non-Fatigue, Transition-to-Fatigue, and Fatigue, although this does not involve a clinical diagnosis of muscle fatigue. Most of the literature in this topic suggests that muscle fatigue has two stages (Non-Fatigue and Fatigue) ^[11,29]; however, other works show that using a three stage approach allows for the development of a prediction system able to detect fatigue before its clinical appearance ^[30]. Therefore, we implemented a type of agglomerative hierarchical clustering based on Ward's method. The clustering procedure was conducted using the contractions obtained from all the tests. In agglomerative clustering, initially every point (in this application, a point is a contraction) in the data set is an individual cluster and clusters are progressively merged according to a similarity criterion. Ward's method states that a cluster should be fused with the cluster that causes the least increase in within-cluster sum of squares ^[31]. Formally, the within-cluster sum of squares is defined as the sum of the squares distance between points in the cluster and the centroid of the cluster. The resulting sum of squares that joins clusters $a$ and $b$ can be computed as ^[32]:

where $n_a$ and $n_b$ are the number of elements in the clusters, $\overline{\bf{x}}_a$, and $\overline{\bf{x}}_b$ are the cluster centroids.

2.4.   Muscle fatigue analysis

Once the fragments of the sEMG signal were obtained and labeled as described in Section 2.3, we computed different sEMG features PE, RMS, ZC, FM and Fmed to confirm that the changes in biomechanical variables were due to muscle fatigue progression. Then, we evaluated if there was a statistically-significant difference between the sEMG features of the three groups. For this purpose, we used an Analysis of Variance (ANOVA) to establish the existence of a statistical difference between the 3 groups (Non-Fatigue, Transition-to-Fatigue, and Fatigue). Subsequently, we performed a post-hoc analysis using the Less Significant Difference Test (LSD) to directly compare pairs of groups, i.e. Non-Fatigue vs Transition-to-Fatigue, Transition-to-Fatigue vs. Fatigue, and Non-Fatigue vs. Fatigue.

Finally, we carried out a Receive Operator Curve (ROC) analysis, which uses logistic regression as classifier, in order to assess the discriminant capacity of the more promising sEMG features to distinguish between Non-fatigue, Transition-to-Fatigue and Fatigue.

2.4.1.   Permutation Entropy

Entropy measures are very common for quantifying the complexity of a time series; however, common approaches, such as Shannon entropy, neglect to take into account effects from the temporal order of the values in the time series. Bandt and Pompe proposed in ^[33], a natural encoding which reflects the rank order of successive ($x_n$) in sequences of length ($N$). In this regard, given an sEMG signal $X = \{ x_1, x_2, \cdots x_N \}$, if we take signal segments of $m$ points (Embedded dimension) with a time delay of $\tau$ as $S_i = \{X_i, X_{i+\tau}, X_{i+2\tau}, \cdots X_{i+\tau(m-1)} \}$. Sorting $S_i$ in ascending order as $X_{i+r_0} < X_{i+r_1} < \cdots < X_{i+r_{(m-1)}}$, we can obtain its ordinal pattern as $\pi_j = \{r_0, r_1, \cdots r_{(m-1)} \}$ ^[34]. Then, the set of ordinal patterns is $\Pi = \{\pi_1, \pi_2, \cdots \pi_{m!} \}$. The permutation entropy (PE) is computed using the Shannon form, but taking the relative frequency of the ordinal patterns as:

where ($p'(\pi_j)$) represents the relative frequencies of the ordinal pattern $\pi_j$. PE essentially measures information based on the occurrence or absence of certain permutation patterns of length $m$, in the ranks of values of a time series ^[35]. A detailed explanation of PE calculation can be found in ^[35].

In this study, after computation, PE values were normalized intra-subject with the aim of reducing the impact of between subject differences, and we conducted an exhaustive search of the embedded dimension $m$ (between 2 and 10) and the time delay $\tau$ (between 1 and 10).

2.4.2.   Other muscle fatigue features

We analyzed features in the time domain, such as the Root Mean Square (RMS) which is commonly used in studies on muscle fatigue ^[9,14,21]. We applied fast Fourier transform to sEMG signals and computed features such as FM and Fmed which are considered traditional indicators of muscle fatigue ^[9,16,36]) in each of the Burst obtained. These features are very common in muscle fatigue analysis, and some studies have shown that during static contractions the FM usually decreases ^[9,15], and Fmed has been shown to have a similar behavior ^[14,29]. All the features described were computed according to table 2 and were normalized intra-subject to reduce the impact of between subject difference.

Table 2.  Features computed from sEMG signals.

Feature	Equation
Permutation Entropy (PE)	$\sum_{j=1}^{m!} - p'(\pi_j) \log_{2}p'(\pi_j)$
Root Mean Square (RMS)	$\sqrt[]{\frac{1}{N}\sum_{n=1}^{N} |x_n^{2}|}$
Mean Frecuency (MF)	$\frac{\sum_{f=F_{min}}^{F_{max}} f S[f]}{\sum_{f=F_{min}}^{F_{max}} S[f]}$
Median Frecuency (Fmed)	$\sum_{f=F_{min}}^{F_{med}}S[f] = \sum_{f=F_{med}}^{F_{max}} S[f]$
$p'(\pi_j)$ is the relative frequency of the ordinal pattern $\pi_j$ $X$ is a sEMG signal as $X=\{x_1, x_2, \cdots x_N \}$ $S[f]$ is the power spectrum density of a signal at frequency $f$

---

 | Show Table

DownLoad: CSV

3.   Results and discussion

3.1.   Fatigue protocol and data acquisition

The data for this study was obtained from only nine subjects; however, this does not seem to be a disadvantage, as several studies have used a similar number of subjects ^[37,38,39] or less ^{[12,21,40,41,42]}. This number of subjects is feasible as several samples are recorded per subject; for example, in isometric exercises long signals can be recorded ^[42] and it is possible to apply windows of time in order to increase the number of samples per subject ^[12]. On the other hand, muscle activity during dynamic contractions can be analyzed by means of a contraction set, as these signals change greatly over time, like in ^[14], where the authors did an analysis of different parameters based on their behaviour in each contraction of a given exercise.

In this study, the nine subjects performed 2 tests, first, with 60 % of the MRV and then with 50 %. This process was repeated for each arm, obtaining 4 signals per subject, for a total of 36 signals. Additionally, signals were treated as a set of contractions, thus obtaining a collection of 950 samples that include sEMG and position signals.

3.2.   Contraction segmentation

3.2.1.   sEMG

The segmentation process is described step by step in Figure 5 and Algorithm 1. (a) presents the original signal obtained in one of the tests. (b) shows the signal after applying the 30–300Hz bandpass filter. (c) is the signal after applying the TKEO. (d) plots the signal of TKEO after passing through a 50-Hz low pass filter and rectification. Finally, (e) describes the way the signal is segmented using the approximate generalized likelihood ratio (AGLR) method.

Figure 5.  sEMG segmentation procedure.

DownLoad: Full-Size Img PowerPoint

The AGLR method using the TKEO seems to be an adequate method for sEMG signal segmentation. In another study ^[25], the inclusion of the TKEO in signal conditioning significantly reduced the detection error of muscle activation, which is reflected in the precision of fragmentation of the sEMG burst. However, in this study, intra-subject preprocessing was necessary since a general model presented significant errors in burst detection. By adopting this method, we were able to segment all the contractions in sEMG signals.

3.2.2.   Biomechanical signals

Figure 6 shows the filtering process of biomechanical signals. Noise can be observed due to the marker tracking process in Kinovea and the involuntary movements of the subject. The applied filters present a suitable correction of the signal under study.

Figure 6.  Pre-processing of biomechanical signals.

DownLoad: Full-Size Img PowerPoint

The proposed segmentation method is based on thresholding using the $25^{th}$ percentile of $\theta$ signal, which has shown to be effective because only one of the 36 tests had to be manually segmented. Since the signals are equal in time, the same mask was applied to the three of them. In Figure 6, we can observe the preprocessing stage applied to the biomechanical signals. (a), (b), and (c) show the signals obtained from the motion analysis system using the methods described in Section 2.3.2. In turn, (d), (e), and (f) are filtered biomechanical signals with low pass filters. Finally, (g), (h), and (i) present the resulting segmentation applied to biomechanical signals.

As a result of the segmentation process, we obtained 944 contractions associated with sEMG bursts with an average duration of 1.34 s (2683 points approximately); however, they were trimmed to a duration of 0.5 seconds, to avoid border effects. Finally, the same number of biomechanical signals ($\theta$, $\omega$ and $\alpha$) were obtained from 36 trials.

3.3.   Fatigue progression labeling based on clustering and biomechanical features

Several biomechanical features were calculated in this study, in order to identify which of them are related to fatigue progression. We identified features with a Pearson's correlation greater than or equal to 70%. This correlation analysis is summarized in Table 3. The results show that all features obtained from $\theta$ are related to fatigue progression, while $Max_\omega$ is the only feature related to angular velocity. On the other hand, angular acceleration features were not related to fatigue progression.

Table 3.  Correlation of biomechanical features with Exercise Progression.

Biomechanical features	Correlation
Max$_\theta$	0.82
Int$_\theta$	0.81
Mean$_\theta$	0.81
Max$_\omega$	0.7
Mean$_{|\omega|}$	0.68
Min$_\omega$	0.55
Int$_\alpha$	0.46
Max$_\alpha$	0.31
Mean$_\alpha$	0.3
Min$_\alpha$	0.23
Int$_\omega$	0.22
Mean$_\omega$	0.15

---

 | Show Table

DownLoad: CSV

Figure 7 shows the results of the clustering process and the evolution of states within the exercise progression. Furthermore, the clustering process identified the 3 following groups:

Figure 7.  Results of agglomerative hierarchical clustering based on biomechanical indices.

DownLoad: Full-Size Img PowerPoint

Group 1 (Non-Fatigue): these are the contractions at the beginning of the test, indicating that muscles are in an optimal state during the contractions and the biomechanical conditions have not been compromised (Non-Fatigue).

Group 2 (Transition-to-Fatigue): these are the contraction in the middle of Group 1 and Group 3, we assume this group of contractions corresponds to a slight change in the biomechanical features of the subject. This stage was called Transition-to-Fatigue.

Group 3 (Fatigue): these are the contractions at the end of the exercise. During these contractions the subject can not continue the exercise and there is a big change in the biomechanical features, indicating the possible presence of muscle fatigue.

In other hand, Figure 7 (b) shows the results of the correct labeling of the signals, while Figure 7 (c) highlights some clustering errors. We identified that 72.3 % of the signals have some mislabeled contractions; however, this represents just 8.06 % of the contractions.

It is important to state that this labeling is not proof of clinical fatigue, but shows a significant variation in the subject's biomechanical variables, which can lead to possible musculoskeletal injuries during sport or work-related movements. Although this labeling process can not confirm the presence of clinical fatigue in the muscle, it is still relevant as it can identify changes in the biomechanical variables, possibly preventing injuries.

Such variability in biomechanical features is difficult to quantify between subjects ^[43], which means that every subject has their own range of motion according to their age and physical condition, among other factors. For this reason, we evaluated the changes in biomechanical features after intra-subject normalization. Figure 8 shows the distribution of 2 biomechanical features $Max_\theta$ & $Max_\omega$) after intra-subject normalization. From both box-plots it is clear that biomechanical features decrease with fatigue progression, however $Max_\omega$ shows a larger variability than $Max_\theta$.

Figure 8.  Biomechanical features distribution among the 3 different fatigue states.

DownLoad: Full-Size Img PowerPoint

Several researches have used strategies for labeling muscle fatigue stages, based on biomechanical or kinematic features, but the majority have not demonstrated it clinically, as it is an invasive and expensive procedure. For instance, some authors have assumed a fatigue stage after a set of intensive exercises ^[44,45] or have carried out the same activity until the participants were exhausted ^[46,47,48]. Others have implemented a threshold when a biomechanical variable such as force, velocity, or distance has decreased enough to assume muscle fatigue ^[6,49]. Our labeling method is more reliable than establishing a threshold for a biomechanical variable or than assuming muscle fatigue after an exercise set, as it not make any assumptions but rather groups the contractions according to their biomechanical similarities.

3.4.   Analysis of sEMG during the exercise prolongation

In order to confirm that the biomechanical changes used to labeled the contractions were due to muscle fatigue progression, we computed sEMG features and compared them between the three muscle fatigue states. Figure 9 shows a Heat-map of the ANOVA test p-values between Non-Fatigue, Transition-to-Fatigue and Fatigue for PE using different values of $m$ and $\tau$ (Embedded dimension and time delay). The greatest statistical difference was obtained with $m = 7$ and $\tau = 6$. Taking into account that the sEMG length was of 1.34 s (i.e. $N = 2683$) and the recommendation $N > > m!$ stated in ^[33], we found that with $m = 7$ ($m! = 5040$) the number of possible ordinal patterns is larger than the number of signal points, indicating that the recommendation of $N > > m!$ is not fulfilled and therefore, the obtained PE values are not a good approximation of sEMG complexity. However, these results confirm that for classification purposes, it is not necessary for this rule to be satisfied as was demonstrated in ^[34].

Figure 9.  Heat-map of p-values for PE at different embedded dimension ($m$) and time delay ($\tau$).

DownLoad: Full-Size Img PowerPoint

On the other hand, Figure 10 shows the behaviour of the PE in different subjects. In this figure the mean value of PE from all the tests applied to each subject was extracted in each of the assessed states (Non Fatigue, Transition-to-Fatigue and Fatigue), and they were superimposed in order to evidence inter-subject trends. It can be observed thaty in all the subjects the PE value decreases progressively, in some subjects with a higher slope than others. The mean PE value (intra-subject normalized) for all subjects was 0.75 in Non-Fatigue state, 0.54 for Transition-to-Fatigue and 0.37 for Fatigue. The PE value decreased around 50.6% from Non-fatigue to Fatigue contractions, indicating a clear reduction in sEMG complexity in the presence of fatigue.

Figure 10.  Behavior of PE during the exercise prolongation in different subjects.

DownLoad: Full-Size Img PowerPoint

Table 4 shows a comparison of the ANOVA test results of each sEMG feature and LSD between the groups defined in Section 3.2, considering p << 0.01 as statically significant difference. Results show that biomechanical changes are due to the progression in muscle fatigue and the retrieved groups from the clustering correspond to different states of the muscle during exercise, as the statistical test revealed that all the sEMG features exhibited changes in at least 2 of the 3 comparisons. Additionally, these findings are comparable to other muscle fatigue studies during dynamic contractions ^[9,14] where the Fmed and ZC shows statistically significant differences between stages of rest and after exercises, with decreasing behavior. This result is contrary to ^[50] who reported a increment of Fmed during the fatiguing exercise.

Table 4.  Statistical comparison of sEMG features between three labeled contractions.

Feauture	ANOVA test	LSD
		G1 vs G2	G1 vs G3	G2 vs G3
PE	p$< <$0.01	p$< <$0.01	p$< <$0.01	p$< <$0.01
Fmed	p$< <$0.01	p$< <$0.01	p$< <$0.01	p$< <$0.01
FM	p$< <$0.01	0.0074	p$< <$0.01	0.25
RMS	p$< <$0.01	0.4856	0.0073	p$< <$0.01

---

 | Show Table

DownLoad: CSV

A ROC analysis was applied to the more promising features (PE and Fmed) in order to observe how they distinguish between the different classes, the Area Under Receive Operator Curve (AUROC) help us to analyze the performance of our algorithm in the classification task. As shown in Figure 11, the performance of PE is slightly higher than Fmed, although both features have a high performance for distinguish between Non Fatigue and Fatigue states. This result was expected as several studies have previously found significant changes between these stages ^[9,51,52].The worst performance was obtained when the algorithm distinguished between Transition-to-Fatigue and Fatigue, possibly due to the fact that the biomechanical changes were not so important in this case; however, we observed that PE performs better than Fmed in these two stages and both are greater than FM (AUROC = 0, 54). The use of both features (PE and Fmed) as inputs does not improve the classification performance between the groups. 11 sEMG muscle fatigue features and 6 mathematical operators

Figure 11.  Sensitivity/specificity report using ROC curve analysis.

DownLoad: Full-Size Img PowerPoint

Authors like ^[14,29] describe the behavior of some frequency and time domain features. We found correspondence in the Fmed feature because (as can be observed from Table 4) it is highly variable with respect to the three groups of muscle fatigue. However, FM and RMS fail to identify the 3 states of muscular fatigue and only differentiate between Non-Fatigue and Fatigue. Another study presents the behavior of PE in the evaluation of muscle fatigue, finding that the PE value decreases with exercise prolongation ^[11]. We obtained similar results; however, we demonstrated a constant change of PE during exercise, while ^[11] only considered the first and last of the six equal segments using time scales. Therefore, this could indicate that PE is the most adequate feature for assessing muscle fatigue in the three states.

Other authors have analyzed the behavior of several biomechanical variables in terms of the evolution of a repetitive task, finding that biomechanical variables can adapt with muscular fatigue progression ^[44,53]. In addition, another works have demonstrated how the muscle fatigue influences directly the biomechanical and kinematic behavior of some joints ^[49,54]. Similar results have been found in this work because changes in variables like angle and velocity enabled to identify the three states of muscle fatigue progression. Thus, these results suggest a relation between the subject's biomechanical variables and muscle fatigue progression.

The classification accuracy of PE between transitions between states can be seen in Table 5, as well as in other reports in the literature. Initially, our method did not perform as well as other approaches, however, there are several reason behind this. Firstly, the highest performances are obtained for isometric contractions, while our method works in dynamics contractions, the latter which are more challenging to analyze due to non-stationary signal and the movement artifacts ^[9]. Secondly, most studies only distinguish between Non-Fatigue and Fatigue, while our study aims to classify between the three states. Thirdly, most muscle fatigue analyses has been conducted in biceps muscle, while our approach is with the deltoid muscle. Finally, it is important to clarify that our method uses a single feature to classify the signals, this means that logistic regression classifier is acting as a simple thresholding approach. In this regard, other methods such SVM or k-NN represent a more robust framework for classifying the signals, but their use is justified when one is using multiple features.

Table 5.  Comparison with other classification methods.

States	Reference	Contraction Type	Muscle	Classification method	Features	Accuracy
Non-Fatigue vs Transition-to-Fatigue	[30]	Isometric	Biceps	LDA	1D spectro-composite feature	90.37%
	[12]	Isometric	Biceps	deep belief networks	Raw Signal	66.32%
	In this Work	Dynamic	Deltoid	Logistic Regression	PE	71.0%
Transition-to-Fatigue vs Fatigue	In this Work	Dynamic	Deltoid	Logistic Regression	PE	68.0%
Non-Fatigue vs Fatigue	[12]	Isometric	Biceps	DBNs	Raw Signal	97.60%
	[30]	Isometric	Biceps	Logistic Regression	1D spectro-composite feature	84.8%
	[55]	Dynamic	Biceps	k-NN	Multifractal features	88%
	[56]	Dynamic	Biceps	SVM	Renyi entropy	86.7%
	[57]	Dynamic	Biceps	Logistic Regression	Multifractal features-Infomation Gain	84%
	[57]	Dynamic	Biceps	k-NN	Multifractal features-Infomation Gain	82%
	[55]	Dynamic	Biceps	Logistic Regression	Multifractal features	82%
	[58]	Dynamic	Biceps	LDA	Pseudo-wavelet	78.45%
	In this Work	Dynamic	Deltoid	Logistic Regression	PE	80.0%

---

 | Show Table

DownLoad: CSV

PE offers a low computational complexity cost compared to other non-linear descriptors like Lyapunov exponents or Correlation dimension ^[33]. However, when the embedded dimension is large (usually above $m > 12$), the number of possible ordinal patterns (Permutation) increases dramatically and the memory consumption required to maintain the counting of those pattern is insufficient in most machines. Despite the above, for this particular application we found that an embedded dimension of 7 is optimal for characterizing the EMG signal in order to detect muscle fatigue. This means that the number of possible ordinal patterns is only 5040 and the required memory to maintain the counting is around 20.1 Kilo bytes (using an integer of 32 bits for counting), such memory capacity is available in low cost devices such as an Arduino UNO (32 KB). Regarding computational complexity, in the literature there are several proposals for computing PE in a fast and robust way, such as Piek et al. that proposed a PE algorithm based on a binary vector representation of ordinal patterns and achieved a computational complexity of $\mathcal{O}(m^2 \cdot N \cdot \log N)$ ^[59], their implementation is based on bit shift operations which are very efficient in execution time, resulting in an attractive possibility for using PE in real time application.

4.   Conclusions

This work presents the analysis of muscle fatigue progression, during dynamic contractions, based on Permutation Entropy and biomechanical features. Results reveal that biomechanical features can be used to distinguish between different stages of muscle fatigue, mainly, features derived from angle ($\theta$) and angular velocity ($\omega$). Additionally, Permutation Entropy (PE) of surface electromyography (sEMG) signal, had the best discriminant capacity between Muscle Fatigue, Transition-to-Fatigue and Non-fatigue, over-performing classical sEMG features as Mean and Median Frequency (Fmed and MF), and Root Mean Square (RMS).

In future work, it would be valuable to look into the use of PE for predicting muscle fatigue from sEMG signals.

Conflict of interest

The authors have no conflicts of interest to report.