Control of Fast-Reaching Movements by Muscle Synergy Combinations

Experimental design and apparatus.

Standing subjects gripped a handle with their right hand and were instructed to move a sphere (diameter, 4 cm) attached to the handle from a start location to a target location and, in experiment 3, from a start location to an intermediate location (via-point) to a target location. The handle (see Fig. 1a) was constructed with two cylinders: one for gripping (diameter, 4 cm; covered with tennis grip tape) and one with the sphere attached on its side. The two cylinders were coplanar, with their symmetry axes at an angle of 16°, so that the axis of the cylinder with the sphere was perpendicular to the forearm axis when the subject gripped the other cylinder while keeping a natural posture of the wrist with a slight ulnar deviation. When the handle was being gripped, the sphere was aligned along the forearm axis with its center at a distance of 12 cm from the center of the gripping area. In experiments 1 and 3, the subjects were instructed to keep the handle vertical, thus imposing a neutral posture with respect to forearm rotation. In experiment 1, the weight of the handle was varied by inserting a tightly fitting cylinder filled with small lead spheres into the gripping cylinder (see Fig. 1b). Two different loads were used so that the total weight of the handle varied among three levels: 180, 630, and 1040 g. In experiment 2, the subjects were instructed to perform reaching movements with the unloaded handle (weight, 180 g) and to vary the forearm posture by rotating the handle into a horizontal orientation after either a clockwise or a counterclockwise rotation (see Fig. 1c). Thus the forearm posture varied among three conditions: neutral, pronated, and supinated.

Nine right-handed subjects (age range, 27–40 years) participated in the experiments after giving informed consent. Three subjects (1, 2, and 3) performed experiment 1, three subjects (4, 5, and 6) performed experiment 2, and three subjects (7, 8, and 9) performed experiment 3. All experimental protocols conformed to the Declaration of Helsinki on the use of human subjects in research. Subjects performed blocks of movements from one start location to a target location, directly (point-to-point movements) or through a via-point (reversal and via-point movements), either in the sagittal or in the frontal plane (see Fig. 1d). Subjects were instructed to reach using only their arm and to minimize trunk movements. The start location was marked by a sphere whose location was adjusted so that the arm was approximately in a parasagittal plane, with the upper arm vertical along the trunk and the forearm horizontal. For the four subjects (1, 2, 4, and 5; apparatus 1; see Table 1) who performed only point-to-point movements, the target locations were marked by a sphere positioned by an experimenter. For the remaining five subjects (apparatus 2), the start, target, and via-point locations were marked by transparent spheres mounted on a supporting structure and illuminated from inside by a light-emitting diode. In all cases, the peripheral targets were positioned approximately on a circle of 30 cm of radius centered in the start location and at eight different locations (see Fig. 1e–g) separated by ∼45°. Because subjects performed both center-out and out-center point-to-point movements, for each movement direction (e.g., upward) there were two movements with a different endpoint (e.g., the center-out movement to target 7 and the out-center movement from target 3).

For experiments 1 and 2, each experiment consisted of repeated blocks of trials (five for all subjects with the exception of subjects 4 and 5, who performed four blocks, all performed in a single session) for each one of the six combinations of two planes of movements and three weights (experiment 1) or three forearm postures (experiment 2). Each block included one center-out movement and one out-center movement for each one of the 8 peripheral targets for a total 48 movement conditions per plane. Each trial started after the subject reached the start position, with a “ready” signal followed by a “start” signal (computer-generated tones; 1 s delay), after which the subject was free to choose when to move but was instructed to move fast and accurately to the target and to remain at the target location for at least 1 s (target “hold” period). A trial was considered successful if the subject performed the movement with a duration (defined as the interval in which the endpoint speed was >10% of its maximum speed) <400 ms, stopped within a distance of 10 cm from the target center, and held the final position for at least 1 s. Subjects were given auditory feedback after an unsuccessful movement, and all unsuccessful movements were repeated. Start and target reference positions were defined as the position of the sphere on the subject’s handle, measured during a calibration at the beginning of each experimental session. In experiment 1, these reference positions were measured with the handle sphere close (<1 cm) to the start or target sphere. In experiment 2, the reference positions were on a plane parallel to the plane containing the start and target spheres and shifted a few centimeters from it to allow for movements with the rotated handle. Subjects sat and rested at will between blocks and showed no sign of fatigue throughout the experiment.

