## Abstract

How the CNS masters the many degrees of freedom of the musculoskeletal system to control goal-directed movements is a long-standing question. We have recently provided support to the hypothesis that the CNS relies on a modular control architecture by showing that the phasic muscle patterns for fast reaching movements in different directions are generated by combinations of a few time-varying muscle synergies: coordinated recruitment of groups of muscles with specific activation profiles. However, natural reaching movements occur at different speeds and require the control of both movement and posture. Thus we have investigated whether muscle synergies also underlie reaching at different speeds as well as the maintenance of stable arm postures. Hand kinematics and shoulder and elbow muscle surface EMGs were recorded in five subjects during reaches to eight targets in the frontal plane at different speeds. We found that the amplitude modulation of three time-invariant synergies captured the variations in the postural muscle patterns at the end of the movement. During movement, three phasic and three tonic time-varying synergies could reconstruct the time-normalized muscle pattern in all conditions. Phasic synergies were modulated in both amplitude and timing by direction and speed. Tonic synergies were modulated only in amplitude by direction. The directional tuning of both types of synergies was well described by a single or a double cosine function. These results suggest that muscle synergies are basic control modules that allow generating the appropriate muscle patterns through simple modulation and combination rules.

## INTRODUCTION

The control of limb movements is challenging because of the kinematic, dynamic, and biomechanical complexity of the musculoskeletal system. Much of this complexity derives from the large number of degrees of freedom in the system, which allows for great flexibility but makes the computations necessary for controlling arbitrary limb movements, in principle, very complicated (Bernstein 1967). How the CNS might perform these computations has been the central topic of many theoretical and experimental investigations and it still remains an open question. One way to simplify the control process is to subdivide it among task-specific control modules, each involving only a few parameters directly selected as a function of task parameters. For example, the task of reaching a target in space, given an initial arm configuration, is specified by a direction in space and a distance (a vector in space). Movement duration could also be a task parameter if temporal constraints are involved. A direct mapping of the task parameters into the appropriate muscle activation patterns could constitute a task-specific feedforward controller adequately solving the control problem for reaching. However, since incalculable numbers of different muscle patterns are potentially associated with different values of even a single task parameter, whether it is possible to establish such a mapping without an explicit computation remains to be demonstrated.

