SYNERGOS: A Multiple Muscle Activation Index

1. Introduction

An important movement control strategy used by the central nervous system (CNS) is the activation of multiple muscles acting in concert with each other to achieve a specific movement [1-6]. The specific temporal pattern of motor unit firing, the number of recruited motor units, and the intensity of the firing units result in a state of Multiple Muscle Activation (MMA) during each movement. We hypothesized that quantifying the MMA would provide information about certain states of muscle activation used by the CNS at any moment of a given movement activity. As neuromuscular disorders disrupt the firing characteristics of the motor units by which the CNS controls human movement, quantifying the degree of MMA might be useful as an assessment tool for these disorders. Additionally, monitoring changes in MMA over time might provide valuable information about the progression of a disease state or conversely, a recovery profile. Currently, no easily implemented screening technique for use in clinical settings exists that allow the changes in the MMA to be measured and tracked over time. Therefore, developing a single value that quantifies the degree of activation among multiple muscles that accounts for temporal and magnitude changes in muscle activity in a combinatorial fashion may be of value to clinicians and researchers interested in evaluating alterations in muscle activation due to different physiological and environmental constraints. In this chapter, a new index, “SYNERGOS” (from the Greek word for “working together”) is introduced to quantify the level of MMA. At its core, SYNERGOS systematically identifies the changes in muscle activity depicted by electromyography (EMG) signals obtained from multiple muscles and summarizes the coactivity among these muscles into a single scalar quantifying the MMA over a predefined period (i.e. in this report, the gait cycle) during a variety of movements.

Surface EMG is commonly used to investigate neuromuscular activation during human movements. EMG can provide insight into the relationship between the underlying activity of the CNS and the resulting muscle contractions that produce the observable movement. Such insight sheds light on the different neuromuscular activation patterns in both healthy movements and disease states [7, 8]. However, several studies identified the inherent stochastic characteristics of EMG signals and recommended applying nonlinear data analysis to investigate the underlying nonlinear features (i.e. determinism) of the signal [9-16]. Recurrence Quantification Analysis (RQA) has shown promising results to detect subtle changes in EMG signals which can be missed by the application of linear data analysis. RQA quantifies several parameters such as percentage of recurrence (% REC) and percentage of determinism (% DET) to detect such subtle changes in dynamical state of EMG signal [13, 15-18]. % DET is an indication of the underlying deterministic patterns in a dynamical system that is associated with the degree that an EMG signal “recurs,” or is predictable [15, 16, 18] as a result of increasing firing rate in recruited motor units [9, 10, 16, 18]. % DET has shown significant sensitivity to the change in the amount of loading and duration of a movement during both static and dynamic movements [13-18]. In this study, the % DET obtained from the RQA for each recorded muscle during predefined duration (i.e. gait cycles within each gait speed) is a single value used as an input for calculation of SYNERGOS index for each gait cycle.

Previous authors have proposed techniques that have demonstrated that individual muscles are not independently controlled by the CNS during movement [2, 4, 6, 19]. These methods of quantifying multiple muscle activities provide several time-variant and time-invariant parameters [2, 4] or several modes [6, 20, 21] demonstrating the level of activation of each muscle during a movement. These techniques share some features in common with SYNERGOS, namely the application of EMG to detect the multi-muscle actions during various activities and the use of sophisticated techniques to explore the underlying dynamics associated with the muscle activation controlling limb motion. While having these features in common with SYNERGOS, these techniques are designed to investigate and parameterize the underlying synergies by quantifying different modes of multiple muscle actions during body movements so that relative time-invariant contributions of different muscles can be identified. Conversely, SYNERGOS is not designed to identify multiple modes of activation synergies but rather SYNERGOS is unique in that a single index value accounting for the simultaneous activity of all muscles during a given movement, or over time in the case of cyclical movements (e.g., one index for each gait cycle), can be conveniently assessed.

In this report, we assess the potential of the SYNERGOS technique to evaluate changes in MMA between walking and running conditions that use different muscle activation strategies due to differences in the complexity of these movements [22]. We hypothesized that the SYNERGOS method can detect significant differences in quantified MMA during walking vs. running due to higher level of coactivities resulted from more frequent and intense simultaneous neuromuscular activities depicted by EMG signals. Additionally we hypothesized that the SYNERGOS method was significantly sensitive to identify the changes in MMA associated with small but incrementally increasing gait speed beginning with a slow walk and concluding with running. These hypotheses were based on previous studies indicating an elevated muscle coactivation associated with higher gait speed (modeled by faster knee extension/flexion) [23, 24]. In this study, the cyclic movement of walking and running was chosen, which by their nature, maintain their general kinematic form as speed increases but the underlying MMA is modulated in concordance with increasing speed. We believe that if SYNERGOS can detect minute changes in MMA associated with slight increases in locomotion speed it may have functional utility to evaluate MMA during a variety of movements.

2. Method

2.1. Subjects

Ten, healthy University of Houston students participated in the study at the Laboratory of Integrated Physiology (5 male and 5 female; age range: 20-25 years, Mean 22, SD 2.3 years, weight: 55-80Kg Mean 67.45, SD 12.90 Kg.). The exercise readiness of each subject was monitored by using a physical activity questionnaire (modified International Physical Activity Questionnaire) in which several wellness aspects such as cardiovascular fitness, discomfort during exercise, history of dizziness, joint problems, pregnancy, diabetes, breathing problems, and history of major surgery were questioned [25]. In addition, the subjects were required to report any history of neuromuscular disorder and lower extremity injury. No subject reported any of the aforementioned exercise readiness risks. These selection criteria were chosen to minimize the known effects of neuromuscular disorders and aging on changes in neuromuscular activation. [26, 27]. All procedures were reviewed and approved by the University of Houston Committee for the Protection of Human Subjects. The subjects were fully informed about the test protocol and provided signed consent prior to participating.

2.2. Experimental protocol

