Effect of Robotic Exoskeleton Motion Constraints on Upper Limb Muscle Synergies: A Case Study

Evidence exists that changes in composition, timing, and number of muscle synergies can be correlated to functional changes resulting from neurological injury. These changes can also serve as an indicator of level of motor impairment. As such, synergy analysis can be used as an assessment tool for robotic rehabilitation. However, it is unclear whether using a rehabilitation robot to isolate limb movements during training affects the subject’s muscle synergies, which would affect synergy-based assessments. In this case study, electromyographic (EMG) data were collected to analyze muscle synergies generated during single degree-of-freedom (DoF) elbow and wrist movements performed by a single healthy subject in a four DoF robotic exoskeleton. For each trial, the subject was instructed to move a single DoF from a neutral position out to a target and back while the remaining DoFs were held in a neutral position by either the robot (constrained) or the subject (unconstrained). Four factorization methods were used to calculate muscle synergies for both types of trials: concatenation, averaging, single trials, and bootstrapping. The number of synergies was chosen to achieve 90% global variability accounted for. Our preliminary results indicate that muscle synergy composition and timing were highly similar for constrained and unconstrained trials, though some differences between the four factorization methods existed. These differences could be explained by higher trial-to-trial EMG variability for the unconstrained trials. These results suggest that using a robotic exoskeleton to constrain limb movements during robotic training may not alter a subject’s muscle synergies, at least for healthy subjects.

time-invariant synergy vector and an associated time-varying neural command. Muscle synergies have been characterized in the human wrist for the purpose of neurorehabilitation [4]; however, the wrist is not as well studied as the shoulder and elbow, for which there is evidence that changes in muscle synergies can be correlated to changes resulting from neurological injury [5] as well as the level of impairment [6].
Robotic exoskeletons are a tool for rehabilitation interventions, and are often used for assessment of recovery throughout a rehabilitation protocol [7]. We want to add synergy analysis as another assessment tool in our robotic rehabilitation protocol. Additionally, we want this assessment to take place while the patient is in the exoskeleton, because getting in and out of the device can be time consuming and challenging for individuals with motor impairments. This raises the question of how the presence of the device acting with the limb might alter synergies. The effect of wearable robotic devices on muscle synergies has been studied in multiple circumstances, including point-to-point arm reaching movements moving freely versus moving inside an exoskeleton [8] and for activities of daily living performed inside versus outside an exoskeleton [9]. All results indicate that the shape of the muscle synergy vectors was preserved across conditions; only the activations of the synergies were changed.
One objective of this paper is therefore to analyze muscle synergies for the elbow and wrist, while a subject is inside the MAHI Exo-II (MEII) exoskeleton, and determine how constraining the limb to certain joint-isolating movements might alter synergies. During robotic training, patients using the MEII perform repetitive, single degree-of-freedom (DoF) movements with variable assistance/resistance. While doing so, the remaining three DoFs may be held stationary by the robot. In order to analyze the effect of those robot constraints on muscle synergies, we have collected EMG from repeated trials of single-DoF movements with and without active robot constraints present. How each anatomical degree of freedom is handled under the two conditions, referred to as constrained (C) and unconstrained (U), is specified in TABLE I.
The use of repeated trials of the same movement raises a methodological question of how to properly combine repeated trials when conducting muscle synergy analysis. Methodological choices are known to have an effect on muscle synergy analysis, but there is still a lack of standardization and guidance [10]. We would like to study how different approaches to analyzing data from repeated trials affects the outcomes of muscle synergy analysis.
Recordings of EMG under different conditions and for different trials are often pooled together before muscle argue that averaging EMG recordings from repeated trials is necessary to boost the signal-to-noise ratio of the EMG signal, which is stochastic in nature [11]. They simultaneously acknowledge, however, that averaging repeated trials also obscures information about any muscle activity that is correcting for the random physiological variations that are known to occur in human motor control. The question of how to combine repeated trials from walking for muscle synergy analysis was examined more directly by Oliveira et al. [12], who compared three approaches: factorization of the trials individually (SNG), concatenated (CNC), and averaged (AVG). Using synergy vectors identified from the three methods, they found that SVs from concatenation (CNC) were able to reconstruct longer periods of muscle activity with the highest accuracy. In particular, they recommend concatenating at least 20 step cycles to accurately capture the variation of individual muscle activities.
A secondary contribution of this paper is therefore to explore methodological differences in synergy analysis when data from repeated trials are available. We use the concatenated (CNC), averaged (AVG), and individual trial (SNG) methods of factorization, as well as a bootstrapping technique (BTS), to explore the impact the choice of method has on our results. This paper helps to lay the groundwork for future studies involving muscle synergy analysis while the human is moving inside of a robotic exoskeleton. First, we show that, for our single subject case study, muscle synergy composition and timing are not significantly affected when a robot is used to constrain human movements, compared to unconstrained movements made while a participant's limb is in the robot. Second, we show that the methodological choice of how to combine repeated trials of arm movement data does have an impact on synergy analysis results.