In experiment 3, for each plane of movement there were 80 distinct movement conditions: 16 point-to-point movements (1 center-out and 1 out-center for each one of the 8 peripheral targets), 16 reversal movements (1 with the first phase center-out and 1 with the first phase out-center for each one of the 8 peripheral targets), and 48 via-point movements (8 out-center movements for the first phase for each of the 8 peripheral targets, each followed by 1 of 6 center-out movements for the second phase, corresponding to angle of movement from first to second phase of ∼135°, 90°, 45°, −45°, −90°, and −135°). These 80 movement conditions were randomly subdivided in 2 blocks with 40 conditions each, and each condition was repeated 5 times for a total of 10 blocks for each plane of movement. Reversal and via-point trials were considered successful if the movement duration was between 500 and 700 ms. This interval was chosen so that the mean duration of each phase of the reversal and via-point movements would match the typical duration (∼300 ms) (see Table 2) of a point-to-point movement.

Data acquisition.

The position of the handle gripped by the subject and the activity of muscles involved in shoulder and elbow movements were recorded during all experiments. In the experiments performed with apparatus 1, the position and orientation of three markers placed on the subject’s wrist, shoulder, and handle and one marker placed on the experimenter’s (target) handle were acquired with an electromagnetic motion-tracking system (Fastrak; Polhemus, Colchester, VT). Each marker was sampled at 30 Hz with a spatial resolution of <4 mm in each direction, as estimated with a calibration performed within the workspace used in the experiment. In the experiments performed with apparatus 2, the position and orientation of one Fastrak marker in the subject’s handle was recorded at 120 Hz, and the positions of five to seven additional markers were recorded at 120 Hz with an optical motion tracking system (Optotrak 3020; Northern Digital, Waterloo, Ontario, Canada) with an accuracy of <0.1 mm in each direction. Two of these additional markers were placed on the handle, one was placed on the wrist (over the styloid process of the ulna), one was placed on the elbow (over the lateral epicondyle), and one was placed on the upper arm (at the proximal end close to the head of the humerus). In the case of subject 6, who performed experiment 2, three markers instead of a single one were placed on the wrist (two over the styloid process of the radius with different orientations and one lateral to the styloid process of the ulna) to record the wrist position with the forearm in the three different postures. In all experiments, the EMG activity of up to 19 muscles (see Table 1) was recorded with active bipolar surface electrodes (DE 2.1; Delsys, Boston, MA). Each electrode consisted of two parallel silver bars (10 mm spacing) and a differential preamplifier (gain, 10; rms noise, 1.2 μV; common-mode rejection ratio, >80 dB) housed in a compact case (41 × 20 × 5 mm). Each electrode was taped on the muscle belly and connected to an amplifier (Bagnoli-16; Delsys) where the EMG signal was bandpass filtered (20–450 Hz) and amplified (total gain, 1000). For each muscle, correct electrode placement was tested by asking the subject to perform a number of maneuvers involving both free movements and isometric contractions (Kendall et al., 1993) and then observing the expected activation patterns.

Data acquisition and experiment control were performed on a workstation with custom software written in LabView (National Instruments, Austin, TX). EMG data were digitized continuously during each block (1 kHz sampling rate; PCI-6035E; National Instruments). Kinematic data were synchronized with the EMG data by logging the time of each sample (indicated by a synchronization transistor-transistor logic pulse generated by the Fastrak and Optotrak control units) on a counter (100 kHz clock; PCI-6602; National Instruments) synchronized with the EMG acquisition clock. Fastrak data were processed on-line to compute movement time and target accuracy and to provide auditory feedback to the subject regarding unsuccessful trials. The experiment control program logged the times of all relevant behavioral events.

Data analysis.

All analyses were performed with custom software written in Matlab (Mathworks, Natick, MA). For each block, the continuous EMG and tracker recordings were subdivided into segments corresponding to successful movements starting at the time of the ready signal (1 s before the start signal) and ending at the end of the hold period (1 s after target acquisition).

Endpoint kinematics.