Each experiment was conducted in a single laboratory session. To alter the gait pattern from walking to running, the participants were instructed to walk at a comfortable initial speed (different between subjects) on a treadmill whose speed increased by 0.045 m/s every five strides up to a maximum duration of 180 seconds. Subjects were encouraged to stop the trial if at any point they became uncomfortable. All subjects successfully completed this protocol by transitioning their gait from walking to running. The transition speed (speed in which gait pattern was changed from walk to run) was recorded for each subject by observation to be used during analyses (see Statistical Analysis below). In addition, to verify the transition between walk and run, the vertical hip velocity and the stride time variability were evaluated. Both of these measures are sensitive to the sudden changes in gait pattern (i.e. walk to run transition). This approach confirmed the observed transitional strides [28]. The mean initial speed across all participants was a slow walk at 1.12 ± 0.13 m/s with final speed at a run of 3.30 ± 0.11 m/s.

2.3. Material and data collection

2.3.1. EMG activity

EMG signals were collected using six preamplifier bipolar active electrodes (EMG preamplifier, Type No: SX230, Biometrics Ltd., Gwent, UK) with a fixed electrode distance of 20 mm from rectus femoris (RF), tibialis anterior (TA), lateral gastrocnemius (GA), soleus (SO), vastus medialis (VM), and biceps femoris (BF) of the right lower limb using double sided tape. The electrodes were connected to a DataLINK base-unit DLK900 of the EMG acquisition system which was connected to a PC using USB cable. To achieve acceptable impedance level the skin over the location of each electrode was shaved and cleaned with alcohol swabs. EMG data were collected at 1000 Hz and passed through an amplifier with the gain set at 1000. The amplification bandwidth was 20–460 Hz (input impedance =100 MV, common mode rejection ratio >96 dB (~110dB) at 60 Hz). A zeroing reference electrode was placed above the right lateral malleolus bone and was secured by elastic wrap and tapes. There was no excessive filtration of the EMG data during collection but a digital filter was applied during data processing (see below). During the collection session, the electrodes were not removed from the subjects until data collection was completed.

2.3.2. Apparatus

Kinematic data were also collected at 200Hz using VICON motion capture device (Oxford Metrics, Oxford, UK) to identify the gait cycles by detecting each heel strike (i.e. right heel strike to heel strike). The kinematic markers were located on the hip, knee, ankle, heel, and toe. An electronic trigger was used to synchronize the EMG and kinematic data and to determine the events in which the treadmill speed increased (one trigger per gait speed increase).

2.4. Data processing

The data were processed by using a customized Matlab script (Mathworks, USA R2007a). The detailed specifications of each are provided below.

2.4.1. Recurrence Quantification Analysis (RQA)

The first step was to calculate the % DET by using RQA (RQASP program [29]; also see details in Appendix-RQA formulation). The EMG signals were band-pass filtered from 10-500 Hz [30, 31] however, no other smoothing techniques were used for the remaining signals to maximize the exposure of the raw signals to the nonlinear analysis method [16, 32]. To minimize the potential effect of any noise in the raw EMG signal on the outcome of RQA, the initial parameters (i.e. time delay and embedding dimension) were selected by following the recommended settings (i.e. embedding dimension, time lag, and radius) [29].

For each walking speed and each muscle, the EMG signal was clustered into five data bins that were defined as the epochs of recorded data between each right foot heel strike (i.e. each gait cycle). For each subject, the percentage of recurrence (% REC) and % DET of the clustered EMG signals was calculated for all muscles within each speed condition (% REC is required for the control during the selection of radius, see Appendix-RQA formulation for more details). As the RQA processes the time-delayed reconstructed space phase of the EMG signal, several parameters were defined prior to performing the analysis. The data were analyzed using an embedding dimension of m=6 based on the False Nearest Neighbor technique [33] and time delay of 5 (τ = 0.005 second) based on the Mutual Information (MI) technique [34]. In the MI technique, the first local minimum of the average mutual information is used to detect the time delay. The proximity radius (see Appendix-RQA formulation) was selected as 2 to 10 units of the rescaled “Maximum” unit of the Distance Matrix to keep the percentage of the recurring points in the recurrence plot (RP) of the signal less than 2% , as has been recommended [18].

2.4.2. Shuffled surrogates tests

While many studies have indicated a nonlinear dynamical pattern for EMG signals [10, 14, 15, 17, 18, 35] the nonlinearity of such signals was tested to justify the application of RQA in SYNERGOS [15, 17, 35-38]. A common practice to test the assumption of nonlinearity in a signal (i.e. EMG) is using surrogate data testing [15, 35-37, 39]. During this test the original EMG data are randomly shuffled using different algorithms namely time shuffling to generate random signals. It is expected that the randomization of the signal has significant effects on the nonlinear characteristics of the signal while keeping the linear characteristics of the signal unchanged which verifies the nonlinear behavior of the signal (i.e. EMG). In this study, 20 series of surrogate data using three algorithms were generated for each set of muscles per gait cycle [15, 35-38]. The approximate entropy (ApEn) and % DET of the original data were used to monitor the changes in the underlying dynamics of the EMG data after shuffling [40-42]. It was hypothesized that the value of % DET would significantly decrease while the value of ApEn would significantly increase for the shuffled data (see Appendix-Shuffled Surrogate Tests of EMG). By rejecting the null hypotheses of this testing procedure the existence of underlying nonlinear dynamics of the EMG signal could be assumed and therefore the application of RQA was justified.

2.4.3. SYNERGOS

We developed the SYNERGOS method to assess the level of MMA based on the activation of each muscle with all possible sets of the other muscles. SYNERGOS employs a two-step method for quantifying MMA. The first step is using RQA to analyze the EMG signals of each recorded muscle separately. The calculated % DET of the EMG signal obtained from each muscle serves as an input variable for the second step of SYNERGOS in which the inputs are combined by using a novel method that quantifies the level of MMA. SYNERGOS accounts for the concomitant activation of all measured muscles rather than only pairs of muscles. This measure results in a single scalar value indicating the overall activity among the set of multiple muscles representing the overall activation of these muscles during the course of the movement. Fig. 1 shows the schematic of the algorithm.