A. Case Study Design
We designed a study to compare muscle synergies identified under two control modes: unconstrained (U) in the stationary DoFs, and constrained (C) in the stationary DoFs. A single participant (male, age 23) provided informed consent according to the protocol approved by the Rice University Institutional Review Board (FY2018-29). For each of the two modes, we collected sufficient data to identify muscle synergies underlying the 4 single-DoF movements of the MAHI Exo-II: elbow flexion/extension (E-Flx/Ext), forearm pronation/supination (F-Pro/Sup), wrist flexion/extension (W-Flx/Ext), and wrist radial/ulnar deviation (W-Rad/Uln). For each degree of freedom, the participant moved from a neutral joint position to a target position at either end of the range of motion, and then back to the neutral position. In this way, they were duplicating typical movements that would occur inside the exoskeleton during a training session with a physical therapist. This combination of an outbound (OB) and inbound (IB) movement was repeated 25 times for each of the 2 possible targets, for each degree of freedom. The result was 400 trials in the unconstrained mode to be compared with 400 trials in the constrained mode.

B. Hardware
The MAHI Exo-II, shown in Fig. 1, is a robotic exoskeleton designed for the rehabilitation of the elbow and wrist joints [13]. It features serially connected joints for elbow flexion/extension (E-Flx/Ext) and forearm pronation/supination (F-Pro/Sup), and a parallel revolute-prismatic-spherical mechanism that achieves wrist flexion/extension (W-Flx/Ext) and wrist radial/ulnar deviation (W-Rad/Uln).
In all phases of the experiment, the robot was position-controlled with proportional derivative feedback. When robot motion was required, the robot followed a pre-programmed reference trajectory.
The Delsys Bagnoli EMG system provides eight channels of surface EMG data. The variable gain for the channel amplification was set to 1000, and analog band-pass filtering was applied by the Bagnoli system at 20 Hz-450 Hz to remove any movement artifacts or aliasing. All robot and EMG data were acquired with the Quanser Q8-USB and sampled at 1kHz.
The height and shoulder abduction angle of the MAHI-Exo II were adjusted so the subject could hold their arm in a natural position with the elbow flexed while seated. The position of the chair relative to the exoskeleton was adjusted to keep both shoulders at equal heights and to keep the shoulder in the scapular plane (∼30 • from the frontal plane). Subjects were instructed not to move their torsos or shoulders during  testing but restraints were not used to enforce this. The wrist handle location was positioned to provide a maximum range of motion while the subject held it in a natural grip.
The MAHI Exo-II is equipped with an adjustable counterweight for passive gravity compensation of the elbow joint. The elbow was set so that the subject was able to be at rest while their elbow was flexed 90 • , which was defined to be the neutral elbow position. The exoskeleton can be configured for left and right handed individuals, so the dominant right arm was used by the subject. Once inside the exoskeleton, the participant was strapped to the robot at the upper forearm and the hand.