We used the position and orientation data from the subject’s handle Fastrak marker and the measured geometric parameters of the handle to compute the position of the sphere attached to the handle used as the endpoint reference position. These position data were low-pass filtered [finite impulse response (FIR) filter; 15 Hz cutoff; zero-phase distortion; Matlab “fir1” and “filtfilt” functions] and differentiated to compute tangential velocity and speed. For each movement, we measured the following: movement onset time, movement end time, movement duration, maximum speed and its time of occurrence, movement vector (in space), and movement direction (in the plane of movement). Movement onset and movement end were identified as the times in which the speed profile superceded 10% of its maximum. Movement duration was defined as the interval between movement onset and movement end. The movement vector was computed as the difference between final and initial positions. For each experimental condition, because the endpoint moved essentially on a plane (either sagittal or frontal), we computed the projection of the movement vector (in space) on the plane of movement identified by the first two principal components of all the movement vectors for that condition. For movements in the sagittal plane, the movement direction was measured as the angle of rotation from the direction of a horizontal, forward movement around the laterally directed axis perpendicular to the movement plane. Thus a movement directed upward had a 90° direction and a movement backward had a 180° direction. For movements in the frontal plane, the movement direction was measured as the angle of rotation from a medially directed horizontal movement around the axis perpendicular to the movement plane directed forward. Thus a movement directed upward had a 90° direction and a movement directed laterally had a 180° direction.

Joint kinematics.

In the experiments performed with apparatus 2, the positions of the shoulder, elbow, and wrist markers were used to estimate shoulder adduction, shoulder flexion, shoulder external rotation, and elbow flexion angles. The coordinates of the Optotrak markers were first rotated in a Cartesian reference frame, with the gravitational acceleration along the z-axis, by means of a calibration of the direction of the gravitational acceleration based on the measurement of two markers attached to a pendulum. The marker position data were then low-pass filtered (fifth-order Butterworth filter; FIR filter; 12 Hz cutoff; zero-phase distortion; Matlab “butter” and filtfilt functions), and the four joint angles were computed (see Fig. 3b) as the angles associated with a sequence of rotations of the shoulder joint (adduction, flexion, and external rotation) and the elbow joint (flexion). Angular velocities and accelerations were computed by numerical differentiation.

EMG preprocessing.

The EMGs for each trial were digitally full-wave rectified, low-pass filtered (FIR filter; 20 Hz cutoff; zero-phase distortion; Matlab fir1 and filtfilt functions), and integrated over 10 ms intervals. All trials in each experimental condition were then aligned on the time of movement onset and averaged (see Fig. 4). In a few cases (Table 1), some muscles showed a change in signal amplitude during the experimental session (indicated by –), likely resulting from a partial detachment of the electrode from the skin, and those muscles were removed from further analysis.

Phasic EMG patterns.

The muscle activities involved in reaching movements in vertical planes are responsible for balancing gravitational forces, for maintaining postural stability, and for accelerating and decelerating the arm to change hand position. Flanders and collaborators (Flanders and Herrmann, 1992; Buneo et al., 1994) have shown that it is possible to dissociate a component of the EMG signal related to holding the arm at specific posture (“tonic”) from one related to the movement (“phasic”). We focused our analysis on the phasic component of the muscle patterns, and we used a subtraction procedure to remove the tonic component from the EMG waveforms. In accordance with Flanders’ data, we modeled the tonic component as a constant muscle activation level before and after the movement and as a linear ramp between two constant levels during the movement. We estimated the activity of each muscle at the initial and final postures by averaging the integrated EMGs from ready time up to 200 ms before movement onset (initial tonic level) and from 200 ms after movement end to the end of the hold period (final tonic level). We then subtracted from the waveform of each muscle the constant initial level from the beginning of the trial up to movement onset, a linear ramp between initial and final levels up to movement end, and the constant final level from movement end to the end of the trial. This procedure was used for all types of reaching movements, including reversal and via-point movements. After subtraction, a phasic waveform could assume negative values, which corresponded to an estimated tonic activation higher than the observed EMG level. Finally, we restricted the time interval considered for further analysis to the samples included between 200 ms before movement onset and 200 ms after movement end, and we normalized the amplitude of each sample for each subject and muscle to the amplitude of the largest sample in that muscle over all conditions.

Muscle synergy extraction.

We modeled the construction of a phasic muscle pattern by the combination of N time-varying muscle synergies as follows: where m(t) is a vector of real numbers, each component of which represents the phasic activation of a specific muscle at time t; w_i(τ) is a vector representing the muscle activations for the i-th synergy at time τ after the synergy onset; t_i is the time of synergy onset; c_i is a non-negative scaling coefficient; and ε(t) is noise.