Figure 1.

The schematic view of the SYNERGOS algorithm. EMG signals are analyzed using the RQA and the output, % DET for each muscle is imported into the SYNERGOS algorithm which eventually provides a single scalar index representing the state of MMA.

In this method, the average of several components is calculated. Each component is an average of the mth roots of the products (m= 2, 3, …, n while n= number of recorded muscles) of % DET of EMG signals for sets of m different muscles (i.e., pairs, triplets, quartets, etc.), where the final component is the nth root of the products of % DET of all n muscles being analyzed (i.e. the geometric mean). The SYNERGOS index was calculated based on the equation (1) by using the % DET derived for each muscle.

$SYNERGOS=\frac{1}{n-1}\left[SY{N}_{bi}+SY{N}_{tri}+SY{N}_{quad}+\dots +SY{N}_{n-muscle}\right]$E1

$SY{N}_{bi}$ represents the paired-muscle coactivities with $\left(\begin{array}{c}n\\ 2\end{array}\right)$ possible elements (see equation 2). While $SY{N}_{tri}$ depicts the contribution of simultaneous coactivity among three-muscle sets to SYNERGOS index with $\left(\begin{array}{c}n\\ 3\end{array}\right)$ possible combinations. The same strategy is used to calculate the higher order muscular coactivity among other multiple-muscle sets. SYNERGOS index is elaborated in equation (2) as:

$SYNERGOS=\frac{1}{n-1}\left[\begin{array}{l}\frac{1}{\left(\begin{array}{c}n\\ 2\end{array}\right)}{\displaystyle \sum }_{i=1}^n{\displaystyle \sum }_{j=i}^n\sqrt{\left(1-{\delta }_{ij}\right)DE{T}_{i}DE{T}_j}+\\ \frac{1}{\left(\begin{array}{c}n\\ 3\end{array}\right)}{\displaystyle \sum }_{i=1}^n{\displaystyle \sum }_{j=i}^n{\displaystyle \sum }_{k=j}^n\sqrt[3]{\left(1-{\aleph }_{ijk}\right)DE{T}_{i}DE{T}_{j}DE{T}_k}+\\ \frac{1}{\left(\begin{array}{c}n\\ 4\end{array}\right)}{\displaystyle \sum }_{i=1}^n{\displaystyle \sum }_{j=i}^n{\displaystyle \sum }_{k=j}^n{\displaystyle \sum }_{l=k}^n\sqrt[4]{\left(1-{\aleph }_{ijkl}\right)DE{T}_{i}DE{T}_{j}DE{T}_k.DE{T}_l}+\dots \end{array}\right]$E2

Where $DE{T}_i$ represents the % DET for each muscle (i, j, k, etc.) from the EMG signal during each gait cycle and delta ${\delta }_{ij}$ is the Kronocker delta defined in equation 3:

${\delta }_{ij}=\{\begin{array}{cc}1& i=j\\ 0& i\ne j\end{array}$E3

While for higher order of SYN (i.e. $SY{N}_{tri}$, $SY{N}_{quad}$, …) the $\aleph$ is defined as follows:

${\aleph }_{ijkl\dots }={\delta }_{ij}{\delta }_{jk}{\delta }_{kl}\dots$E4

Therefore for any combination including a pair or multiple similar muscles the ${\aleph }_{ijkl\dots }$ is zero. δ and $\aleph$ are used to ignore the coactivation of each muscle with itself so the result of Equation (1) represents only the coactivity of sets of two or more separate muscles. Following the strategy of the equation (2), $SY{N}_{n-muscle}={\left({\displaystyle \prod }_{i=1}^{n}DE{T}_i\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$n$}\right.}$ and represents the geometric mean of the % DET of the recorded muscles.

This calculation represents the interaction of each muscle’s activity with the other muscles’ activity throughout a dynamic task since these muscles are all active at various times during the task. Each term in Equation 2 calculates the possible coactivation of each muscle with the other muscles in the limb. SYNERGOS basically summarizes the multidimensional correlation tensors of muscle activity into a scalar value while explicitly removing the unidimensional activity of single muscles from such tensors. Consequently, the SYNERGOS method calculates the interaction of all of the muscles, in all possible combinations. The magnitude of SYNERGOS can vary between 0 to 100 indicating the lowest and the highest level of MMA. A SYNERGOS of zero indicates that none of the muscles being measured during a task are simultaneously active regardless of the magnitude of activity in each muscle. A SYNERGOS of 100 indicates that all muscles are simoultanously active and each muscle is activated to its potentially maximum contraction level during each movement cycle (e.g., gait cycle). While theoretically possible, the index could only reach 100 if electrical stimulation was used to achieve tetanus in all monitored muscles. Values between 0 and 100 represent the average simultaneous activation scaled by the magnitudes of activity and temporal sequencing of the respective muscles. During movements a certain coactivity level among several agonist and antagonist muscles exists therefore, a SYNERGOS value between 0 and 100 will be obtained by analyzing the measured muscles.

For instance, in the current study the SYNERGOS algorithm for measuring the level of multiple coactivation for pairs of muscles (i.e., the first component in Equation 2) uses the upper triangular matrix of % DET indicated in equation (5) (where D represents % DET of each EMG signal).

$\left[\begin{array}{cccccc}{D}_{SO}{D}_{SO}& D_{SO}{D}_{GA}& D_{SO}{D}_{TA}& D_{SO}{D}_{VA}& D_{SO}{D}_{RF}& D_{SO}{D}_{BF}\\ & D_{GA}{D}_{GA}& D_{GA}{D}_{TA}& D_{GA}{D}_{VA}& D_{GA}{D}_{RF}& D_{GA}{D}_{BF}\\ & & D_{TA}{D}_{TA}& D_{TA}{D}_{VA}& D_{TA}{D}_{RF}& D_{TA}{D}_{BF}\\ & & & D_{VA}{D}_{VA}& D_{VA}{D}_{RF}& D_{VA}{D}_{BF}\\ & & & & D_{RF}{D}_{RF}& D_{RF}{D}_{BF}\\ & & & & & D_{BF}{D}_{BF}\end{array}\right]$E5

The delta function negates the elements of the matrix located on the diagonal (equation 6) for the EMG from the following muscles: rectus femoris (RF), tibialis anterior (TA), lateral gastrocnemius (GA), soleus (SO), vastus medialis (VM), and biceps femoris (BF).

$\left[\begin{array}{cccccc}0& D_{SO}{D}_{GA}& D_{SO}{D}_{TA}& D_{SO}{D}_{VA}& D_{SO}{D}_{RF}& D_{SO}{D}_{BF}\\ & 0& D_{GA}{D}_{TA}& D_{GA}{D}_{VA}& D_{GA}{D}_{RF}& D_{GA}{D}_{BF}\\ & & 0& D_{TA}{D}_{VA}& D_{TA}{D}_{RF}& D_{TA}{D}_{BF}\\ & & & 0& D_{VA}{D}_{RF}& D_{VA}{D}_{BF}\\ & & & & 0& D_{RF}{D}_{BF}\\ & & & & & 0\end{array}\right]$E6

Finally, the technique calculates the sum of the square roots of the Equation 7 matrix elements. The outcome is averaged over the number of combinations (equation 7).

$SY{N}_{bi}=\frac{1}{\left(\begin{array}{c}6\\ 2\end{array}\right)}{\displaystyle \sum }_{i=1}^6{\displaystyle \sum }_{j=i}^6\sqrt{\left(1-{\delta }_{ij}\right)DE{T}_{i}DE{T}_j}=\dots$E7

Other components of the SYNERGOS requires the calculation of combinations of % DET of three muscles ($SY{N}_{tri}$), four muscles ($SY{N}_{quad}$), five muscles ($SY{N}_5$), and six muscles ($SY{N}_6$) while controlling for the number of combinations ($\left(\begin{array}{c}6\\ 2\end{array}\right)=15$ for two muscles, $\left(\begin{array}{c}6\\ 3\end{array}\right)=20$ for three muscles, $\left(\begin{array}{c}6\\ 4\end{array}\right)=15$ for four muscles, $\left(\begin{array}{c}6\\ 5\end{array}\right)$=6 for five muscles, $\left(\begin{array}{c}6\\ 6\end{array}\right)$=1 for six muscles).

Finally, for each subject, to obtain a single SYNERGOS index for the EMG signals during each gait speed, the root mean square of the five SYNERGOS indices obtained from the clustered EMG signals (five gait cycles per gait speed; see Recurrence Quantification Analysis) were calculated (equation 8). This single value represented the quantified MMA during each gait speed.

$SY{N}_{GS,j}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^{5}SY{N}_{i,j}^2}{5}} , j=1 to number of gait speed conditions$E8

Where $SY{N}_{GS,j}$ is the SYNERGOS index calculated for the jth gait speed condition and $SY{N}_{i,j}$ is the SYNERGOS index calculated for the data clustered in the ith gait cycle ($i=1,2,\dots ,5$) within the jth gait speed condition.

2.5. Statistical analysis

To analyze the efficiency of the proposed method, a restricted maximum likelihood linear mixed model was employed to identify changes in MMA measured by SYNERGOS associated with gait pattern (i.e. walk or run) and with changing gait speed. The model included three fixed effects, speed, pattern (walk or run), and speed-by-pattern interaction, and two random effects, subjects and measurement error (i.e., random within-subject variation). This analytical approach is similar to repeated measures analysis of variance in that it accounts for dependency resulting from multiple measures per subject but unlike analysis of variance does not require the same number of measures for each subject. The fixed effects were used to test the study hypotheses. The random effects were used to compute intraclass correlation coefficients (ICC) of type (2,1) (i.e., degree of consistency among measures) [43] and the corresponding standard errors of measurement (SEm) as relative and absolute reliability estimates, respectively (i.e. indicators of the consistency and precision of the SYNERGOS measure). The significance level was set at p≤ 0.05. The analysis was conducted by using SPSS 16.0.1 (SPSS Inc., Chicago, Illinois, USA).

3. Results

3.1. Surrogate testing

The results of the discriminant statistics (see Appendix-Shuffled Surrogate Tests of EMG). for each muscle and algorithm are shown in Table. 1. For all subjects, muscles, and gait cycle, the results of three different surrogate tests rejected the null hypotheses of equal or more determinism in surrogate data compared to the original data (φ > 2 and p< 5% ). The rejection of this null hypothesis indicated significant change in the nonlinear behavior of the surrogate signals (i.e. reduction of determinism) compared to the collected EMG signals which justified the application of higher order nonlinear data analysis techniques such as RQA to investigate the underlying dynamical pattern of the EMG signals specifically in SYNERGOS.

In Fig.2(a) an example of soleus EMG activity obtained during a single gait cycle (right heel strike to right heel strike) is depicted. Soleus contributes to the ankle planterflexion during body propulsion in stance phase of the gait. Muscle activity dramatically increases during midstance and peaks during the terminal stance phase with a rapid decrease in muscle activity in the pre-swing phase. The muscle remains fairly quiet during the swing phase of gait. Fig.2(b) displays the recurrence plot (RP) of the soleus activity in which recurrent points are positioned along several parallel diagonal recurrent lines demonstrating the existence of a specific deterministic muscular activity in the soleus during the gait cycle. Fig.2(c) represents the randomized shuffling of the EMG signal using surrogate testing algorithms. The RP of the shuffled data are presented in Fig.2(d). No particular deterministic pattern can be recognized by observing several recurrent points scattered on the plot as these points do not generate long diagonal recurrent lines that would indicate determinism of the signal. The analyses of the original and surrogate signals also confirmed the results displayed in Fig.2(b) and Fig.2 (d). The reduction in % DET and increase in ApEn indicates that the surrogate data follows different dynamical patterns than the original EMG signal therefore the EMG signal included nonlinear behaviors (i.e. determinism) that were altered by the randomization of the original signal.

			SO	GA	TA	VA	RF	BF