We have recently provided support for the hypothesis that muscle synergies constitute task-specific control modules that generate all muscle patterns necessary for adequately performing reaching movements through simple rules involving just a few parameters (d'Avella et al. 2006). We have shown that the muscle patterns for fast point-to-point movements in different directions on vertical planes, with different loads and forearm postures, as well as for reaching movements with a reversal and through a via-point, are well captured by the combinations of four or five time-varying muscle synergies. By scaling in amplitude and shifting in time these synergies, it has been possible to account for the complex changes in the activation waveforms of individual muscles across movement conditions, reported in previous studies of reaching movements in vertical planes (Buneo et al. 1994; Flanders 1991, 2002; Flanders and Herrmann 1992; Flanders et al. 1996; Soechting and Lacquaniti 1981), as the result of simple modulation and combination rules. However, our study focused on fast movements and on the phasic components of the muscle activation waveforms responsible for accelerating and decelerating the arm. In contrast, natural reaching movements occur at a range of different speeds, depending on environmental constraints. Moreover, the arm is maintained in stable postures before and after a reaching movement, requiring, for an unsupported arm reaching targets in a vertical plane, adequate muscle activation. How a controller based on muscle synergy combinations might generate muscle patterns appropriate for movements with different speeds and for maintaining stable arm postures is at the core of important questions that must be addressed.

Theoretical considerations (Hollerbach and Flash 1982) suggest a straightforward strategy for generating the joint torque profiles appropriate for movements with different speeds. The equations of motion for an articulated arm are invariant for changes in the timescale of the movement if the joint torques are scaled in amplitude according to a simple rule. After separating the total torque into an antigravity component and a dynamic component and after scaling in time both components by a factor *r*, if the dynamic torque component is scaled in amplitude by *r*^{2}, and the amplitude of the antigravity component is not changed, the joint motions will be scaled in time by the same factor *r*. Thus if a torque profile adequate for reaching a given spatial target at one speed is known, a simple scaling rule allows one to generate the torque profiles for reaching that target, along the exact same path, at different speeds, consistent with the invariant dynamic characteristics of reaching movements observed in humans (Soechting and Lacquaniti 1981). Such a feature of the dynamic equations for an articulated arm suggests that the CNS might implement a simple amplitude-scaling mechanism for mapping movement speed into muscle patterns. We hypothesized that such a scaling mechanism is realized by scaling in amplitude and in time a small number of time-varying muscle synergies.

Maintaining the arm in a specific posture requires generating static joint torques to balance the posture-dependent torques arising from gravity and generating the joint stiffnesses appropriate for maintaining stability in the presence of unexpected perturbations. Specifying the muscle activations that satisfy these postural requirements at any arm posture along the path of a reaching movement—because of the multiarticular nature of the arm and the redundancy of the musculoskeletal system—presents similar computational challenges to specifying the phasic muscle pattern necessary for accelerating and decelerating the arm. We then also hypothesized that a combination of a small number of muscle synergies is sufficient to generate the muscle activations necessary for postural maintenance, thus simplifying the postural control problem to the selection of a few synergy combination coefficients.

In this study, we investigated the muscle patterns for reaching in the frontal plane in eight different directions over a wide range of movement speeds. We characterized the organization of both the postural muscle patterns, responsible for maintaining the hand stable at the end of the reaching movement, and the entire time-course of the muscle patterns, before, during, and after the movement, including both phasic and tonic components. The postural muscle patterns have been decomposed as combinations of time-invariant synergies with an iterative nonnegative matrix factorization algorithm (Lee and Seung 1999; Tresch et al. 2006). The reaching patterns have been decomposed as combinations of phasic and tonic time-varying synergies using a modified version of the iterative algorithm that we have previously used for the analysis of phasic patterns of fast movements (d'Avella et al. 2006). We found that simple synergy amplitude and timing modulation rules explain the variation in the muscle patterns across different directions and speeds, providing further support to the hypothesis that muscle synergies are basic building blocks for constructing simple, task-dependent feedforward controllers.

## METHODS

### Experimental setup and protocol

We investigated the changes in the activity patterns in shoulder and arm muscles across movement direction and speed during point-to-point arm reaching movements in the frontal plane. The experimental apparatus, described in detail in our previous report on fast reaching movements (d'Avella et al. 2006), consisted of one central sphere, indicating a start location, eight peripheral spheres, positioned on a circle at 30 cm from the start on a vertical plane at approximately 45° from each other and indicating target locations, and a handle, gripped by the subject and provided with a reference sphere to be displaced from start to target. Start and target spheres (4-cm diameter, made of transparent plastic) could be illuminated from inside by a light-emitting diode (LED) and were supported by a structure that allowed adjustment of the height of the start sphere to match the height of the elbow of a subject standing with the upper arm along the trunk. Five right-handed subjects (four males and one female, age range 19–36 yr) participated in the experiments after giving informed consent. All experimental protocols conformed to the Declaration of Helsinki on the use of human subjects in research.

Subjects performed blocks of reaching movement trials while standing with the start and target spheres in the frontal plane at a distance that allowed positioning the reference sphere on the handle close to the start sphere while keeping the upper arm vertical and the forearm horizontal. Each trial started, when the reference sphere was positioned at the start location, with a *ready* signal followed by a *start* signal (computer-generated tones, 1-s delay) after which subjects were free to choose when to start moving, although they were instructed to reach the target within a given time interval from movement start and to remain at the target location for ≥1 s (target *hold* period). The target sphere LED was turned on at the time of the ready signal and the start sphere LED was turned off at the time of the start signal. To vary the speed of movement, the instructed movement time interval varied across blocks and subjects were informed at the end of the trial whether their movement was too fast or too slow (*fast* and *slow* signals, computer-generated tones, or high and low pitch, respectively). Five different intervals were used spanning movement times from 240 to 820 ms. Targets in all eight directions were presented in pseudorandom order in each block, until at least four (subject 5) or five (subjects 1–4) successful trials were performed in each direction. One subject performed two blocks for each movement time interval (subject 1, total 10 blocks), whereas all other subjects performed three blocks for each movement time interval (total 15 blocks).

### Data collection and preprocessing

We collected kinematic and electromyographical (EMG) data during the experiments. The position and orientation of a marker inserted in the handle gripped by the subjects were recorded using an electromagnetic motion-tracking system (Fastrak, Polhemus, Colchester, VT) at 120 Hz. The position of the reference sphere attached to the handle, the actual arm endpoint, was computed using translation and rotation matrices that were determined with a calibration performed with a second marker attached to the handle in place of the reference sphere. EMG activity from ≤18 muscles (see Table 1) was recorded using active bipolar surface electrodes (DE 2.1, Delsys, Boston, MA), band-pass filtered (20–450 Hz), and amplified (total gain 1,000, Bagnoli-16, Delsys). Correct electrode placement was verified by observing the activation of each muscle during specific movements known to involve it (Kendall et al. 2005). Possible contamination of the EMG recordings by electrical cross talk from adjacent muscles was assessed by performing a cross-correlation analysis between all pairs of channels. The peak of the normalized cross-correlation at lags close to zero was >0.2 in only seven pairs of channels (four for subject 2 and three for subject 3). Because of the difficulty in distinguishing cross talk from synchronous recruitment of motor units in different muscles (Kilner et al. 2002), we did not remove these muscles from the set used for further analyses. However, we verified that the removal of the muscles potentially affected by cross talk did not change any conclusion drawn from those analyses.

Data acquisition and experiment control were performed on a workstation using 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 kinematic sample with a counter (100-kHz clock, PCI-6602, National Instruments) synchronized with the EMG sampling clock.

Kinematic and EMG data were then digitally low-pass filtered (15-Hz cutoff for kinematic data, 20-Hz cutoff for EMG data after rectification; finite-impulse response filter with zero-phase distortion) and segmented into individual trials. The EMG waveforms were further integrated over 10-ms intervals to reduce the size of the data set. Data for all trials, including those with a movement time outside the instructed interval, were visually inspected and trials with either irregular movement or with artifacts on some EMG channel were removed from the data set (subject 1: 7/672; subject 2: 47/766; subject 3: 56/833; subject 4: 3/748; subject 5: 21/612).

### Data analysis

##### ENDPOINT KINEMATICS.

We characterized the kinematics of the endpoint by measuring: movement onset time, movement end time, movement duration, maximum speed and its time of occurrence, and movement direction in the frontal plane. Movement onset and movement end were identified as the times in which the speed profile crossed a threshold equal to 10% of its maximum. Movement duration (or movement time [MT]) was defined as the interval between movement onset and movement end. The movement direction was computed as the angle of rotation of the movement vector on the frontal plane from a vector pointing medially around the axis perpendicular to the movement plane directed forward. Thus a movement directed upward had a 90° direction and a movement directed laterally a 180° direction.

##### TRIAL GROUPING.

For each subject, all trials were ordered according to their movement duration and the trials with a movement duration included from the 5th to the 95th percentiles were subdivided into five groups (see Table 2). By removing the tails of the movement duration distribution of each subject we ensured a better degree of consistency in the kinematic and EMG characteristics within each group. Since the movement distance was the same in all conditions, grouping according to movement duration was equivalent to grouping according to speed.

##### POSTURAL MUSCLE PATTERNS AND POSTURAL SYNERGIES.

To characterize the postural activity of individual muscles before and after the movement, we computed the mean integrated EMG over a 400-ms interval, between 900 and 500 ms before the movement onset and between 500 and 900 ms after movement end. We then studied the effect of the movement parameters on these postural muscle patterns of individual trials with two-way ANOVA (averaged EMG activity of each muscle vs. 8 directions and 5 speed groups) and multiple linear regression (averaged EMG activity of each muscle vs. cos θ, sin θ, and maximum movement speed, where θ is the movement direction angle).

The variations of the postural muscle patterns after movement end across the experimental conditions were characterized by identifying muscle synergies from the mean postural EMG activity, averaged over repetitions in each direction and speed group. Thus for each subject, we computed a set of 40 vectors, each representing the average activity of all recorded muscles in one condition (8 directions × 5 speeds). We then used a nonnegative matrix factorization algorithm (d'Avella and Bizzi 2005; Lee and Seung 2001; Ting and Macpherson 2005; Tresch et al. 1999) to decompose each one of these muscle activity vectors (**m**) as the combination of a unique set of *N* time-invariant synergies (**w**_{i}) multiplied by condition-specific scaling coefficients (*c*_{i}) These muscle synergies capture the coordinated synchronous recruitment of groups of muscles with specific amplitude balances. Thus they capture the spatial structure of the postural muscle patterns. We refer to these synergies as *postural*. To extract a set of *N* synergies, the decomposition iterative algorithm was initialized with random values for synergies and coefficients and it stopped when the reconstruction *R*^{2} value increased by <10^{−4} over 10 consecutive iterations. At each iteration the algorithm performed two steps: *1*) it updated the synergies given the data and the coefficients; *2*) it updated the coefficients given the synergies and the data. Since both synergy and coefficients are constrained to be nonnegative, the updates are performed efficiently using a multiplicative procedure. For each *N*, the extraction was repeated 10 times and only the synergy set with the highest *R*^{2} was retained. We then selected a specific number of synergies, as a compromise between parsimony and accuracy, according to the curve of the reconstruction *R*^{2} as a function of *N*. We chose the number of synergies at which the curve showed a change in slope (a knee), indicating that adding additional synergies did not significantly improve the accuracy of the reconstruction (d'Avella et al. 2006; Tresch et al. 2006).

##### Directional tuning of postural synergy coefficients.

We characterized the directional tuning of the synergy amplitude coefficients with a cosine function. For each subject, synergy, and speed, we performed a linear regression of the synergy coefficients (*c*) on movement direction (θ) (1) To increase statistical power, we performed the regression on individual trial data. To obtain the synergy coefficients for individual trials, we fit the synergy extracted from averaged data on the muscle patterns for individual trials with the same algorithm used for the synergy extraction without updating the synergies. We then expressed the regression coefficients as preferred direction, modulation amplitude, and offset of a cosine function where θ^{PD} = tan^{−1} (B/A) is the preferred direction. We tested the significance of the regression with an *F* test and quantified the goodness of the fit with an *R*^{2} value. We also quantified the variability of the preferred direction across speeds by computing its angular deviation (Batschelet 1981), defined as the square root of 2(1 − *r*), where *r* is the length of the vector resulting from the sum of unit vectors directed as the preferred directions divided by the number of vectors. Finally, to test for dependence of the synergy coefficients for both direction and speed, we also performed the linear regression where sp_{max} is the maximum of the tangential velocity during the movement.

##### Extraction of time-varying synergies from phasic muscle patterns for fast movements.

To study the organization of the phasic muscle patterns across speeds, we first normalized in time the integrated EMG waveforms to equal movement duration. The waveform of each muscle in each trial was aligned to the time of movement onset and resampled, using linear interpolation, with 50 samples per movement duration, from 1 MT before movement onset to 1 MT after movement end. Time-normalized EMG waveforms were then averaged over trials in the same direction and speed group.

We first identified time-varying synergies, as we did previously (d'Avella et al. 2006), from the phasic muscle patterns for the fastest movements. We found a decomposition of the time-varying muscle activation vectors [**m**(*t*)] as a combination of *N* time-varying synergy vectors, or time-varying synergy [**w**(*t*)] each scaled by an amplitude coefficient (*c*_{i}) and recruited with a specific onset delay (*t*_{i}). Thus a time-varying synergy constitutes a collection of waveforms, each one specific for a muscle, not necessarily synchronous. We extracted time-varying synergies from the phasic waveforms. We estimated the phasic waveforms with a subtraction procedure. We first fit a linear ramp, made by a constant EMG level from 1 MT before movement onset, a constant level 1 MT after movement end, and a constant slope segment from movement onset to movement end, to the EMG waveforms to estimate their tonic components (Buneo et al. 1994; Flanders and Herrmann 1992). The constant levels before and after the movements were computed by averaging the EMG waveforms from 1 to 0.5 MT before onset and from 0.5 to 1 MT after offset. The phasic waveforms were then obtained by subtracting the tonic components from the total EMG waveforms. As for the postural synergies, for each subject, we extracted sets with one to eight synergies using an iterative optimization algorithm and we selected the number of synergies according to the reconstruction *R*^{2} curve. The extraction algorithm was initialized to random values for time-varying synergies 75 samples long (1.5 MT) and it stopped when the reconstruction *R*^{2} value increased by <10^{−4} over 10 consecutive iterations. At each iteration the algorithm performed three steps: *1*) it found the onset delay times, for each synergy and movement condition, with an iterative search (matching pursuits; Mallat and Zhang 1993), given the data and the synergies; *2*) it found, for each synergy and movement condition, a nonnegative scaling coefficient by nonnegative least squares, given the data, the synergies, and the onset times; and *3*) it updated the synergies, given the data and the timing and amplitude coefficients for all conditions, by gradient descent of an error function weighting the reconstruction error and the presence of large negative components in the synergy waveforms (d'Avella et al. 2006). For each *N*, the extraction was repeated 10 times and only the synergy set with the highest *R*^{2} was retained.

##### SIMULTANEOUS EXTRACTION OF PHASIC AND TONIC TIME-VARYING SYNERGIES.

Inspection of the EMG waveforms for slow movements (see Fig. 2) indicated that the tonic components were in most cases much larger than the phasic components and that their shape differed across muscles. Thus the subtraction procedure used for the fastest movement, based on fitting the tonic components with a linear ramp between movement onset and end, although providing a reasonable approximation in case of fast movements with large phasic components, was not adequate for all movement speeds. We then decided to estimate the tonic waveforms directly from the data, with the same approach used for the phasic components, that is, identifying a set of phasic and tonic time-varying components whose combinations would match the entire muscle activation waveforms. Moreover, since the tonic components included—before and after the movement—postural EMG activations, we modeled the generation of the tonic waveforms as linear combinations of time-varying muscle synergies, in accordance with the model of postural activations as a combination of time-invariant synergies. Thus we decomposed the averaged EMG data for all directions and speed, as a combination of phasic and tonic time-varying synergies where the superscripts *p* and *t* refer, respectively, to phasic synergies and tonic synergies. The iterative algorithm used for the extraction of phasic synergies (d'Avella et al. 2006) was modified to allow the simultaneous extraction of tonic synergies, defined as time-varying synergies with a fixed timing relationship with the movement kinematics and with slowly varying activation waveforms. Thus there was no timing coefficient associated with the tonic synergy, whose duration was equal to the duration of the time-normalized EMG waveforms, 3 MT, and their activation waveforms were low-pass filtered after each synergy update step during the optimization. We chose a cutoff frequency equal to the inverse of 2 MT, corresponding to the frequency of an oscillatory sinusoidal movement from start to target with the same normalized average speed as that of the actual movements. For each subject, we extracted the same number of phasic synergies as the number of phasic synergies adequately capturing the phasic muscle patterns of the fastest movements and the same number of tonic synergies as the number of postural synergies selected to describe the postural muscle patterns at the end of the movement. We initialized phasic and tonic synergies with random values and iterated the extraction algorithm until the reconstruction *R*^{2} value increased by <10^{−4} over 10 consecutive iterations. The phasic synergy duration was chosen at 1.5 MT. For each subject, we repeated the extraction 10 times and selected the synergy set with the highest *R*^{2} value.

We compared the phasic time-varying synergies extracted from the muscle patterns for all speeds with the time-varying synergies extracted from the phasic patterns for the fastest movements by computing, as a similarity measure between two synergies, the maximum of the normalized scalar product over all possible relative time shifts of the two synergies (d'Avella et al. 2003). We also estimated a chance-level similarity by computing the distribution of similarities between random time-varying synergies generated with the same amplitude distribution as that of the data sets from which the actual synergies were extracted and with similar smoothness (d'Avella et al. 2006). Analogously, we compared the spatial structure of the tonic synergies after movement end with the spatial structure of the postural synergies extracted from the postural muscle patterns. To this aim we computed the normalized scalar product between the vector constructed with the averaged tonic synergy waveform in the interval from 0.5 to 1 MT after movement end and the vector representing a postural synergy. Chance similarity levels were computed generating random synergies with the same amplitude distribution as that of the data.

##### DIRECTIONAL TUNING OF TIME-VARYING SYNERGY AMPLITUDE COEFFICIENTS.

We characterized the directional tuning of the time-varying synergy amplitude coefficients with cosine functions. We fit the directional dependence of the coefficients with both a single cosine function, using the same linear model used for the postural synergies (*Eq. 1*), and with the sum of two nonnegative, offset cosines where A_{i} is the amplitude modulation of the *i*th cosine, B_{i} is its offset, and *θ _{i}^{PD}* is its preferred direction, with the positive part function [

*x*]

^{+}is

*x*for

*x*>0 and 0, otherwise. The use of a threshold nonlinearity is necessary to describe a bimodal cosine tuning since the sum of two cosines is still a cosine (Flanders and Soechting 1990a). To determine the fit parameters, (A B θ

^{PD}) for a single cosine function and (A

_{1}B

_{1}

*θ*A

_{1}^{PD}_{2}B

_{2}

*θ*) for a double cosine function, we used a nonlinear least-square iterative curve-fitting algorithm (Matlab function

_{2}^{PD}*lsqcurvefit*). We initialized the amplitude and offset parameters to half the mean coefficient value and the preferred direction to the values of the largest local maxima of a periodic cubic spline approximating the directional dependence of the coefficients on the movement direction. We performed this nonlinear fit, for each subject, synergy, and speed on individual trial coefficients, determined by fitting the synergies extracted from averaged data on the individual trial data using a single iteration of the first two steps of the synergy extraction algorithm. We quantified the goodness of the fit with an

*R*

^{2}value and the variability of the preferred direction across speeds by computing its angular deviation. For each condition, we fit both a single and a double nonnegative cosine function and we compared the results with the Akaike information criterion (AIC; Akaike 1974) and the Bayesian information criterion (BIC; Schwarz 1978). We then selected the number of cosines according to BIC, which favors models with smaller number of parameters than AIC (Zucchini 2000).

## RESULTS

We investigated the organization of the muscle patterns for reaching in different directions and with different speeds. Subjects moved a reference marker attached to a handle from one central location to eight peripheral targets in the frontal plane with a movement duration ranging, on average across subjects, from 285 to 846 ms. For all movement directions and speeds the endpoint paths were approximately straight and the tangential velocity profiles bell-shaped (Fig. 1). For each subject, trials were subdivided into five groups according to movement duration (MT, Table 2) or, equivalently, given that movement distance was constant, to movement speed. The rectified, low-pass filtered, and integrated EMG waveforms recorded from 17–18 shoulder and arm muscles (Table 1) were aligned on the time of movement onset and averaged across trials in the same direction and speed group. These EMG waveforms showed a tonic activity, responsible for stabilizing arm posture against gravity during movement, and a phasic activity associated with the generation of the forces necessary to accelerate and decelerate the arm (Fig. 2). The postural activity after the movement appeared modulated by movement direction. The phasic activity during the movement changed with movement direction and, for a given direction, increased in amplitude with increasing speed. Our goal was to relate the changes in the activation waveforms of individual muscles with direction and speed to the modulation of a few muscle synergies, the coordinated recruitment of muscle groups.

### Postural muscle synergies

We first analyzed the postural muscle activity recorded before and after the movement. There was no significant effect of movement direction on the EMG activity before movement onset in most cases (79.1%; 68/86 muscles in five subjects; *P* > 0.05; two-way ANOVA; averaged EMG vs. speed group, five levels, and direction, eight levels) but there was a significant effect of speed in half of the cases (51.2%; 44/86 cases; *P* < 0.01). The dependence on speed and direction of the postural activity at the start location is likely due to movement preparation, since trials with the same instructed movement duration range were performed in blocks and the time interval analyzed followed the delivery of the target cue (see methods). The averaged EMG activity after movement end depended, as expected, on movement direction (100%; 86/86 cases; two-way ANOVA; averaged EMG vs. speed group and direction; *P* < 0.01) but also depended, in most cases, on speed (82.6%; 71/86 cases). This dependence on speed is probably related to the control of endpoint stiffness at the target location through an increase in cocontraction for faster movements. In fact, the EMG activity after movement end increased with the movement maximum speed in most cases (68.6%; 59/86 cases; significant positive regression coefficient for maximum speed in the multiple regression of averaged EMG vs. maximum speed and movement direction; *P* < 0.01).

We tested whether the postural muscle patterns recorded after movement end, i.e., with the endpoint at eight different locations in the frontal plane, could be reconstructed by the combination of a small number of muscle synergies. For each subject, we identified sets of one to eight muscle synergies from the postural muscle patterns averaged over trials in the same direction and speed, using a nonnegative matrix factorization algorithm (see methods). Three postural muscle synergies were selected in all five subjects as the number of synergies at which the curve of the reconstruction *R*^{2} had a knee (Fig. 3). The *R*^{2} value of the reconstruction ranged from 0.90 to 0.94 across subjects, indicating that the variations in the postural muscle patterns were well captured by the combination of three synergies.

Each postural synergy expressed a specific balance in the activation of the muscles and was modulated in amplitude across movement conditions. For example, distinctive features of the first postural synergy identified in subject 4 (W_{1}, Fig. 4 *A*) were the strong activations of the short head of biceps brachii (BicShort), anterior deltoid (DeltA), and trapezius, especially the inferior (TrapInf) portion. This synergy was recruited maximally for upward targets, as indicated by the unimodal directional tuning of the synergy amplitude coefficient (C_{1}_{,} Fig. 4*B*). The amplitude coefficient for the upward target was also modulated by the speed of the movement. The structure of the second synergy (W_{2}) was characterized by a strong activation of the long head of biceps brachii (BicLong) and pectoralis major (PectClav and PectInf) and by null activation of superior (TrapSup) and middle (TrapMed) trapezius. This synergy was maximally recruited for medial targets and modulated in the preferred direction by speed as well (C_{2}). Finally, the third synergy (W_{3}) showed a strong activation of the long head of biceps brachii, middle deltoid (DeltM), and no activation of the short head of biceps brachii and the clavicular portion of pectoralis major (PectClav). This synergy was maximally recruited for lateral targets (C_{3}).

The directional tuning of the recruitment amplitude of the postural synergies was well captured by a cosine function. To perform a statistical analysis of the dependence of the synergy amplitude coefficients on the movement direction and speed, we fit the data for individual trials with the synergies extracted from averaged data (see methods). We performed a linear regression of the synergy amplitude coefficients for individual trials, separately for each subject and speed group, on the movement direction (Fig. 5). This regression, corresponding to fitting a cosine function (see methods), was significant for all subjects, speeds, and synergies (75/75 cases; *F* test; *P* < 0.01). The median of the distribution of the *R*^{2} values of all the regressions was 0.84. In comparison, the median on the *R*^{2} values of all the regressions of the EMG activity of individual muscles (*n* = 435, 5 subjects, 16–17 muscles, 5 speeds) was 0.60, significantly lower than the median of the *R*^{2} values for the synergy coefficients (*P* < 10^{−6}; Wilcoxon rank-sum test). The angular deviation of the preferred direction of the cosine tuning across speeds, for each subject and synergy, was very small. The median of the distribution of angular deviation values for all subjects and synergies was 7.9° and the maximum 26.7°, indicating that the directional tuning of the synergy amplitude coefficients did not change substantially with movement speed. Finally, we performed a linear regression of the synergy amplitude coefficients, for each subject, on movement direction and maximum speed. All regressions were significant (15/15; *F* test; *P* < 0.01) and in 8/15 cases there was a significant maximum speed regression coefficient (*P* < 0.01; 6/15 positive; 2/15 negative). Thus some of the synergies were also modulated in amplitude with movement speed.

In summary, the variations of the postural EMG activity in many shoulder and arm muscles at the end of point-to-point reaching movements with different directions and speeds were well captured by the combinations of three postural muscle synergies. The synergy recruitment amplitude was modulated with direction as a cosine function with a preferred direction that did not change with movement speed.

### Phasic and tonic time-varying muscle synergies

We then considered the muscle activity recorded before, during, and after the arm movements. Our goal was to identify phasic synergies (i.e., synergies responsible for accelerating and decelerating the arm during the movement) and tonic synergies (i.e., synergies counteracting gravity and stabilizing posture during movement) that could explain the complex changes observed in the muscle activation waveforms across movement directions and speeds. To compare the muscle patterns for movements with different durations, we normalized in time the muscle waveforms to the movement duration of each trial.

We first determined the number of phasic muscle synergies required to adequately reconstruct the muscle patterns for the fastest movements. For these movements, we estimated the tonic waveform of each muscle by fitting a linear ramp and we subtracted it from its total waveform to obtain a phasic waveform (d'Avella et al. 2006). For each subject, we then extracted sets of one to eight time-varying synergies from the phasic waveforms for movements in the group of trials with the shortest duration in all eight directions, time-normalized to unit movement duration and averaged over trials in the same direction (see methods). Each time-varying synergy spanned a time interval equal to 1.5 MT. We selected the sets with three synergies, in all subjects, as the number of synergies at which the curve of reconstruction *R*^{2} had a clear change in slope (Fig. 6). The *R*^{2} value for three synergies ranged from 0.84 to 0.90 across subjects.

We then identified, for each subject, three phasic and three tonic time-varying synergies from the entire set of EMG waveforms of all muscles for all movement directions and speeds, time-normalized and averaged across repetitions. The optimization algorithm used to extract phasic synergies was modified to allow the identification of additional tonic synergies, defined as time-varying synergies with a fixed delay with respect to movement onset and slowly varying activation waveforms (see methods). We chose the same number of phasic synergies as the number necessary to adequately reconstruct the phasic muscle patterns for the fastest movements and the same number of tonic synergies as the number of postural synergies necessary to capture the variations in the postural muscle patterns after the movement end. The duration of the phasic synergies was set to 1.5 MT and the duration of the tonic synergies to 3 MT, corresponding to the entire duration of the muscle pattern analyzed in each condition. The *R*^{2} value for the reconstruction of the muscle patterns with three tonic and three phasic synergies ranged from 0.86 to 0.89 across subjects. The phasic synergies extracted with this procedure were in most cases similar to the phasic synergies extracted from the phasic EMG waveforms of the fastest movements. For each subject, we matched the three phasic synergies extracted from all speeds to the three synergies extracted from the fast movement only and we compared their similarity with the similarity between random synergies (see methods). In 10 of 15 comparisons, the similarity was significantly higher than that by chance (*P* < 0.05). The mean similarity between the two sets of synergies was 0.78 ± 0.13 (SD). Similarly, in most cases the vectors obtained averaging the tonic synergy waveforms over the final 0.5 MT interval, representing the spatial structure of the tonic synergies, were similar to the postural synergy vectors extracted from the postural EMG activity at the movement end. In 9 of 15 comparisons, the similarity was significantly higher than that by chance (*P* < 0.05). The mean similarity between these pairs of synergies was 0.90 ± 0.16 (SD). In sum, the organization of the muscle patterns for reaching movements in different directions and speeds was well captured by three phasic and three tonic time-varying synergies extracted simultaneously from the entire set of EMG waveforms of all conditions.

The spatiotemporal organization of the phasic and tonic synergies showed distinctive features. In each synergy different muscles were recruited with waveforms of different shapes and different amplitudes. The synergies identified in subject 4 (Fig. 7) illustrate many of the features common to all subjects. Some of the waveforms of the phasic synergies had a single positive peak, whereas other waveforms had more complex shapes with two or three positive and negative peaks. In some cases the peaks in the waveforms of different muscles were roughly aligned but different muscle groups peaked at different times. For example, W_{1}^{p} (Fig. 7, *first column*) showed an early synchronous burst in elbow flexors (BicShort, BicLong, Brac, and BrRad), shoulder flexors (DeltA, PectClav, and PectInf), and two portions of the trapezius muscles (TrapSup and TrapInf) followed by a second synchronous burst of elbow extensors (all three heads of triceps brachii: TrLat, TrLong, and TrMed) and latissimus dorsi (LatDors). Similarly, in W_{3}^{p} (*third column*), the waveforms of middle and posterior deltoid (DeltM and DeltP), middle and inferior trapezius (TrapMed and TrapInf), and infraspinatus (InfraSp) were recruited synchronously. In other cases the timing of the different muscle waveforms varied more gradually. For example, in W_{2}^{p} (*second column*), there was an initial negative (deactivation) peak in anterior deltoid (DeltA) and the two portions of pectorialis, followed, in sequence, by a positive peak in DeltP, LatDors, and teres major (TeresMaj); a negative peak in TrapMed and TrapInf; a large burst in the three heads of triceps; a negative peak in the two heads of biceps; a positive peak in DeltA, PectClav, and PectInf; and, finally, a positive peak in TrapMed, TrapInf, and InfraSp. Thus in general, the phasic synergies showed a rather complex spatiotemporal organization with specific waveform amplitude and timing for each synergy. Some features of this spatiotemporal organization, such as the synchronous recruitment of shoulder and elbow flexors in the initial phase of W_{1}^{p}, are in accordance with some of the well-established kinematic and dynamic characteristics of reaching movements, such as the tight coupling of shoulder and elbow motions and torques (Lacquaniti et al. 1986; Soechting and Lacquaniti 1981). Finally, the different tonic synergies (Fig. 7, *columns 4*–*6*) were characterized by specific amplitude balances among muscles, whereas the shapes of the waveforms, constrained to contain only low frequencies, were mostly slowly varying increasing or decreasing ramps.

### Muscle pattern reconstruction by synergy modulation

The example of the reconstruction of the integrated, time-normalized, and averaged EMG patterns for upward and downward movements at different speeds of subject 4 (Fig. 8) illustrates a remarkable finding of our analysis. The EMG patterns for movements in a given direction at different speeds, once decomposed into phasic and tonic synergies, appear generated by a simple rule. A specific balance of tonic synergies was combined with a specific balance of phasic synergies, with amplitude increasing with movement speed and invariant timing. Even for the slowest movements, with EMG waveforms made up for the most part by the tonic activity, the small phasic components were well captured by the same balance of phasic synergy waveform highly recruited in the fastest movements. For example, the phasic Brac and TrapMed waveforms during the slowest upward movements (Fig. 8, *column 5*, peaking close to movement onset time, first dashed vertical line), despite their small amplitude, had a shape and a timing remarkably similar to those for the fastest movements (Fig. 8, *column 1*).

The modulation of the recruitment of the phasic and tonic synergies across all movement directions and speeds fully characterized the rule for the generation of the muscle patterns as synergy combinations. Phasic synergies were modulated both in amplitude and in timing by direction and only in amplitude by speed. Tonic synergies were modulated in amplitude by movement direction but not substantially by movement speed. This modulation pattern is illustrated by the polar plot of the synergy amplitude and timing coefficients as a function of movement direction for different speeds of subject 4 (Fig. 9). The amplitude coefficients (Fig. 9*A*) for both phasic (*columns 1*–*3*) and tonic (*columns 4*–*6*) synergies had a similar directional tuning for all speeds, but the maximum of the tuning curves depended on speed (more saturated colors corresponding to higher speeds) only for phasic synergies. Moreover, the onset time (Fig. 9*B*) of the phasic synergies strongly depended on the movement direction, although it had a subtle yet significant dependence on direction and speed in the range of directions in which each synergy had a large amplitude coefficient (Fig. 9*B*, *bottom*, bar plots).

### Directional tuning of synergy coefficients

The directional tuning of the amplitude coefficients was captured in most cases by either a single or a double cosine function. We performed a statistical analysis of the synergy coefficients for individual trials. We determined the coefficients by fitting the synergies extracted from averaged data on the individual trial data (see methods). We first performed a linear regression of the synergy amplitude coefficients, separately for each subject and speed, on the movement direction. Such linear regression amounts to fitting a single cosine with an offset to the data. The regression was significant in all cases for tonic synergies (75 cases: 5 subjects × 5 speeds × 3 synergies; *F* test; *P* < 0.01) and in 71 of 75 cases for phasic synergies. The median value of the distribution of the *R*^{2} values of all the regressions for tonic synergies (*n* = 75) was 0.82, whereas the median value for phasic synergies (*n* = 75) was substantially lower, 0.41. Visual inspection of the tuning curves indicated that in many cases the tuning curve for phasic synergies had a single large peak but also a smaller peak roughly in the opposite direction. We then performed a nonlinear curve fit of the amplitude coefficients of the phasic synergies with either a single cosine or the sum of two nonnegative offset cosine functions (see methods). The median of the distribution of the *R*^{2} values for single nonnegative cosine fit was 0.46 and the median for double cosines was 0.68. In almost all cases, a double cosine better described the directional tuning of the phasic synergy than a single cosine, even taking into account the larger number of fit parameters of the double cosine fit (66/75 cases according to the BIC; 70/75 cases according to AIC; see methods). In comparison, the median of the *R*^{2} values of the linear regression of the amplitude coefficients of the postural synergies extracted from the postural EMG activity after the movement (see earlier text) was 0.84, not significantly different from the median value for the tonic time-varying synergy amplitude coefficients (Wilcoxon rank-sum test, *P* = 0.24). In contrast, the median of the *R*^{2} values of the cosine fit of EMG activity of individual muscles after the movements (see earlier text) was 0.60, significantly lower than the median value for the tonic time-varying synergies (*P* < 10^{−6}). We also compared the *R*^{2} values for the fits with a single or a double cosine function of the amplitude coefficients of the phasic synergies with the *R*^{2} values of the fits with single or double cosines of the largest burst of EMG activity, integrated over a 0.5 MT interval, of individual muscles. For both synergies and muscles we chose one or two cosines according to the BIC criterion. The median of the *R*^{2} values for the synergies was 0.67, significantly higher than the median value for the individual muscles, 0.48 (*P* < 10^{−6}). Finally, in most cases, when the directional tuning of a phasic synergy was well characterized by a double cosine, the onset time of the synergy for movements in the direction of peak amplitude of one cosine was different from the onset in the direction of peak amplitude of the other cosine. The mean difference, in absolute value, over all conditions and with amplitude tuning best described by a double cosine (according to BIC, *n* = 66), between the mean onset time of all trials with a movement direction within 30° from the preferred direction of the first cosine and the mean onset time of all trials with a movement direction within 30° from the preferred direction of the second cosine was 0.35 MT, significantly different from 0 (*t*-test, *P* < 10^{−6}).

Since large variations of the synergy onset time occurred in most cases for directions with a large angular deviation from the synergy preferred direction, whereas in directions close to the preferred direction the timing modulation was more subtle, we assessed the role of synergy timing modulation by evaluating the effect of fixing the onset time on the reconstruction *R*^{2} values. For each subject and synergy, we fixed the onset time for all conditions as the mean onset time weighted by the synergy amplitude coefficient. We found that the *R*^{2} values on average decreased by 11.2 ± 3.9% (Δ*R*^{2}/*R*^{2}, mean ± SD), which indicated a small but significant contribution of the timing parameter in the muscle pattern reconstruction. In addition, we performed a two-way ANOVA to test the effect of direction and speed on synergy onset time, for individual trials in directions within 60° of the synergy preferred direction of the single or first cosine. We found a significant (*P* < 0.01) main effect of speed in 13/15 cases (3 phasic synergies × 5 subjects) and direction in 14/15 cases.

### Synergy amplitude modulation with speed

The preferred direction of the cosine functions describing the directional dependence of the amplitude coefficient of tonic and phasic synergies did not vary substantially with movement speed. For tonic synergies (Fig. 10 *A*), the angular deviation of the preferred directions of the single cosine fit across speeds, for all subjects and synergies (*n* = 15), ranged from 0.6 to 6.5°, with a median value of 3.7°. In contrast, the angular deviation of a cosine fit of the EMG activity at the end of the movement, across speeds, for all subjects and muscle (*n* = 87) ranged from 1.0 to 69.5°, with a median value of 6.0°, significantly higher than that for the tonic synergy coefficients (Wilcoxon rank-sum test, *P* = 0.0071). For phasic synergies (Fig. 11 *A*), the angular deviation of the preferred direction of the first cosine (see methods) ranged from 0.5 to 41.4°, with a median of 3.7°. In contrast, the angular deviation of the largest integrated EMG burst ranged from 1.3 to 71.5°, with a median of 7.7°, significantly higher (*P* = 0.014) than the median for the amplitude coefficients of the phasic time-varying synergies.

The amplitude of the cosine functions depended on speed in most cases for phasic synergies (Fig. 11*B*) and in a few cases for tonic synergies (Fig. 10*B*). In particular, the maximum amplitude of the cosine functions for phasic synergies appeared to increase more than linearly with the speed. To evaluate the effect of speed on the synergy amplitude coefficients, we performed a linear regression, for each subject and synergy, of the logarithm of the amplitude coefficients for all trials, in the direction with the highest mean coefficient, on the logarithm of the maximum speed (Fig. 12). For phasic synergies, all regressions were highly significant (*n* = 15, *F* test, *P* < 10^{−6}), the median *R*^{2} value was 0.84 (range 0.66–0.92), and the median regression slope was 2.03 (range 1.40–2.73). For tonic synergies, 7 of 15 regressions were not significant (*P* > 0.05) and the median slope of the 8 significant regressions was 0.32 (range 0.12–0.69), whereas their median *R*^{2} value was 0.50 (range 0.10–0.58). Thus the amplitude of the phasic synergies in the direction of maximal recruitment increased roughly as the square of the maximum speed. In contrast, in many cases, the amplitude of the tonic synergies did not vary significantly with speed and, when it did, it increased with speed less than linearly.

## DISCUSSION

Analysis of the muscle patterns recorded during reaching movements in different directions and with different speeds led to a number of findings that support our hypothesis of a modular control architecture based on a direct mapping of task parameters into muscle synergy recruitment parameters. First, the complex changes observed in the muscle activation waveforms of several shoulder and elbow muscles across movement directions and speeds are captured by scaling in time, scaling in amplitude, and shifting in time a few time-varying muscle synergies. Each individual muscle waveform is reconstructed by the sum of component waveforms belonging to different synergies and its changes in shape and timing across movement conditions are systematically related, through the synergy structure, to the changes of the other muscle waveforms. Thus simple rules allow one to specify the entire muscle pattern by selecting a small number of synergy recruitment parameters. Second, the muscle patterns for reaching are generated by modulating two different types of time-varying synergies: phasic and tonic. Phasic synergies are modulated both in amplitude and in onset timing with movement direction and speed. Tonic synergies are modulated in amplitude with movement direction and not substantially with speed. The directional modulation of the amplitude coefficients of the tonic synergies is well described by a cosine function, whereas the directional modulation of the amplitude coefficients of the phasic synergies is captured, in most cases, by the sum of two nonnegative cosine functions. In those cases, synergies are recruited, in the preferred directions of different cosines, with different onset times. The dependence of the amplitude coefficients on the phasic synergies on the movement speed is more than linear and, on average, close to quadratic. Thus the rules mapping movement direction to muscle patterns can be expressed through simple synergy coefficient tuning functions. Third, the muscle activity recorded after the movement, responsible for maintaining a stable posture at the target location, is generated by the combination of a few postural, time-invariant synergies. The spatial structure of these synergies is similar to the spatial structure of the tonic time-varying synergies at the end of the movement, compatible with the idea that tonic synergies are responsible for antigravity and postural control and phasic synergies are responsible for overcoming inertia and accelerate and decelerate the arm.

### Synergy model formulation and identification

We have used two types of synergy models to characterize the coordinated recruitment of groups of muscles underlying the control of reaching. Postural muscle activity after the movement end was modeled as the combination of time-invariant synergies, whereas muscle activity before, during, and after the movement was modeled as the combination of time-varying synergies. A time-invariant synergy expresses a relationship in the activations of different muscles that does not depend on time. This was a natural choice for capturing the muscle patterns during the maintenance of a static arm posture. In contrast, a time-varying synergy expresses a spatiotemporal (across muscles and time) relationship in the muscle activations. Such a relationship is described by a collection of activation waveforms that can be different across muscles and can capture both synchronous and asynchronous muscle recruitment. Thus a time-varying muscle pattern can be described by just two recruitment parameters—amplitude and onset time—for each synergy. However, a time-varying pattern can also be described by time-varying combinations of time-invariant synergies. We chose to characterize the regularities in the time-varying patterns before, during, and after the movement using a time-varying model because of the parsimony of this model with respect to the time-invariant model. In fact, the latter requires specifying the entire activation waveform for each synergy to describe a time-varying muscle pattern (d'Avella et al. 2006). Moreover, when we extracted time-invariant synergies from the same time-normalized muscle patterns used to extract phasic and tonic time-varying synergies (data not shown), we found that it was difficult to distinguish phasic and tonic time-invariant synergies because most synergies showed both phasic and tonic components in their activation waveforms. This is not surprising since, with a time-invariant model, any temporal regularity in the muscle patterns can only be captured by the synergy activation waveforms, and time-invariant synergies can be separated into phasic and tonic only if the phasic and tonic muscle activities are generated by orthogonal synergies, a condition that does not appear to occur in practice. In contrast, tonic time-varying synergies can be separated from phasic synergies according to their spatiotemporal structure by requiring that their waveforms have frequencies not higher than the movement frequency and that their onset is fixed with respect to movement onset. In summary, we believe that parsimony and the possibility of separating phasic and tonic synergies according to their spatiotemporal structure support our choice of a time-varying synergy model. However, our data and our approach based on characterizing the regularities of the muscle patterns during unperturbed reaching movements do not allow us to distinguish experimentally between a time-varying and a time-invariant model. Such a distinction might be possible by testing the predictions of the two models on the effect of experimental manipulations. Although the specific conclusions may not be relevant for the control of reaching movements in primates, such an approach has recently been pursued on spinalized frogs by perturbing proprioceptive feedback (Kargo and Giszter 2008).

Time-invariant postural synergies were identified with a nonnegative matrix factorization algorithm (Lee and Seung 1999, 2001; Tresch et al. 1999, 2006). This algorithm iteratively finds a set of nonnegative synergy vectors, common to all conditions, and a set of nonnegative combination coefficients, specific to each condition, that minimize the reconstruction error. This decomposition of muscle patterns into time-invariant synergy vectors and time-varying combinations of coefficients is related, and in many cases similar (Tresch et al. 2006), to that obtained with other algorithms such as principal component analysis (PCA) (Mardia et al. 1979), factor analysis (Basilevsky 1994), and independent component analysis (Bell and Sejnowski 1995; Makeig et al. 1997). These algorithms have been used to analyze the muscle patterns recorded in different species and in different behaviors in many studies (Cheung et al. 2005; Hart and Giszter 2004; Ivanenko et al. 2003, 2004; Krishnamoorthy et al. 2003; Olree and Vaughan 1995; Patla 1985; Saltiel et al. 2001; Ting and Macpherson 2005; Torres-Oviedo et al. 2006; Tresch et al. 1999; Weiss and Flanders 2004). In contrast, time-varying synergies were identified with an iterative optimization algorithm developed specifically for identifying spatiotemporal components in the muscle patterns (d'Avella and Tresch 2002) and, to date, used only in a few studies (d'Avella and Bizzi 2005; d'Avella et al. 2003, 2006). A time-varying decomposition approach based on PCA has been recently used (Klein Breteler et al. 2007), although this approach does not allow for independent time shift of the muscle synergies.

Synergy models are useful only if the regularities in the motor output expressed by the synergies provide information on the functional architecture of the controller. Although the assumption that regularities in the output of a system reflect its internal organization is common in science, there are some concerns about the application of these approaches to the analysis of biological motor systems. First, it must be ensured that the regularities are not due to poor sampling of the motor output. For example, in examining movements in only one condition one might erroneously take a stereotypical muscle pattern as a centrally organized synergy simply because there is not enough variability in the task specification. By sampling eight reaching directions and a broad range of movement speeds, we believe we have found an adequate compromise between task variability and duration of the experimental sessions. Muscle patterns change significantly across directions and speeds (see Fig. 2), thus indicating that the structure captured by the synergies is not simply due to a lack of variability in the data. Second, if the synergies are organized by the CNS and not mere statistical description of the regularities in a specific data set, they should be able to capture the regularities in the muscle patterns recorded in conditions not included in the data used for the synergy identification. For reaching movements, we have previously tested the robustness of the time-varying synergies extracted from point-to-point movements by showing that they can adequately reconstruct the muscle patterns in novel dynamic and postural conditions as well as patterns underlying more complex via-point and reversal movements (d'Avella et al. 2006). For postural control, Torres-Oviedo and collaborators (2006) have shown that the time-invariant synergies extracted from the hindlimb muscle patterns of cats during support surface translations in different directions at a normal fore–hind paw distance capture the muscle patterns recorded at different fore–hind paw distances as well as during support surface multidirectional rotations. Finally, a direct validation of the synergy models, as for any model based on the parsimonious description of the regularities of the motor output, can come only from testing causal predictions of the model, such as the effects of manipulations of the controller. One exciting possibility is to perform such a test in the context of motor learning and skill acquisition.

### Number and structure of synergies

For all subjects, three phasic time-varying synergies were selected to characterize the organization of the phasic muscle patterns across movement directions and speed according to the curve of the reconstruction *R*^{2} for time-varying synergies extracted from the phasic muscle patterns of the fastest movements (Fig. 6), as in our previous study of fast reaching movements (d'Avella et al. 2006). The point on the curve at which the slope markedly decreased was taken as an indication of the correct number of synergies. The rationale for this selection is that, until all synergies capturing structural variation in the data are included, the reconstruction *R*^{2} increases rapidly and that, after that point, only a small amount of variation due to noise is captured by additional synergies. Similarly, three time-varying tonic synergies were selected because the curve of the reconstruction *R*^{2} for postural time-invariant synergies showed a clear change in slope (Fig. 3).

We found that the structure of the time-varying synergies extracted from the phasic muscle patterns of the fastest movements was in most cases similar to the structure of the phasic time-varying synergies extracted together with the tonic synergies from the entire set of muscle patterns. This result suggests that the structure of the phasic synergies is robust for changes in speed, supporting their role as basic modules for the generation of dynamic control signals. Moreover, the extraction of time-varying synergies from the subset of the data collected in our previous study (d'Avella et al. 2006) containing the same movement conditions of the fast movements of the present study (center-out point-to-point movements to eight targets in the frontal plane with an unloaded handle and neutral forearm posture) resulted in the selection of the same number of three synergies for all nine subjects of that study. In 125 of 135 (92.6%) comparisons between the structure of the three best-matching pairs of synergies in the sets of synergies extracted from the phasic muscle patterns of the nine subjects of that study and from the five subjects of the present study, the synergies were significantly more similar than by chance (*P* < 0.05, with the procedure described in methods for comparing time-varying synergies). Thus it appears that the difference between the number of synergies selected in our previous study (four to five) and that of the present study (three), depends on the larger number of experimental conditions of our previous study, which was nonetheless restricted to fast movements only, including movements in both frontal and sagittal planes and in both center-out and out-center orientations. In fact, the structure of each one of the three phasic synergies identified in this study (Fig. 7) appears related to the structure of five synergies extracted from fast movements in our previous study (see Figs. 6 and 10 in d'Avella et al. 2006) through simple combinations. For example, the first synergy in this study (W_{1}^{p}) is roughly matched by a combination of the first and the third synergy of our previous study (W_{1} and W_{3}). Consistently, all three synergies have a similar directional tuning on the frontal plane (Fig. 9 in both papers), whereas the two synergies extracted in our previous paper have different tuning on the sagittal plane. In sum, extracting synergies from a variable but not exhaustive movement repertoire might result in a partial resolution of the fundamental modular unit when these units are recruited in a correlated fashion across the conditions tested. Stronger evidence for a synergistic organization in human arm movements will therefore come from an investigation of a broad range of natural movements, as done previously for human locomotion (Ivanenko et al. 2005) and in lower vertebrates (d'Avella and Bizzi 2005).

### Synergy modulation with direction and speed

The directional modulation of the synergy amplitude coefficients is characterized by a cosine tuning for postural (Fig. 5) and tonic (Fig. 10) synergies and by a single or double cosine tuning for phasic synergies (Fig. 11). Cosine tuning is characteristic of motor cortical neurons (Caminiti et al. 1991; Georgopoulos et al. 1982) and might represent an optimal strategy for force generation with redundant actuators in the presence of signal-dependent motor noise (Todorov 2002) or a general tuning property of motor variables encoding time derivatives of a position-dependent function (Zhang and Sejnowski 1999), such as muscle shortening rates (Mussa-Ivaldi 1988). At the muscle level, cosine tuning might also derive from a minimum-effort criterion (Fagg et al. 2002). Cosine tuning has been observed in the muscle activations during step-tracking movements of the wrist in different directions (Hoffman and Strick 1999). However, the tuning observed in shoulder and elbow muscles during reaching movements (Flanders et al. 1996) or isometric force production (Flanders and Soechting 1990b; Pellegrini and Flanders 1996) is described by more complex functions than a single cosine function or even a double cosine function. A key issue here is the time interval considered for estimating the directional tuning. When it is possible to clearly separate an *agonist* and an *antagonist* activation time interval, as in most cases for the wrist muscles (Hoffman and Strick 1999), the rectified and integrated EMGs show a single cosine tuning. However, when the timing of muscle activation gradually shifts with movement direction and the directional tuning is computed integrating the EMGs over the time interval with the largest activity (Flanders et al. 1996), the tuning functions have multiple peaks. Thus a possible explanation of these differences in directional tuning is that independently modulated components of the EMG activation waveforms may be associated with a single cosine tuning but, because these components can be shifted in time and scaled in amplitude differently across direction, the combined tuning can be quite complex and can vary over the time course of the movement. This explanation is in accordance with the existence of fixed relationships among these individual muscle components expressed by muscle synergies and with our observation of single and double cosine tuning of the synergy amplitude coefficients. In particular, when the phasic synergies are recruited according to a double cosine, the preferred directions of the two cosine functions are associated with different onset times, suggesting that each cosine tuning function is associated with a different biomechanical role in the control of the movement.

Phasic synergies are strongly modulated in amplitude with movement speed (Figs. 11 and 12), whereas tonic synergies, in most cases, are not (Fig. 10). As mentioned in the introduction, scaling in amplitude and in time the dynamic torques, by an amplitude scaling factor equal to the square of the time-scaling factor, and in time only the antigravity torques generates invariant endpoint paths (Atkeson and Hollerbach 1985; Hollerbach and Flash 1982). Since our EMG data have been time-normalized to equal movement time, amplitude scaling of the phasic synergies with speed suggests that these synergies are involved in the generation of the appropriate dynamic torques for controlling movement trajectories with invariant paths. Similarly, the lack of significant amplitude modulation of the tonic synergies with speed suggests that they are mainly involved with the generation of the appropriate antigravity torques along the movement path.

Our results on the amplitude scaling of phasic and tonic time-varying synergies with speed are compatible and extend previous results on the scaling of individual muscle waveforms based on the normalization of the EMG time base (Flanders and Herrmann 1992). After scaling the time base of the EMG waveforms and endpoint kinematics so that the movement times at different speeds are equal, the waveforms associated with phasic synergies scale in amplitude with speed but the waveforms associated with tonic synergies do not. However, when individual muscle activation waveforms are not time-normalized, the phasic components of individual muscle waveforms appear to correlate best with templates whose time base is scaled in time differently for different muscles (Buneo et al. 1994; Flanders 2002). For example, during reaching movements in one direction on the sagittal plane, Buneo and collaborators reported that the time base of the phasic component of anterior deltoid scaled with movement time more than the time base of the phasic component of biceps brachii. We found that the waveforms of most muscles, including anterior deltoid and biceps brachii (see Fig. 8), are well captured by the combination of the synergy components after time normalizing to equal movement time, i.e., with the same time base factor for all muscles. These observations are not incompatible. Muscle waveforms generated by combining different components, each scaled in time by the same factor but each shifted in time and scaled in amplitude independently, might correlate best with a template whose time base is scaled by a factor different from the common scaling factor of the components. However, the synergy model used in this analysis, which assumes equal scaling for all synergies and captures a time-scaling mechanism common to all synergies, would miss any additional synergy-specific scaling. Thus it will be interesting to refine our model to allow for independent time-scaling of the synergies. To this aim, we are developing a novel synergy extraction algorithm capable of identifying a synergy timescale in addition to the amplitude scale and time shift identified by the current algorithm.

### Postural and tonic synergies

Our results suggest that muscle synergies underlie both the control of movement and the control of posture. Both control tasks must deal with the complexity of a multiarticular arm and both tasks might benefit from a low-dimensional representation of the musculoskeletal characteristics provided by muscle synergies. Maintaining a stable static arm posture requires coordinating the activity of many muscles spanning shoulder and elbow joints to generate the appropriate net torques and endpoint stiffness (Perreault et al. 2004). Postural time-invariant synergies (Fig. 4*A*) might allow one to find an approximated solution to the problem of finding an arm posture with the hand in a desired position in space and the appropriate tonic muscle activations for maintaining that posture against gravity and against unexpected postural perturbations by selecting a few synergy amplitude coefficients. Thus postural synergies might allow one to implement a direct mapping of high-level postural goals into muscle patterns.

Movement involves transitions between static postures and requires appropriate muscle torques and arm impedance (Burdet et al. 2001) for balancing gravity and maintaining stability in each intermediate posture as well as muscle torques for accelerating and decelerating the joints. Tonic and phasic time-varying muscle synergies appear naturally related to the generation of antigravity and dynamic torques, respectively. Such a decomposition of the torque generation process might be a consequence of the different scaling rules for the two torque components necessary for achieving speed-invariant endpoint paths that the CNS might have exploited to simplify its implementation of a controller.

### Neural implementation of a synergy-based controller

The neural control of arm movements is implemented in a distributed cortical and subcortical network involving many different structures of the CNS. However, the anatomical organization and the physiological characteristics of the motor cortex suggest that it might have a specific role in the neural implementation of a controller relying on muscle synergy combinations. The divergent connectivity patterns of the axons of corticospinal neurons in the motor cortex (Schieber 2001), making synaptic contacts onto spinal interneurons (Shinoda et al. 1979, 1981) and motoneurons (Cheney and Fetz 1985; Fetz and Cheney 1980; McKiernan et al. 1998) across several spinal segments, is an adequate anatomical substrate for the organization of muscle synergies. Furthermore, the extensive horizontal, intrinsic axon collaterals in the motor cortex (Huntley and Jones 1991) together with the widespread and overlapping representation of individual muscles on the cortical surface (Rathelot and Strick 2006) might underlie the spatial and temporal patterning of the activity across the population of corticospinal neurons. In humans, focal transcranial magnetic stimulation of the motor cortex shows increased excitability of elbow and wrist muscles during a reaching movement involving coactivation of shoulder muscles, suggesting the involvement of common cortical circuits in the synergistic recruitment of these muscles (Devanne et al. 2002). Finally, recent experimental results in anesthetized cats indicate that the EMG outputs of two cortical points simultaneously stimulated are additive (Ethier et al. 2006, 2007). These results suggest that the corticospinal circuitry is capable of implementing a physiological mechanism allowing for linear summation of the muscle activations associated with the neural representation of different synergies. Thus we speculate that a time-varying muscle synergy might be driven by the dynamic activation of a specific population of neurons in the motor cortex. The spatial structure of the synergy might be determined by the selective recruitment of spinal interneuronal and motoneuronal populations by the set of corticospinal neurons active at each specific time in the time course of the synergy. As time flows, the activity distribution in the cortical population will evolve along an attractor trajectory in neural space, thus determining the spatiotemporal structure of the synergy.

### Conclusions

Our results suggest that the combination of muscle synergies is a general strategy that the CNS uses for implementing the sensorimotor transformations necessary for reaching. A direct mapping of goals and initial states into synergy modulation coefficients might allow for constructing an inverse internal model without an explicit representation of the dynamic behavior of the musculoskeletal system. Such a mapping might exploit the structure of the equation of motion of the arm to generate the torques appropriate for controlling reaching movements with invariant kinematic features through simple scaling rules expressed as the amplitude and timing modulation of a small number of time-varying synergies. Thus the knowledge of the dynamics of the system would be implicitly incorporated in the structure of the synergies. Such a control architecture—based on a low-dimensional representation of motor output provided by the synergies and low-dimensional mapping of task-relevant sensory information into synergy recruitment parameters—might also simplify motor learning and motor adaptation.

## GRANTS

This work was supported by the Italian Ministry of Health, Fondo per gli Investimenti della Ricerca di Base and Progetti di Ricerca di Interesse Nazionale from the Italian Ministry of University and Research and a Disturbi del Controllo Motorio e Cardiorespiratorio grant from the Italian Space Agency.

## Acknowledgments

We thank Y. Ivanenko and V. C. K. Cheung for comments on the manuscript.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2008 by the American Physiological Society