We used a modified version of the optimization algorithm that we introduced recently (d'Avella et al., 2003) to identify a set of muscle synergies and their onset times and scaling coefficients that reconstructed the entire set of phasic muscle patterns in a dataset with minimum error. The algorithm is initialized by choosing random values for N synergies of a specific duration (500 ms corresponding to 50 samples of integrated EMGs), and it proceeds by iterating the following three steps: (1) given a set of synergies, find the synergy onset times for each trial by a matching pursuit procedure (Mallat and Zhang, 1993; d'Avella and Bizzi, 2005); (2) given a set of synergies and their onset times, find the non-negative scaling coefficients for each trial by non-negative least squares (Matlab function “lsqnonneg”); and (3) given the onset times and the scaling coefficients for all the trials, update the synergies according to the gradient descent on the error function as follows: where w_ij = w_i(τ_j) and θ ⁻(w) sets to zero all positive components of w; s runs over trials (movements); k_s is number of time samples in trial s; t_k is the time of the k-th sample in each trial; and c^s_i and t^s_i are the amplitude and timing coefficients for the i-th synergy in trial s.

The matching pursuit procedure essentially consists of an iterative search for a set of time-shifted synergies that best match the data. For each trial, the instances of each synergy obtained by using all possible time-shifts are compared with the data, and the instance with the highest normalized scalar product is selected and subtracted from the data. This search is then repeated among the remaining synergies by using the data residual.