Time Shuffled	${\phi }_{ApEn}$	Mean	37.54	37.02	24.76	38.99	30.46	27.35
SD	6.52	7.12	7.30	7.81	7.00	3.42
${\phi }_{\%DET}$	Mean	454.65	518.50	344.20	292.60	42.17	213.00
SD	165.09	93.84	111.21	89.51	26.36	68.31
FT	${\phi }_{ApEn}$	Mean	42.16	40.60	32.22	43.79	24.18	25.61
SD	7.43	7.38	12.72	8.92	5.82	4.83
${\phi }_{\%DET}$	Mean	61.26	94.09	232.65	22.43	69.42	48.27
SD	23.08	57.94	46.55	8.34	26.57	9.20
IAAFT	${\phi }_{ApEn}$	Mean	16.63	14.77	18.56	12.47	17.22	15.78
SD	4.16	4.70	6.01	4.59	4.44	1.45
${\phi }_{\%DET}$	Mean	49.24	66.33	53.88	25.84	35.58	11.03
SD	4.95	22.98	41.35	23.47	11.35	4.57

Table 1.

The values of surrogate testing for three different algorithms; ${\phi }_{ApEn}$ and ${\phi }_{\%DET}$ represents the value of statistics calculated from ApEn and % DET of the EMG signals and surrogate data series. (Fourier transform (FT), Iterated amplitude adjusted Fourier transform (IAAFT), rectus femoris (RF), tibialis anterior (TA), lateral gastrocnemius (GA), soleus (SO), vastus medialis (VM), and biceps femoris (BF))

Figure 2.

(a) The neuromuscular activities of soleus muscle during a gait cycle are depicted by EMG activities in in which soleus is mostly active during the propulsion of body during stance phase while showing little activity during swing phase. (b) The Recurrence plots (RP) generated based on the original data is shown in the graph (a) indicating the existence of a specific pattern in the soleus activity during the gait cycle. This pattern is depicted by several recurrent points located along particular diagonal lines which are parallel to the main diagonal line (see Appendix-RQA formulation). The outcome of the RQA also verified the existence of the aforementioned pattern (% REC = 1.98; % DET = 51.80, radius=2.86, ApEn=0.72). Figure (c) contains the randomized shuffled data of the signal shown in 1(a). Figure (d) is the RP of the randomized signal which shows no significant determinism in the shuffled data. The time delayed dimensional data in RP are randomly scattered around the main diagonal line and the recurrent points are positioned along very short length. In addition, the outcome of RQA has shown significant reduction in determinism in the randomized data (% REC = 1.99 radius=8.70, ApEn=1.46). The drastic drop in the determinism of the signal detected by decreasing % DET and increasing ApEn verified the nonlinear dynamics of the EMG signal.

3.2. SYNERGOS analysis

The SYNERGOS indices during walking, running, and the whole protocol are summarized in table2.

Condition	Average	SD
Walking	17.61	11.57
Running	33.86	10.89
Overall	25.71	13.76

Table 2.

The SYNERGOS indices averaged across all subjects and during walking, running, and the full protocol. The average index increased significantly during running compared to walking.

As expected, the changes in the gait mode (walking vs. running) altered the activation of muscles contributing to gait. Fig.3 shows the normalized EMG activities in one of the subjects whose SYNERGOS indices increased significantly from 2.58% during the slowest walking speed to 41.56% during the fastest running speed indicating an increase in cooperation among multiple muscles.

This change was detected by a significant alteration in the SYNERGOS index with the index being significantly greater during running than during walking (F(1,382.149) = 54.067, p <0.001) suggesting increased muscle multiple coactivity during more vigorous movements. In addition, the hypothesis that SYNEGOS could detect changes in MMA associated with slight increases of 0.045 m/s in gait speed was confirmed (F(1,382.537)= 675.85, p< 0.001). A significant interaction effect of gait speed and movement pattern ( F(1,382.082) = 48.075, p <0.001) was also detected indicating that increases in SYNERGOS values occurred at a lower rate during increases in running speed than during increases in walking speed. A high degree of reliability was detected during both walking (ICC=0.92, SEm=3.14) and running (ICC=0.91, SEm=2.71).

Figure 3.

The EMG activity of the lower extremity muscles in soleus (SO), gastrocnemius (GA), tibialis anterior (TA), vastus medialis (VM), rectus femoris (RF), and biceps femoris (BF) comparing walking and running stage. The data are time normalized using a linear length normalization method to convert different each gait cycle into equally scaled units (each unit represents 1% of gait cycle). A gait cycle was defined as the duration between right heel strike and the next right heel strike (gait cycle = 100% of scaled unit). For demonstration purposes, the amplitude for each muscle was normalized to the maximum EMG obtained during the data collection. This maximum was always reached during running.

Fig.4 displays the 95% confidence intervals (95% CI) of SYNERGOS index across all subjects for each of the gait speeds. In addition, the final fitted model slopes both during walk and run in the SYNERGOS indices across treadmill speeds (solid lines) is depicted. The model consistently increased in response to increasing gait speed and had larger absolute values during running, indicating greater magnitude and coactivation (SYNERGOS walk 16.517 ± 3.14 vs. SYNERGOS run = 32.84 ± 2.71). Furthermore, the fitted slopes for increases in the SYNERGOS index during walking (slope = 18.30) and running (slope=15.35) were significantly (p<0.001) different.