D. Kinematic Data Processing
During an experimental session, robot encoders measured joint angles while EMG electrodes recorded activity on all 8 selected muscles. A real example of the data generated over a short window containing several trials is given in Fig. 2.
The first step of our data analysis was to segment the data into individual trials based on the joint velocity at the moving joint. A 4th-order low-pass Butterworth filter with a cutoff frequency of 30Hz was first applied to the velocity signal, as can be seen in Fig. 2. Then, a fixed velocity threshold of 0.01rad/s was used to mark the beginning and end of a movement. Rohrer et al. [15] used a velocity threshold of 2% of the peak speed of a given trial, but we chose a fixed value to avoid any sensitivity to variation in velocity peak magnitude from trial to trial. The selected threshold of 0.01rad/s was typically around 0.5%-1% of the peak velocity of a trial. Finally, for each trial, we identified the contiguous segment of movement above the velocity threshold with the largest velocity integral-which was also moving in the correct direction-as the motion segment to be analyzed.
Next, we screened for movements that were not kinematically similar to the rest, since we wanted to characterize repeatable patterns of EMG activity. Shapes of velocity profiles were compared to a typical "bell curve", or "minimum jerk trajectory". Highly skewed segments were trimmed so that they were better aligned with more symmetric segments, and trials that did not pass all shape-and magnitude-based criteria were removed from the data set. The remaining trials were grouped based on all factors (mode, DoF, target, direction) and then checked for outliers in both normalized position and non-normalized velocity. Any outliers were also removed from the data set.
Finally, the position trajectories of the stationary DoFs were checked for significant movement. For each trial, if the user moved more than 50% of the distance from neutral position to target position in a stationary degree of freedom, then that trial was considered to be a multi-DoF movement and was removed from the data set.
All remaining trials were then counted in their groupings based on the experimental factors. The smallest group was composed of 16 trials (unconstrained mode, outbound elbow flexion, i.e., going from elbow neutral to elbow flexed). In order to have an equal number of trials in each group, 16 trials were selected for each group based on which had the highest variance across the EMG signals.

E. EMG Data Processing
All EMG data were acquired at a sampling rate of 1kHz. After data collection, the following signal processing steps were applied to each EMG channel to convert the raw EMG signal to a non-negative envelope: 4th-order Butterworth band-pass filter at 20 Hz-450 Hz (for removal of movement artifacts and aliasing), full-wave rectification (absolute value), 4th-order Butterworth low-pass filter at 8.05Hz, full-wave rectification (absolute value). Fig. 2 shows an example of the raw EMG recordings after band-pass filtering in light gray, with the non-negative EMG envelope overlaid in black. The low-pass filter cutoff frequency of 8.05Hz was derived from the estimated average duration of a movement cycle, using the formula f cut o f f = 7/dur. This formula has been used for synergy analysis of muscles in the lower limb involved in walking [10]. A movement cycle was defined here as a pair of outbound and inbound movements and was calculated to be two times the average duration of all movement segments.
The start and end of each trial, as determined from the kinematic data, were used to segment the EMG envelope after all filtering was applied. The start of each EMG segment was shifted backward by 200ms relative to the onset of movement, due to the physical delay between an EMG signal and the actual joint motion that it generates [1]. The end of each EMG segment was not shifted in order to capture all relevant EMG activity.
Normalization of EMG data is generally required for muscle synergy analysis. The maximum value recorded on each of the 8 channels was found within the 16-trial data set, and each channel was then divided by this maximum value so that it would reach a maximum value of 1 within the data set being analyzed. This is equivalent to the "MaxOver" method of normalization used in [10].
Lastly, the data were resampled in time using linear interpolation to be 101 data points long at equal spacing. Hence, every trial now contained samples at normalized discrete times of 0 to 100 [10].