As a consequence of the tonic activity subtraction, the phasic EMG waveforms can have negative values, and thus we could not use the original error minimization procedure based on a non-negative matrix factorization approach (d'Avella et al., 2003). We used instead a gradient descent minimization to update the synergies (step 3) with an error function that penalizes both inaccurate reconstructions (first term in Eq. 2) and large negative values in the identified synergies (second term). The denominator in the first term represents the squared norm of the entire dataset, and it was used to simplify the setting of the gradient descent parameters across different datasets (i.e., different subjects) by normalizing the contribution of the reconstruction error to the total error. The penalty term for negative synergy components was motivated by two considerations. First, as for other factorization methods, the decomposition of a set of muscle patterns as combinations of time-varying muscle synergies is not necessarily unique. In fact, there might be transformations of the basis elements, i.e., synergies, and associated transformation of the coefficients, which together provide alternative solutions that explain the data as well as the original bases. For example, in the case of factor analysis (Basilevsky, 1994), the basis vectors (factor loadings) can be rotated arbitrarily to obtain a new, equivalent factor structure, and this arbitrariness is usually removed by choosing some additional criterion such as maximizing the variance of the squared normalized factor loadings (varimax rotation). In the case of the time-varying synergy decomposition, because the basis elements are sequences of muscle activation vectors that can be shifted in time by an arbitrary number of samples, the solution is not invariant under a rotation of the synergies; however, a synergy rotation can lead to a new, different solution that approximates the data with a level of accuracy similar to that of the original solution. Thus the penalty term for negative synergy components in the cost function reduces this arbitrariness by constraining the algorithm to identify a set of synergies with minimal negative components among the sets that explain the data equivalently well. Second, large negative activation waveforms do not have a simple physiological interpretation. A negative deflection in the muscle waveform of one synergy can be interpreted as inhibitory influence acting on the motoneuronal pool excitation deriving from other synergies and from postural control; however, the amount of inhibition sufficient to completely silence a muscle, being equal to the amount of excitation driving the tonic activity of the muscle, is typically much smaller than the peak phasic activity. Thus one would expect a maximum negative deflection in the synergy waveforms to be a small fraction of the waveform peak amplitude.

We used a convergence criterion of five consecutive iterations for which the decrease in the error was <10⁻⁴. To minimize the probability of finding local minima, for each N, we repeated the optimization 10 times and selected the solution with the lowest error. In one subject (subject 3), the step size in the gradient descent update rule was set to 2 after testing values of this parameter ranging from 0.5 to 10 and choosing the value associated with the lowest error. The trade-off parameter (λ) between the two terms of Equation 2 was set to 0.05 after testing values from 0 to 0.5 in the same subject and choosing the minimum value for which the minimum of the normalized synergy waveforms was above −0.3, i.e., representing a moderate amount of inhibition.

Measures of synergy reconstruction goodness.

Because the EMG patterns and the residuals of the reconstruction of the patterns by synergy combinations are multivariate time-series, a measure of the goodness of the reconstruction, typically a ratio of two variances, must be defined using a multivariate measure of data variability. We used the “total variation” (Mardia et al., 1979), defined as the trace of the covariance of the muscle activations, to define a multivariate R² measure as follows: where SSE is the sum of the squared residuals, and SST is the sum of the squared residual from the mean activation vector (m̄), i.e., the total variation multiplied by the total number of samples (K = Σ_s k_s). Thus R² represents the fraction of total variation accounted for by the synergy reconstruction. Although this measure is a global indicator of the goodness of reconstruction (operating over all trials), it is also useful to define a goodness of reconstruction measure for individual trials. If we use the same R² value restricted to one trial, we are comparing the residual of that trial with the data variation over that trial, which may vary considerably across trials. Thus the same residual error would have a very different R² value in a trial with a small amplitude modulation of the muscle waveforms and in a trial with a large modulation. We therefore defined a measure of the goodness of the reconstruction of a single trial based on the total variation (across all trials) per unit sample as follows: With this definition, the global R² is the mean of the R_trial² values weighted by their duration.

Selection of the number of synergies.

The number of extracted synergies N is a parameter of the model, and thus we repeated the extraction with a number of synergies ranging from 1 to 8. We selected a specific number of synergies for further analysis according to a procedure based on the dependence of the amount of total variation explained (R²) on N. If the data were generated by combinations of N* synergies and no noise was added, a set with N > N* synergies would not explain more variation in the data than the set with N* synergies, because it already explains 100% of the data variation. When a certain amount of noise is added, the set with N* synergies would explain <100% of the data variation, and the synergy sets with N > N* synergies would explain larger fractions of the data variation; however, the variation captured by the model with N > N* in addition to the variation captured by the model with N = N* represents only variation attributable to noise. Thus, under the assumption that the random variation attributable to noise is smaller that the structured variation attributable to the synergy combinations, we expect the R² curve to change slope at N = N*. Therefore, we examined the R² versus N plot to estimate the correct number of synergies. This procedure is analogous to the “scree test” commonly used for the selection of the number of factors in factor analysis (Cattell, 1966) and has been shown to be adequate for identifying the correct number of non-negative synchronous muscle synergies (Tresch et al., 2006). We also used a quantitative criterion for the selection of the number of synergies. Under the further assumption that for N > N* each of the additional synergies captures an equal amount of noise-generated variation, we expect the R² curve to follow a straight line for N ≥ N*. A similar assumption underlies the Bartlett’s test, another commonly used test for the selection of the number of factors in factor analysis (Bartlett, 1950). We used a linear regression procedure (Cheung et al., 2005) to identify the value of N after which the R² curve is essentially straight. We performed a series of linear regressions, starting from a regression on the entire R² curve and progressively removing the smallest N value from the regression interval. We then compared the mean square residual errors of the different regressions and selected as the optimal number of synergies N* the first N value corresponding to a regression line from N to N^max with a mean square error <10⁻⁴.

Significance of extracted synergies.

We verified that the structure of the identified synergies did not result from a bias built into the extraction method with a simulation. We compared the R² value for the reconstruction of the actual data by the combinations of the identified synergies with the R² value for the reconstruction of structureless simulated data by the combinations of synergies extracted from those simulated data. We generated structureless data with the same empirical amplitude distribution and smoothness of the observed data by randomly reshuffling the samples independently for each muscle and then low-pass filtering the reshuffled samples (10 Hz cutoff). The reshuffling procedure consisted of generating a random permutation of an index vector (running from 1 to the total number of samples; Matlab function “randperm”) for each muscle and then reordering the samples using the permutated index. We then constructed a simulated dataset with the same number of trials and with the same duration of each trial as the observed data. For each observed dataset, we simulated 50 datasets and repeated the same synergy extraction procedure used for the observed data. We estimated the significance (p < 0.05) by computing the 95th percentile of the R² distribution for simulated data.

Synergy fit to point-to-point, reversal, and via-point EMGs.

We tested the robustness of the synergy model by reconstructing the muscle patterns in one dataset with the synergies extracted from a different dataset. To determine the amplitude and timing coefficients that best reconstruct a set of point-to-point EMG patterns given a set of synergies, we ran one iteration of steps 1 and 2 of the synergy extraction algorithm after initializing the synergies to the given set. We used this procedure in the analysis of experiments 1 and 2 for reconstructing the muscle patterns recorded during movement with additional loads and with pronated or supinated forearm postures with the synergies extracted from the pattern for point-to-point movements with the unloaded handle and a neutral forearm posture. We also used this procedure to fit data from individual trials with the synergies extracted from averaged muscle patterns to perform a statistical analysis of synergy modulation across conditions. For each subject, plane of movement, and synergy, we studied the effect of direction (eight levels), load or posture (three levels), and endpoint (two levels) and their interactions on the synergy amplitude coefficients with a three-way ANOVA (Matlab function “anovan”). In the case of reversal and via-point movements, in which there are two distinct kinematic phases, each associated with a peak in the endpoint tangential velocity, we used a slightly different fitting procedure. Assuming that in these more complex reaching movements each synergy may be recruited more than once during each movement, we again ran one iteration of the first two steps of the extraction algorithm to fit the synergies to the novel dataset but we modified step 1 to allow for multiple instances of the same synergy to be recruited in each movement. To prevent excessive temporal overlap of the same synergy, we introduced a refractory period of 100 ms (10 samples) in the matching pursuit procedure used to determine the synergy onset times (d'Avella and Bizzi, 2005). In summary, we reconstructed each phasic EMG pattern observed during reversal or via-point movements with a sequence of two or more instances of each one of the time-varying muscle synergies extracted from point-to-point movements; each instance was independently scaled in amplitude and shifted in time.

Synergy set comparison.

The cosine of the angle between two synergies was used as a measure of their similarity. This was defined as the maximum of the normalized scalar products between the two time-varying synergies shifted by k₁ and k₂ samples over all possible relative delays (k₁ − k₂) [for details, see d'Avella et al. (2003)].

To compare synergies extracted from different subjects, we grouped them using a hierarchical cluster analysis. We used the similarity between a pair of synergies (s_ij), computed with the subset of muscles common to all subjects, to define a distance measure (d_ij = 1 − s_ij) and to create a hierarchical cluster tree from all synergy pairs (Matlab function “linkage” with the “average” distance method, i.e., using as distance between two clusters the average distance between all pairs of objects across the two clusters). We then partitioned the tree with the minimum number of clusters for which there was no more than one synergy from the same subject in each cluster. Finally, in each cluster, we identified the element with the highest mean similarity with all other elements and used this mean similarity to quantify the mean similarity in a cluster. We also defined a similarity index (S) between two synergies, as a value ranging from 1, in the case of identical synergies, to 0, for a similarity at chance level, as follows: where s_data is the similarity between the two synergies and s_chance is the mean similarity between 500 pairs of random synergies. We simulated random synergies by sampling the empirical distribution of the activation amplitude of each muscle in the dataset from which the synergies were extracted, constructing sequences of random data with the same length as the synergy duration, and low-pass filtering the resulting sequences (10 Hz cutoff) to match the smoothness in the actual synergies. We then estimated the significance (p < 0.05) of the similarity index from the distribution of the similarities between random synergies.

Directional tuning of synergy amplitude coefficients.

To characterize the directional dependence of the synergy amplitude coefficients in each load or posture and endpoint condition, we used multiple linear regression (Matlab function “regress”) to fit the following model: where θ is the movement direction, c is the synergy amplitude coefficient, and β₀ is an offset parameter. We then rewrote Equation 6, in terms of amplitude and preferred direction of the cosine tuning function, as follows: The goodness-of-fit was quantified with an R² value, and its significance was tested with an F test. Finally, to compare the angular dispersion of the preferred direction of the synergy amplitude coefficient tuning across conditions with the angular dispersion of the preferred direction of the activation of individual muscles, for each muscle and each trial we computed the maximum integrated EMG activity in a 100 ms window sliding across the entire trial (Flanders et al., 1996), and we fit these data with the same procedure used for the synergy amplitude coefficients. Then, for each plane of movement, synergy, or muscle, we compared the angular deviation (Batschelet, 1981) of the preferred direction across the six different conditions (three loads or postures and two endpoints).

Fit of point-to-point shifted EMGs to reversal and via-point patterns.

We compared the reconstruction of reversal and via-point muscle patterns by the combinations of the muscle synergies extracted from point-to-point patterns with the reconstruction obtained by fitting the muscle patterns of the two point-to-point movement corresponding to the two separate phases of the reversal and via-point movements, each shifted in time to align the tangential velocity peaks. For each subject and plane of movement and for each averaged, phasic reversal and via-point muscle pattern, we first shifted in time the averaged, phasic patterns of point-to-point movement between the start location and intermediate location of the reversal or via-point movement and of the point-to-point movement between the intermediate location and the target location so that the times of peak tangential velocity in the point-to-point movements matched the times of the first and second peak in the tangential velocity of the reversal or via-point movement. We then found two non-negative amplitude scaling coefficients for the two shifted point-to-point patterns that minimized reconstruction error, using a non-negative least squares minimization algorithm (Matlab function “lsqnonneg”).