Figure 4.

The shaded area indicates the 95% CI of the SYNERGOS index for each gait speed. The vertical dotted line indicates the average walk to run transition speed. The linear fit model for walk and run are shown by solid line while the extrapolation of the linear fit model into the opposite gait pattern was shown by dotted lines for the SYNERGOS indices. In both walking and running conditions the SYNERGOS indices increased consistently in response to increasing gait speed. The index had larger absolute values (>20) during running, indicating greater magnitude and activation of muscle activity, but the slope for walking (slope = 18.30) was significantly (p<0.001) higher than the slope in running condition (slope=15.35) indicating a greater rate of change in MMA when increasing walking speed relative to increasing running speed.

4. Discussion

We have introduced a new analysis method, SYNERGOS that provides an index for quantifying the state of muscle multiple coactivation during a given movement task and demonstrated that it could successfully discriminate between muscle coactivation patterns associated with changing gait mode and speed as a particular combination of the % DET for multiple muscles. As mentioned previously a SYNERGOS index of 100 would represent 100% contraction and simultaneous activation of all measured muscles. Such a possibility is extremely unlikely when considering any voluntary contraction but definitely would not occur during a dynamic movement as rigidity would result.

4.1. SYNERGOS considerations

Several studies investigated the effect of increasing gait speed on the neuromuscular activities of the lower extremities [22, 44-49]. These studies reported increasing average and peak EMG voltage associated with increasing speed in the soleus, gastrocnemius, tibialis anterior, vastus medialis, rectus femoris, and bicep femoris muscles [22, 44-49]. In addition, significant increases in musle cocontraction were reported during increasing gait speed [48]. The increases observed in the SYNERGOS indices are consistent with the previous studies verifying greater EMG activity during faster gait speeds.

Previous studies have indicated that the stability of the human body decreases during higher gait speed which might be correlated with higher risk of fall and injury [50-52]. Muscular cocontraction is a strategy used to stiffen the joints resulting in the reduction of kinematic degrees of freedom to enhance stability during dynamic movements that may threaten postural stability [48]. During faster movements kinematics (velocity and accelerations) and kinetics (i.e. forces, torque, and momentum) parameters alter with higher rates therefore more reliable postural and movement strategy is required to ensure relatively quicker response to the variations in the stability of the system. Thus, increasing the level of MMA provides an effective way for the CNS to reduce the numerous DOF during more demanding movements such as running to provide more stability of human body in faster gait speeds. Increasing SYNERGOS indices are compatible with the aforementioned observation of the motor control strategy.

Several methods for quantifying the coactivity of a group of muscles have been described previously, including muscle cocontraction and various linear techniques. Muscle cocontraction studies are limited to evaluation of only two antagonistic muscles at a time [53, 54]. While analysis of muscle cocontraction may be valuable for various clinical assessments, evaluation of only two muscles does not adequately represent the complex control required to evaluate full-body motion; SYNERGOS overcomes this limitation by offering the potential to represent the combination of EMG activity from all monitored muscles. Linear data analysis methods such as principal components analysis or other factorization techniques have been used to identify the time-invariant patterns of multiple muscle coactivation [2, 4-6, 19, 21, 55], but the nonlinear patterns of information embedded in EMG signals have received little attention [4, 6, 18, 56]. In contrast, SYNERGOS quantifies the changes in MMA of a potentially unlimited number of muscles within the constraints imposed by the practical considerations of the number of muscles EMG can reasonably be collected from during a given movement and analyzes the signals using a powerful nonlinear technique. The SYNERGOS method detects the changes in muscle coactivity states by accounting for both time dependent and time invariant characteristics of EMG signals assessed by % DET EMG without assuming linearity (i.e. stationarity) of EMG signals [13, 14, 16-18].

The SYNERGOS index is an overall estimation of MMA during a specific cycle. The simplicity of the single quantity will come with a price of losing some temporal aspects of muscular activation. Although SYNERGOS algorithm in the first step captures subtle changes in the temporal and magnitude characteristics of each EMG signal by using the % DET in the second step it calculates the overall MMA by averaging the muscular activities using equation (16). Therefore the single quantity cannot demonstrate the exact simultaneous multiple activities of each muscle with others in every single EMG data point. The time-unit of each SYNERGOS index can be set to the duration of the epochs. However this limitation (single SYNERGOS index per epoch) provides simplicity to monitor and track the multiple muscle activations over a longer period of time.

As the CNS likely uses optimized MMA strategies to control different movement tasks, a quantity indicating an overall multiple muscle coactivity may be valuable, particularly in clinical settings, to assess the changes in the performance of the CNS of different patient populations. After further evaluation and validation using data collected during a variety of activities from a variety of patient populations, SYNERGOS may enable clinicians to screen the effectiveness of treatments on neuromuscular activities and could potentially be used as a diagnostic tool to detect abnormal activation in the neuromuscular system.

4.2. RQA considerations

Several previous studies have investigated the application of RQA and % DET to provide insight into the state of muscle activation in various activities during both isometric and dynamic movements. These studies have demonstrated the benefits of such nonlinear techniques to study the neuromuscular activities quantified by EMG signals [13, 14, 16-18, 32]. To perform various activities, muscles are required to generate different forces to satisfy the task related goals that may result in variation of motor unit recruitment and ultimately changes in motor units synchronization [9, 10]. The great sensitivity of RQA to the subtle changes in dynamical systems has increased the use of this analysis for understanding various procedures in motor control, specifically in analyzing EMG signals [13, 16, 18, 29, 35, 40]. However, RQA should be conducted with careful selection of initial parameter settings. % DET has been shown to have high sensitivity to the interaction of noise and embedding dimension if the time lag is more than 8 samples (τ > 8 samples) [29]. Therefore, the selection of the embedding dimension and time lag was conducted using False Nearest Neighbor and Mutual Information techniques for all EMG signals during each gait cycle as the artificial changes in % DET caused by noise may alter the outcome of SYNERGOS.