F. Synergy Factorization
Given a time series of EMG sample vectors x(t) concatenated into a N-by-M matrix X, where M is the number of time samples, the approximation of the data by K synergies can be rewritten as X ≈ W H . The number of synergies, K , is the inner dimension of the matrix product and is less than The algorithm most commonly used for finding the elements of W and H that best approximate X is called non-negative matrix factorization (NNMF), and the details of its implementation are given in [16]. After factorization, the product of W and H can be considered a reconstruction of the original data, X, and metrics of comparison between the two can be used to determine the quality of reconstruction.
We have chosen to compare four different methods of factorization of this data set. All four approaches make use of the MATLAB function nnmf to perform non-negative matrix factorization, and only differ by the input data given to the function.
The EMG envelopes of individual trials that are created by our data segmenting process are each stored as an 8-by-101 matrix X trial , illustrated in Fig. 3. For all four methods, the trials were first grouped by experimental condition, meaning each group of 16 trials has the same control mode, degree of freedom, target, and direction (e.g., unconstrained, outbound elbow flexion, abbreviated U OB E-Flx). This first step has also been illustrated in Fig. 3. There are 32 such groups in total, 16 for each control mode. In all four methods, we wish to end up with SVs and NCs from the unconstrained and constrained group that can be directly compared, and we want all 16 groups of same-condition trials to contribute equally to the synergy decomposition.
The first method of factorization is concatenation (CNC) of all trials (Fig. 3). Trials from each condition are horizontally concatenated, resulting in a 8-by-101·16 data matrix. Then, for each control mode, all 16 groups are horizontally concatenated, resulting in an 8-by-101·256 data matrix. Factorization of a single X C NC for each mode results in a single set of synergies for the unconstrained mode to be compared to a single set of synergies for the constrained mode.
The second method is averaging (AVG) of trials from the same condition (Fig. 3). The EMG envelopes from each channel are averaged for each group of trials, resulting in an 8-by-101 data matrix. Then, for each control mode, all 16 groups are horizontally concatenated, resulting in a 8-by-101·16 data matrix. Again, factorization of a single X AV G for each mode results in a single set of synergies for the unconstrained mode to be compared to a single set of synergies for the constrained mode.
The third method is repeatedly concatenating single (SNG) trials from each condition, for all 16 trials (Fig. 3). More specifically, for each group we assign each trial a number 1 through 16 according to the order in which they were recorded, and then horizontally concatenate all trials of the same number. This yields 16 data matrices X n S NG of size 8-by-101·16, each containing a single trial from each condition. Each X n S NG is then factored into a set of synergies, meaning there are 16 sets of synergies for the unconstrained mode to be compared with 16 sets of synergies for the constrained mode. Grouping trials based on the order in which they were recorded could be considered a naïve way to sample the repeated trials, given there is no meaningful relationship between the trials that have been concatenated together for factorization. However, it is a simple and fast way to generate a distribution of synergies from data sets with repeated trials, which is sometimes desirable.
The fourth and final method is using bootstrapping (BTS) to resample the trials from each condition over many repetitions, and then perform factorization as in the first method (CNC) with concatenation of all trials (Fig. 3). Bootstrapping is a strategy commonly used in statistics to obtain information about the uncertainty of a parameter, and has been used to resample repeated trials of EMG data for muscle synergy analysis [17]. The purpose of this method is to obtain an estimate of the distribution of synergies seen from a population similar to the data sample from this single subject [18]. When applying the bootstrapping technique, it is important that the data are always resampled at the same size as the original sample, and that it is sampled with replacement. This means that for each condition, one of the trials was randomly The bootstrapping procedure must be repeated many times for the estimate to be accurate, so we chose to use 100 repetitions. Each repetition yields a data matrix X n BT S of size 8-by-101·256 to be factored into a set of synergies, resulting in 100 sets of synergies for the unconstrained mode to be compared with 100 sets of synergies for the constrained mode. Because factorization needs to be performed at every repetition, the procedure is computationally costly, which is why 100 repetitions were chosen instead of 1000 as in previous analyses of muscle synergies [17].
All major results will be reported for each of the four methods described, when appropriate.
1) Factorization Repeatability: The NNMF algorithm is not guaranteed to find a global minimum in the residual between W ·H and X [16]; therefore, the algorithm must be repeated a sufficient number of times with randomized initial conditions in order to calculate a likely global minimum.
Using nnmf function in MATLAB to perform NNMF, the parameter settings listed in TABLE II were necessary to obtain sufficiently repeatable results for this data set. We used the multiplicative update ('mult') version of nnmf instead of alternating least-squares ('als'), as is recommended by MATLAB documentation when using multiple replicates. Values were tuned for X C NC , assuming they would also apply to the X n BT S . Likewise, values were also tuned for X AV G , assuming those values would also apply to the X n S NG . Some of the parameters given in TABLE II could be adjusted for better repeatability, but this choice is ultimately a tradeoff between precision and computation time. It is also noteworthy that MATLAB documentation suggests that for most optimization problems, setting 'TolFun' and 'TolX' any smaller than 1 × 10 −6 will not have much effect. This was not found to be the case for synergy factorization, particularly when working with the smaller data set, X AV G .
2) Synergy Normalization, Matching, and Sorting: There is no one correct way to scale the factorization output, and various methods have been used successfully [10]. We have chosen to normalize the synergy vectors to have unit magnitude, thus leaving the neural commands to be scaled accordingly without other constraints. Normalizing SVs to have unit magnitude focuses on the relative weighting of the specific muscles within a synergy as the most important characteristic, and then interprets NCs as conveying the level of recruitment of an individual synergy.
After normalization, SVs from different data sets can be matched based on their similarity in shape, which we have calculated using the cosine angle-the normalized inner-product of two vectors. Therefore, SVs to be compared between the unconstrained and constrained conditions were provided to the K-means clustering algorithm to identify individual clusters (synergies) based on the cosine distance metric. Checks on the clustering algorithm's output were added to ensure that for each set of synergies coming from an individual factorization, one synergy was assigned to each cluster. These checks were passed in all data analysis, indicating that the synergies were easily differentiated from each other. Finally, synergies were sorted by maximum weighted muscle and magnitude (while preserving matching) for consistency of presentation.
3) Data Dimensionality: The dimensionality of a data set is used in this context to refer to the minimum number of independent variables-synergies-needed to accurately reconstruct the original data. The 8-channel EMG data that we collected can be factored into any number of synergies between 1 and 8, 8 being the smaller dimension of the data matrix X. Choosing 8 synergies simply returns the original 8 EMG channels unaltered. Therefore, it is not considered in our analysis. Choosing a number smaller than 8 will generate synergy vectors that are weighted combinations of the different EMG channels, similar to a principal component analysis.
In order to choose the optimal and parsimonious number of synergies, we examined metrics for the quality of reconstruction at each number of synergies. The dimensionality of the data was then estimated to be the minimum number of synergies that met our criteria for goodness of reconstruction. To avoid over-fitting, the data were split into "training" and "testing" subsets when estimating the number of synergies needed-a technique known as cross validation [19]. Factorization was performed on the training data to obtain a set of synergy vectors, which were then used to reconstruct the testing data. The process of randomly splitting and factoring the data was repeated 100 times, as in [20], with checks to ensure no permutations were repeated.
Cross validation had to be carried out differently for each of our factorization methods. For the CNC and AVG methods, each data set was randomly split into a training set and a testing set, with 50% of the data used for training (8 out of 16 trials) and the remaining 50% used for testing. The trials were split such that for each type of trial, 8 went to the training set and 8 went to the testing set. For the SNG method, only one trial of each type can be used at a time; therefore, a single trial number (the n index in Fig. 3) was randomly assigned to training and another single trial number was randomly assigned to testing, for each iteration. Finally, for the BTS method, trials were randomly split into 50% training and 50% testing, as in the first two methods. Then, the bootstrapping procedure of repeated random sampling with replacement was applied to the training and testing data. Due to the computational cost of this procedure, it was only done once. After each factorization, we compute the quality of reconstruction of the original data using the variance accounted for (VAF). The global VAF has the simple formula 1− SS E/SST , where SS E is the sum of squares due to error, and SST is the total sum of squares [21]. VAF has become a commonly used metric for muscle synergy analysis [22], [23]; this statistic is generally more stringent than the Pearson's correlation coefficient (r-squared), because it measures similarity in both the shape and magnitude of the signals [22].
When the VAF for a certain data set are plotted versus the number of synergies used in the reconstruction, all VAF values begin to approach 1 as the number of synergies nears the number of EMG channels, 8-at which point all values must equal 1. The optimal number of synergies was chosen to be the number at which the VAF crossed above a threshold of 0.9, or 90%, similar to threshold criteria used in previous studies [2], [10], [24].
Finally, the cosine similarity metric comparing two vectors that was used in the synergy vector matching algorithm was also used to compare synergy vectors between the unconstrained and constrained modes. A similarity of 1 indicates that the vectors are identical, and a similarity of 0 means the vectors are orthogonal. Cosine similarity was calculated between synergies from the unconstrained and constrained modes for each method of factorization.

A. Data Dimensionality
After cross-validation, the mean of the global VAF curve for the unconstrained data set (shown in Fig. 4) crosses the chosen threshold of 90% at 3 synergies, regardless of the method It is easiest to directly compare muscle synergies between two data sets when the same number of synergies have been used in factorization. Therefore, in accordance with the VAF curves found from the CNC, SNG, and BTS methods of factorization, 4 synergies were extracted from the data sets for further analysis and visualization, as it is the number at which the VAF for both data sets has passed the threshold.

B. Synergy Factorization
Factorization of the concatenated (CNC) data into 4 synergies (Fig. 5) reveals a high degree of similarity between the unconstrained and constrained modes. The synergy vectors, represented by the bar graph in Fig. 5(a), each have significant contributions from only one or a few of the muscles recorded. For each of the four synergies extracted from the unconstrained data set (see Fig. 5(c)), there is an obvious counterpart from the constrained data set (Fig. 5(d)) that was found through our automated synergy matching procedure.
The first synergy is responsible for the activation of BB (biceps brachii, elbow flexion), the second synergy is mainly responsible for the activation of PT (pronator teres, forearm pronation) and FCR (flexor carpi radialis, wrist flexion), the third synergy is mainly responsible for the activation of S (supinator, forearm supination) and ECU (extensor carpi ulnaris, wrist etension), and the fourth synergy is responsible for the activation of ECRL (extensor carpi radialis longus, wrist radial deviation). Thus, the first synergy is associated purely with elbow flexion, and the remaining three synergies are blended in their contributions to the movements of the wrist, but with clear preference for direction.
The neural commands shown in Figure 5(b) include only one trial from each of the 16 subsets (see Fig. 3 Key), as all 256 trials used in this data set could not be visualized. Light gray vertical dotted lines separate the individual trials. The plots shown are, however, representative of the full data set, and also show similarity between the unconstrained and constrained data sets. There is one notable difference found in the neural commands of the first synergy, which is effectively the activity of the biceps, during trials of forearm pronation/supination. Here, the biceps activity is near zero in the constrained case because the robot is supporting the weight of the arm, while the biceps is actively maintaining the flexed arm position in the unconstrained case.
The plots in Figures 5(c) and 5(d) illustrate how well the 4 synergies can reconstruct the EMG input data, with separate plots showing the unconstrained (darker) and constrained (lighter) data sets. A gap between the EMG envelope and the summed synergy contributions indicates EMG activity that is not described by the muscle synergy factorization. For example, it can also be seen that with just 4 synergies we cannot describe the activity of the triceps, which is present only during elbow extension movements. Otherwise there are only a few small gaps between the recorded and reconstructed signals, as is expected when the global VAF is above 90%. The individual EMG channels are each composed of mainly one synergy.
These factorization results can be compared with those of the averaging (AVG) method in Fig. 6, the single trial (SNG) method in Fig. 7, and finally the bootstrapping (BTS) method in Fig. 8. The muscle synergies found by these other methods were mostly identical, with the third synergy being the exception.
In Fig. 6, the results of the single AVG factorization show that the third synergy heavily weights the supinator (S) in the unconstrained case and the flexor carpi radialis (FCR) in the constrained case. It can also be seen that the activity of the supinator is absorbed into the first synergy for the unconstrained case, unlike in the CNC results.
The SNG and BTS methods generated multiple factorizations for the two conditions being compared, and so the   Fig. 7 and Fig. 8 are shown as means and standard deviations. Focusing on the distributions of the synergy vectors, both the SNG and BTS methods reveal large variations in the primary muscles of the second and third synergiesmuscles contributing to wrist pronation/supination and wrist flexion/extension. The trial-to-trial variation is, however, much more significant for the data taken in the unconstrained condition.
The visual comparison of the synergy vectors in Fig. 5 -8 is quantified in TABLE III with the cosine similarity metric. The first, second, and fourth synergies are highly similar after all four methods of factorization, but the third synergy (mainly involved in wrist supination and extension) was shown to differ, depending on the method of factorization.