4.3. EMG considerations

EMG signals depict the overall presentation of action potentials from motor units. Low frequency noises such as power line noise and cable movement are removed using the current technology in EMG data collection devices [57]. Two major sources of noise that may affect the integrity of the EMG signal are baseline noises and movement artifact noise. The baseline noise is generated in electrode amplification process during data collection. In addition, skin movement artifacts are generated during dynamical movement of the muscles resulting in the relative change in the location of the electrode to the targeted muscles. These artifacts can also be generated during highly demanding movement activities in which the impulse of the forces may travel through the muscles and approach the electrodes. In this study a 10-500Hz bandwidth was used to filter the collected EMG signals during the movements. Although in more vigorous activities the corner frequency of 20Hz was shown to remove some additional noises [57], in this study the corner frequency of bandwidth was chosen based on the recommendations of International Society of Electrophysiology and Kinesiology (corner bandwidth frequency of 10 Hz) [30]. As further filtering of data might remove some portions of the ‘true’ EMG signal generated by neuromuscular activation, which might result in limited exposure of actual muscular activities to the RQA, hence dismissing ‘true’ subtle changes in the EMG signal [16, 32] Thus no further filtering was applied on the EMG signals used in this study. Additionally the effect of noise on the % DET as the inputs of SYNERGOS was also minimized by the careful selection of initial parameters (i.e. embedding dimension, delay, and radius) to ensure the integrity of the algorithm (see “RQA consideration”).

In conclusion, we have proposed a nonlinear multiple-muscle coactivation quantification tool, “SYNERGOS”, that is sensitive to changes in both the magnitude and the timing of muscle activity caused by environmental or task related changes. In the future, this method may have application as a diagnostic tool for the evaluation of the therapeutic interventions in individuals with neuromuscular disorders, or those in rehabilitation settings. Further development, validation, and application of the SYNERGOS measure in clinical populations are currently being explored. Additionally, assessment of SYNERGOS’s intra- and inter-day reliability is also underway.

Acknowledgments

We would like to extend our sincere thanks to all the University of Houston students who participated in the study. Additionally, we are thankful to Mr. Chris Arellano, Mr. Marius Dettmer, and Ms. Azadeh Khorram for their help with the data collection.

Apendix

RQA formulation

RQA is a nonlinear data analysis tool, that quantifies the recurrent data in the phase space trajectory of a time series data set, here EMG signals [15, 18]. For RQA analysis each time series is defined as the EMG activity recorded for each muscle during each gait cycle depicted by “$\overrightarrow{d}$”:

$$\overrightarrow{d}={\left[d_1,d_2,\dots ,d_N\right]}^T$$E9

Where $\text{N}$ represents the total number of collected data points in each bin (i.e. EMG data). In our study, N decreased as the gait speed increased, owing to the shorter gait cycles. In the next step the phase space vector is calculated based on the embedding dimension $(\text{m}=6)$ found by False Nearest Neighbor method and time delay ($\text{τ}=0.005\text{ second}$ equal to 5 EMG data samples) found by Mutual Information technique [33], [34]. The equation of the phase space can be constructed as [15]:

${\widehat{\overrightarrow{s}}}_i={\left[d_i,d_{i+\tau },d_{i+2\tau },\dots ,d_{i+\left(m-1\right)\tau }\right]}^T$E10

resulting in a vector with ${\text{N}}_{\text{s}}=\text{N}-\left(\text{m}-1\right)\text{τ}$ elements. Based on the above phase space vector a Distance Matrix (DM) is defined. The elements of the DM are the Euclidian norm of the distance of each of two generated elements of the phase space vector [13].

$D{M}_{i,j}^m=\Vert {\overrightarrow{s}}_i-{\overrightarrow{s}}_j\Vert {\overrightarrow{s}}_i\in {ℝ}^m i,j=1,\dots ,N_s$E11

Here $ℝ$ indicates the real numbers. In the next step, RQA assesses the proximity of each element in the DM with other elements. This proximity is tested based on a predefined threshold radius (${\text{ε}}_{\text{i}}$). In this study, the threshold radius is found by an algorithm to keep the recurrence rate less than 2 percent [18] resulting in $2<{\epsilon }_{\text{i}}<10$ units of the normalized DM by the maximum element in the original DM. Next, the outcome of this assessment is converted into a binary matrix as representing the approximately close elements while 0 indicates the “not-close” elements:

$R_{i,j}^{m,{\epsilon }_i}=H({\epsilon }_i-D{M}_{i,j}^m)$E12

In which H represents the Heaviside function defined as