IV. DISCUSSION
For this data set collected from a single subject, our estimate of the data dimensionality using VAF showed small differences in the number of synergies to be selected for the unconstrained vs. constrained condition. The global VAF curves shown in Fig. 4 were very similar across three of the methods of factorization, while the averaging method (AVG) yielded a slightly different pair of curves that pointed toward the same number of synergies being needed for the unconstrained and constrained conditions. This result highlights the fact that there is a more fundamental difference between the AVG method and the other three methods of factorization-in AVG the EMG data of multiple trials was averaged before factorization, while in the other methods EMG data were never averaged across trials. Averaging EMG data across trials obscures variations in the shape and timing of muscle activation for individual trials, which is not necessary for the NNMF algorithm to uncover the synergistic activation of different muscles. Therefore, we do not recommend this method of factorization when looking to extract synergy vectors and compare them across conditions.
Regarding the factorization results from the CNC method in Fig. 5, all four synergy vectors extracted from the two data sets are very similar in shape; the cosine similarity metric reported in TABLE III confirms this. Taken by itself, this is preliminary evidence that the active robot constraints did not alter the muscle synergies in this task.
In contrast, the factorization results from the AVG method, shown in Fig. 6 and TABLE III, report less similarity in the two sets of synergies. The third synergy, originally contributing to wrist supination and extension, appears to take on a different role for the unconstrained data set. However, the results of the SNG and BTS methods in Fig. 7 and Fig. 8 offer a clearer explanation for this difference. The large error bars on the weighting of certain muscles in the synergy vectors-mostly for the unconstrained data set-indicate the variability in synergy vector shape based on the trials included in factorization.
The larger trial-to-trial variation in synergy vector shape for data collected in the unconstrained case is expected, as the user is not being aided by the robot to move on a precise trajectory. Ultimately, the synergy vectors generated by the CNC, SNG, and BTS methods show a close similarity between the unconstrained and constrained data sets, TABLE III. The AVG method distorts the information contained in the trialto-trial variation, making it seem that the third synergy is very different between the two conditions, while the SNG and BTS methods clarify that though there is variability in the trials, synergies vectors from the two conditions are much more similar than not.
Analyzing the factorization results across the four methods used, the data from this single subject show that when the robot is actively constraining the user to move a single degree  of freedom of elbow or wrist, the variability in the muscle synergies from trial to trial is reduced, but the underlying muscle synergies of the task are not altered by the robot constraints.
Finally, the CNC, SNG, and BTS methods of factorization provided results that were consistent with each other. The CNC method, however, could provide no information on inter-trial variation in synergies. The BTS method was included in this analysis as a more appropriate (and more time consuming) way to repeatedly sample the distribution of trials and estimate the distribution of the synergy vectors, when compared to the SNG method. Still, the results of the SNG and BTS methods were very similar, which supports the possible usage of simpler methods such as SNG when approximating synergy distributions from repeated trials.
A major limitation of this study is its inclusion of only a single participant, which limits our ability to draw more general conclusions from this data set. Due to the significant computational algorithm development required for this work, we have chosen to focus the present study on providing a detailed explanation of the methods developed to answer the questions being considered. This work is also limited by its ability to only investigate a part of the robotic assistance necessary in a rehabilitation protocol, namely, the constraint of stationary degrees of freedom. In a typical rehabilitation setting, these constraints would be combined with active robotic assistance or resistance at the moving DoF, depending on the patient's ability. However, both forms of robotic interaction affect the physical task being performed by the patient, and it is reasonable to study them in isolation to understand how each component of the human-robot interaction can potentially alter muscle synergies.

V. CONCLUSION
In this paper, we explored the identification of synchronous muscle synergies in the muscles controlling the elbow and wrist, and the possible effects of robot-imposed task constraints on the "neural constraints" represented by muscle synergies. For a single subject, we collected EMG from repeated trials of single-DoF movements with and without active robot constraints present, and compared the resulting synergy vectors calculated from the two conditions. We also examined how different approaches to analyzing data from repeated trials affects the outcomes of muscle synergy analysis.
The results of this single-subject case study suggest that active robot constraints for isolating single-DoF movements of the elbow and wrist do not alter the muscle synergy composition, and that the question warrants further investigation. We extracted 4 muscle synergies from the data collected in the unconstrained and constrained conditions, and found the shapes of these synergies to be very similar (>0.94 cosine similarity) when all of the trials were concatenated (CNC).
We found that the methodological choice of how to combine repeated trials of arm movements did have an impact on the results in the case of averaging trials versus concatenation. In particular, the method of averaging trials before factorization (AVG) obscured the fact that all four synergies extracted from the two data sets were similar in shape, despite there being some variability in the synergy shapes introduced by the repeated trials. The two methods introduced for producing a distribution of synergy factorizations from the repeated trials (SNG and BTS) allowed for more in-depth analysis of the results. Future work involving more participants is needed to verify this conclusion as well.
Most importantly, our results suggest that constraining the stationary degrees of freedom during a single-DoF movement inside the exoskeleton has little to no effect on the underlying muscle synergies used in the task, while simultaneously reducing trial-to-trial variability in the muscle synergies extracted from the data. If these conclusions hold for a greater number of participants, then this knowledge can be applied to designing future experiments that investigate the robot's ability to change a user's muscle synergies for the purpose of neurorehabilitation.