$H\left(x\right)=\{\begin{array}{cc}0& x<0\\ 1& x\ge 0\end{array}$E13

Equation 12 can be summarized into equation 14 as

$R_{i,j}=\{\begin{array}{cc}0& d{m}_i\ncong d{m}_j\\ 1& d{m}_i\cong d{m}_j\end{array} i,j=1,\dots ,N_s$E14

in which $d{m}_i$ and $d{m}_j$ are two elements on the DM matrix [15]. All elements compared with itself (i.e. i = j) results in the recurrence matrix element of 1. Recurrence plots which visualize the recurrence matrix can be generated based on the frequency distribution of the recurrent points (non-zero elements in recurrence matrix). To calculate the % DET the noncumulative frequency distribution of the constructed diagonal lines (recurrent points) in the Recurrence matrix is defined as

$P^{\epsilon }\left(l\right)= \left\{l_i;i=1,2,\dots , N_l\right\}$E15

Where $N_l$ represents the number of diagonal lines with the length of $l_i$ [17]. Due to the increase in the deterministic pattern of EMG signal resulted from increasing motor unit firing rate during more intense activities (i.e. running) $N_l$ increased for longer diagonal lines and decreased for shorter lines. Finally, % DET was calculated based on equation (16).

$\%DET=\frac{{{\displaystyle \sum }}_{l=l_{min}}^{N}l{P}^{\epsilon }\left(l\right)}{{{\displaystyle \sum }}_{i,j}^N{R}_{i,j}^{m,{\epsilon }_i}}$E16

in which $l_{min}$ is the minimum number of recurrent point in a diagonal line required to define a line [15]. $l_{min}=1$ represents the condition generated by the tangential motion of phase space trajectory [15] that is not indicating the systematic determinism of the recorded EMG signal. In this study, $l_{min}=3$ was chosen to demonstrate the deterministic pattern of the space trajectory based on the recommendation of previous studies [16, 18]. As the denominator of the equation (16) indicates the total recurrent points, % DET measures the proportion of recurrent points that define recurrent diagonal lines longer than lmin representing the determinism or predictability of the dynamical system, i.e. EMG signal.

Shuffled Surrogate Tests of EMG

To obtain a one-tailed significance level of α = 0.05, nineteen (M=19) surrogate data series out of 20 should reject the null hypothesis ($M=\frac{K}{\alpha }-1$ in which K=1). The results of the surrogate testing were evaluated using a one-tail significance level because we had a directional hypothesis that predicted a reduction in the determinism of the surrogate signals after shuffling the original EMG [35, 36, 40, 42]. In the first series, temporally independent surrogate data were generated by random shuffling of the time ordering of EMG data which destroyed any time synchronization and correlation in the original data while saving statistical properties such as the mean and standard deviation [35, 38]. In this step, rejecting the null hypothesis that the EMG signals are originated from white noise is evidence of the existence of a dynamical system in the EMG signal [37]. Next, a phase randomized surrogate algorithm [37] was used to shuffle the original data by these three steps: 1) determining the Fourier transformation of the EMG signals 2) randomization of the phase of Fourier transform 3) applying the inverse Fourier transform to obtain a surrogate time series. The goal of this test is to reject the null hypothesis that the EMG signal has a linearly correlated Gaussian noise pattern. Iterated amplitude adjusted Fourier transform (IAAFT) was the third algorithm to generate surrogate data series [36, 37, 39]. IAAFT algorithm generates surrogate data, which resemble the rank ordering and power spectrum of the original EMG signals. Rejecting the null hypothesis by using IAAFT algorithm can be an indicator of deterministic chaos in the original time series [36, 37, 39, 40], therefore the latter algorithm has the advantage of demonstrating the nature of nonlinearity in the signal.

To test the null hypotheses, discriminating statistics were applied to investigate the changes in the dynamical pattern of surrogate data. These statistics should be sensitive to higher order nonlinearity of the signal [36, 37]. In this study, the null hypotheses were generated based on the changes in the approximate entropy (ApEn) and % DET of the surrogate data. ApEn is a single quantity by which the regularity of a signal can be measured [58, 59]. ApEn has been shown to classify underlying complexity of the signals while no significant changes in frequency and amplitude parameters were detected [40-42]. ApEn rages from 0 to 2 for completely predictable signals (i.e. sine wave) to white Gaussian noise respectively. A value of ApEn=2 indicates complete uncertainty in prediction of future behavior of a dynamical system. Therefore, in surrogate testing the expected outcome is a significant increase in ApEn value while a significant drop in % DET as a result of random shuffling procedure [35, 38].

ApEn value for both original EMG signals and surrogate data were calculated by applying the parameter settings of embedding-dimension of m=2 and r= 0.2 × standard deviation of the signal. In addition, after conducting the RQA analysis on the original EMG signals, a Matlab script was used to perform the RQA on the surrogate data using the same parameter settings ( m=6 and τ= 0.005 second). The script increased the radius to obtain similar level of % REC (the percentage of recurrent points in the recurrence plot graphs) for each set of muscles per gait cycle and the % DET was calculated base on the modified radius.

Finally, the ApEn and % DET of surrogate data were statistically compared to the ApEn and % DET of original EMG signal by defining the following statistics [37, 40]:

${\phi }_i=\left|Q_{original}-{\overline{\mu }}_{{\left(surrogate\right)}_i}\right|/S{D}_{{\left(surrogate\right)}_i} i=1 to n$E17

Where n is the number of muscles and $Q_{original}$ indicates ApEn or % DET statistic from original EMG for each muscle per gait cycle. $${\overline{\mu }}_{{\left(surrogate\right)}_i}$$ and $S{D}_{surrogate}$ denote the average and standard deviation of computed ApEn and % DET of the surrogate series. ${\phi }_i$ indicates the amount of change in the ApEn and % DET of the original data in the scale of standard deviations. To reject the null hypothesis a minimum ${\phi }_i>2$ is required to obtain a 5% significance level (normality assumed) [39].

Nonlinear data analysis techniques are capable of revealing subtle changes in dynamical systems that may be ignored during linear data analysis. However they require more sophisticated analysis procedures and are generally more time consuming and costly, therefore the application of such techniques should be justified especially during clinical measurements. In the current study, three algorithms were used to test the state of nonlinearity in the EMG signals. The first two algorithms (time shuffled and FT) confirmed the fact that the recorded EMG signals contain higher order nonlinear dynamics. The use of the third surrogate testing method (i.e. IAAFT algorithm) expands our understanding from the nature of the dynamical nonlinearity. In our investigation, the existence of deterministic patterns measured by % DET is a key to the SYNERGOS equation. Therefore the use of third surrogate algorithm, IAAFT was justified. In previous studies amplitude adjusted Fourier transform (AAFT) algorithms were used to investigate the nature of the EMG signal [40] however it has been argued that for short and highly correlated data series the AAFT algorithm may result in flatness of power spectrums [36] therefore IAAFT algorithm was introduced to overcome such a bias by iteratively correcting the deviations in the power spectrum. [36, 39].