Dimensionality of joint torques and muscle patterns for reaching

Introduction

How the central nervous system (CNS) coordinates a large number of muscles to generate complex motor behaviors is an open question. The dynamic complexity of the skeletal system with its many degrees of freedom (DoF), the versatility of the motor system, capable of accomplishing many different tasks, and the redundancy and non-linearity of the muscular apparatus all pose challenging control problems. A modular architecture has been proposed as a way for the CNS to tackle the complexity of motor control. In a modular architecture control is subdivided among basic building blocks, allowing for an efficient yet flexible task decomposition. In particular, a modular generation of the muscle patterns might allow for a low-dimensional representation of the motor output incorporating knowledge on the dynamic behavior of the musculoskeletal system into a small set of basic functions shared across tasks and conditions. Recently, the modular control hypothesis has been supported by observations of low-dimensionality in the muscle patterns underlying a variety of motor behaviors in different species. Using multidimensional decomposition techniques such as principal component analysis (PCA), factor analysis (FA), independent component analysis (ICA), and non-negative matrix factorization (NMF) it has been possible to reconstruct the muscle activation patterns as the combination of a small number of components (Tresch et al., 2006; Giszter et al., 2007; Ting and McKay, 2007; Bizzi et al., 2008; Tresch and Jarc, 2009; Lacquaniti et al., 2012; d'Avella and Lacquaniti, 2013). These components may capture different features of the muscle patterns shared across task conditions, such as specific relationships in the strength of activation of groups of muscles, i.e., muscle synergies (Tresch et al., 1999; Ting and Macpherson, 2005) or M-modes (Krishnamoorthy et al., 2003), specific time-courses of the activation waveforms for all muscles, i.e., temporal components (Ivanenko et al., 2004; Dominici et al., 2011), and specific collections of muscle activation waveforms, i.e., time-varying muscle synergies (d'Avella et al., 2003, 2006) but they all construct muscle patterns by linear combinations of a small number of generators. However, even if muscle patterns can be accurately described by such generators, task accomplishment depends on the actual joint torques and the consequent joint motions produced by muscle contractions. Thus, to better understand how motor tasks may be accomplished by the combination of a few muscle pattern generators it is necessary to assess the relationship between the organization of muscle patterns and that of joint torques.

While joint torques underlying many different motor behaviors have been investigated extensively, a characterization of their dimensionality with multidimensional decomposition approaches such as those recently used to analyze muscle patterns is missing. Joint torque generators for a two-joint arm have been identified before using NMF from simulated data (Chhabra and Jacobs, 2006) but not from experimental data. Focusing on reaching movements in vertical planes, as in many previous studies (Soechting and Lacquaniti, 1981; Lacquaniti et al., 1982; Flanders et al., 1994, 1996; d'Avella et al., 2006, 2008, 2011), our aim was to investigate the dimensionality of joint torques and to compare it with the dimensionality of the muscle patterns. Moreover, we wanted to explore systematically the dimensionality of different types of generators, i.e., generators capturing shared structure in the spatial–across joints or muscles–, temporal, and spatiotemporal dimensions. Planar point-to-point reaching movements, for which joint torques can be estimated using a simplified dynamical model of the arm with two joints, are normally associated with bell-shaped velocity profiles and biphasic torque profiles (Morasso, 1981; Soechting and Lacquaniti, 1981). The shape of such profiles is invariant with respect to movement speed (Soechting and Lacquaniti, 1981) or load (Lacquaniti et al., 1982) and the relationship between shoulder and elbow dynamic torques is almost linear (Soechting and Lacquaniti, 1981; Lacquaniti et al., 1986; Gottlieb et al., 1997). Similar observation were made for reaching movements in three-dimensional space (Lacquaniti et al., 1986). These observations indicate that joint torques for reaching have remarkable regularities suggesting that their dimensionality is also low. One might hypothesize that there is a one-to-one relationship between muscle pattern generators and torque generators. However, biomechanical characteristics and constraints must be taken into account.

To generate the joint torques τ(t) required to move the arm along a given joint trajectory q(t), i.e., torques for which the trajectory is a solution of the arm motion equations [see Equation (3) in Material and Methods], the CNS, according to the modular control hypothesis, combines a set of N_m (spatial, temporal, or spatiotemporal) muscle pattern generators:

$m(t)={\displaystyle \sum _{n=1}^{N_m}{a}_n}{v}_n(t)$

where v_n(t) is the n-th spatiotemporal generator or the product of the n-th temporal component times the n-th spatial weighting vector for spatial and temporal generators (see “Dimensionality of motor commands” in Materials and Methods). The tension generated by the activation of the each muscle is determined by the dynamics of the musculotendon unit, which depends non-linearly on muscle length, velocity, and muscle activation. Muscle length and velocity depend on joint angles and joint velocities via a matrix of moment arms. Then, muscle torque depends on joint angles, joint velocities, and muscle activations

$\tau =\tau (q,\dot{q},m).$

Thus, the required torque profile can be generated by appropriate combination of the muscle pattern generators, i.e., it can be expressed as a function of the combination coefficients a:

$\tau =\tau (q,\dot{q},{\displaystyle \sum _n{a}_n}{v}_n).$

The torque does not depend in general linearly on the muscle activations and, consequently, on the combination coefficients a. When linearity is an adequate approximation, muscle torque can be expressed as a linear combination of the “force fields” associated to each generator, φ_n = φ_n(q, $\dot{q}$) = τ(q, $\dot{q}$, e_n), where e_n is the unit vector along the n-th dimension in coefficient space:

$\tau ={\displaystyle \sum _n{a}_n}{\phi }_n(q,\dot{q}),$

i.e., limb control can be achieved by combination of force-field primitives (Bizzi et al., 1991; Giszter et al., 1993, 2007; Mussa-Ivaldi et al., 1994; Kargo and Giszter, 2000a,b; Kargo et al., 2010; Giszter and Hart, 2013). However, torques profiles observed across different task conditions can also be expressed as (or approximated by) a linear combination of N_τ torque generators:

$\tau (t)={\displaystyle \sum _{n=1}^{N_{\tau }}{b}_n}{u}_n(t).$

Even if there is a one-to-one relationship between muscle pattern generators and force-field primitives, the number of muscle pattern generators must be larger than the number of torque generators because of the non-negativity constraint on muscle activations. As muscles can only pull, muscle pattern generators are combined with non-negative combination coefficients and, even considering a linear muscle-to-torque mapping, to generate torques spanning a N_τ-dimensional space at least N_τ + 1 non-negative generators are required (Davis, 1954; Valero-Cuevas, 2009). Moreover, because of the redundancy of the muscular system, different muscle patterns can generate the same torques and thus the muscle pattern dimensionality can be larger than the minimum imposed by non-negativity, i.e., N_m ≥ N_τ + 1. Thus, the minimum number of muscle pattern generators depends on the actual dimensionality of the joint torque required to perform all conditions of a specific task and the actual number of muscle pattern generators can be larger than the minimum and must be determined experimentally. Importantly, we consider here tasks whose conditions can be described by a set of parameters, such as, for example, the position of a target of a point-to-point reaching movement. Then, since the skeletal system is also redundant for the performance of many tasks, e.g., a specific position of the wrist in space can be achieved with many different joint angle configurations, the actual dimensionality of the joint torques may be lower than the number of joints (i.e., DoF) involved and it must also be determined experimentally. Finally, while the minimal number of muscle pattern generators might guarantee an optimal solution in terms of computational complexity, it might be suboptimal in terms of other costs such as muscular effort. Thus, comparing the torque and muscle pattern dimensionality can provide new information on the control strategy employed to perform a specific task.

We analyzed EMGs data recorded from 19 muscles and kinematic data collected from markers positioned on the arm of subjects performing fast reaching movements from one starting position to 8 targets on the sagittal plane and eight targets on the frontal plane. We used a dynamic model of the arm with four rotational joints (three at the shoulder and one at the elbow) and three translational DoF (the position in space of the shoulder) to estimate joint torques from joint angles with an inverse dynamics computation (Corke, 1996). We then considered the dynamic component of the torques, i.e., the total torques with the gravitational components removed (Gottlieb et al., 1997), and the phasic component of the muscle activity waveforms, i.e., the total rectified and filtered EMG waveforms with the tonic, anti-gravity components removed (Flanders and Herrmann, 1992; d'Avella et al., 2006). Spatial, temporal, and spatiotemporal torque generators were identified by performing PCA on different arrangements of the data matrix. Similarly, spatial, temporal, and spatiotemporal muscle pattern generators were identified with NMF. We first determined the dimensionality of generators according either to a threshold on the fraction of data variation explained (Tresch et al., 1999; Ting and Macpherson, 2005; Torres-Oviedo et al., 2006; Roh et al., 2012) or to the detection of a “knee” in the curve of the variation explained as a function of the number of generators (d'Avella et al., 2003; Cheung et al., 2005; d'Avella et al., 2006; Tresch et al., 2006). We used the former criterion for the torque data and the latter for the EMG data. However, to directly compare the dimensionality of torques and muscle patterns, we then also used a single criterion which took into account the different intrinsic variability of the two datasets when determining their dimensionality (Cheung et al., 2009).

Materials and Methods

Participants, Experimental Apparatus, and Task

Four right handed subjects (aged between 27 and 40) gave their written informed consent to participate in the study, which conformed with the Declaration of Helsinki and had been approved by the Ethical Review Board of the Santa Lucia Foundation. The experimental apparatus and reaching task has been described in details in a previous report (d'Avella et al., 2006). Briefly, standing subjects gripped with their right hand an handle (weight 180 g) which had a sphere (diameter 4 cm) attached to one extremity. The center of sphere was aligned with the axis of the forearm at a distance of 12 cm from center of the palm. Participants were instructed to move the sphere between a central position and 8 targets uniformly arranged on a circle at 15 cm of distance on either the frontal or sagittal plane while minimizing shoulder and wrist movements. The central position was adjusted for each subject so that it required maintaining the upper arm vertical and aligned to the trunk and the elbow flexed at 90°. The targets were indicated by transparent spheres lighted from inside by an LED. In each trial, after holding the sphere at the start position for at least 1 s, subjects were instructed to move after a go signal, to reach the target with a movement of a duration (defined as the interval in which the speed of the sphere was above 10% of its maximum) shorter than 400 ms, and to hold there for at least 1 s. Unsuccessful trials were repeated. Each subject performed each movement successfully five times in different blocks of trials for a total of 160 point-to-point movements (2 planes × 8 targets × 2 directions -from the center to the target and from the target back to the center- × 5 repetitions).

Data Acquisition

The motion of the arm was recorded using an optic motion-tracking system (Optotrack 3020, Nothern Digital, Waterloo, Ontario, Canada) with a sampling frequency of 120 Hz and spatial resolution below 0.1 mm. Active optical markers were positioned on the shoulder (acromion), the upper arm (at the proximal end close to the head of the humerus), the elbow (epicondylus lateralis), the wrist (one over the styloid process the radius and one on the styloid process of the ulna). The motion of the sphere on the handle (end-point) was recorded with an electromagnetic motion-tracking system (Fastrak, Polhemus, Calchester, VT) with sampling frequency of 120 Hz and spatial resolution below 4 mm, as estimated by a calibration process performed within the workspace used in the experiment.

EMG activity was recorded with active bipolar surface electrodes (DE 2.1; Delsys, Boston,MA) from the following muscles: biceps brachii, short head (BicShort), biceps brachii, long head (BicLong), brachialis (Brac), pronator teres (PronTer), brachioradialis (BrRad), triceps brachii, lateral head (TrLat), triceps brachii, long head (TrLong), triceps brachii, medial head (TrMed) deltoid, anterior (DeltA), deltoid, middle (DeltM), deltoid, posterior (DeltP), pectoralis major, clavicular portion (PectClav), pectoralis major, lower portion (PectLow), trapezius superior (TrapSup), trapezius middle (TrapMid), trapezius inferior (TrapInf), latissimus dorsi (LatDors), teres major (TeresMaj), infraspinatus (InfraSp). EMG signal was band-pass filtered (20–450 Hz) and amplified (total gain 1000, Bagnoli-16, Delsys Inc.). EMG data were digitized at 1 KHz (PCI-6035E, National Instruments, Austin, TX).

Data acquisition and experiment control were performed on a workstation with custom software written in LabView (National Instruments, Austin, TX). Fastrak data were processed on-line to compute the movement time and target accuracy and to provide auditory feedback about unsuccessful trials. The experiment control program logged the time of all relevant behavioral events.

Data Analysis

End point kinematics

All analyses were performed with custom software written in Matlab (Mathworks, Natick, MA). Position and orientation of the handle and the measured geometric parameters of the handle were used to compute the position of the end-point. The data were low-pass filtered (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 computed the onset time and the end time, defined respectively as the time in which the speed profile crossed 10% of its maximum value, and the movement duration (MT), defined as the interval between the movement onset and the movement end.

Arm model

A kinematic and kinetic model of the arm, incorporating geometrical and inertial parameters of the upper arm and forearm segments, was used to estimate joint angles and joint torques from the recorded spatial position of the shoulder, the elbow, and the wrist markers. The kinematic model was developed using the Denavit-Hartenberg (D-H) notation (Hartenberg and Denavit, 1955), i.e., as chain of articulated links with four parameters for each link (a: length, α: twist, d: offset, ϑ: joint angle) describing the position and orientation of a Cartesian reference frame fixed on each link with respect to the reference frame fixed on the preceding link of the chain according to the 4 × 4 homogeneous transformation matrix T:

$\begin{array}{cc}T=\left[\begin{array}{cccc}\cos \vartheta & -\sin \vartheta \cos \alpha & \sin \vartheta \sin \alpha & a\cos \vartheta \\ \sin \vartheta & \cos \theta \cos \alpha & -\cos \theta \sin \alpha & a\sin \theta \\ 0& \sin \alpha & \cos \alpha & d\\ 0& 0& 0& 1\end{array}\right].& \left(1\right)\end{array}$

The rotation axis of each joint coincides with the z axis of the preceding link in the chain. The x axis in each frame is directed as the normal between the z axis of that frame and the z axis of the next frame. In this way the joint angle is the angle between the x axes of the frames of the two links connected by the joint. We modeled four rotational degrees-of-freedom (DOFs) of the arm—three rotations at the shoulder, i.e., adduction, flexion and external rotation, and one rotation at the elbow, i.e., elbow flexion (see Figure 1)—and three translational DOFs of the shoulder. We assumed that shoulder was a spherical joint (i.e., the rotation axes of the three joints intersect at a single point). Lengths of upper arm, forearm, and hand of each subject were estimated as a function of the subject's weight and height according to regression equations (Winter, 1990). Forearm, hand, and handle were considered a single link (7th) of length equal to the sum of the forearm length and the length of the opened hand, thus approximating the total length of the closed hand and the handle along the direction of the forearm axis with the length of the opened hand.

FIGURE 1

Figure 1. Joint angle definition for the arm model. The four joint angles included in the model (shoulder adduction, shoulder flexion, shoulder external rotation, and elbow flexion) are illustrated by a sequence of postures in space of a two-link arm. The four rotation axes are indicated by colored arrows and correspond to the z-axes of the four reference frames defined according to the D-H notation (see text) and labeled z₀–z₃.

The kinetic model of the arm was developed adding to each link its inertial parameters (mass, center of mass, inertia tensor) also estimated as a function of the subject's weight and height according to regression equations (Zatsiorsky and Seluyanov, 1983). No mass was associated to the first three links required to represent the spatial position of the shoulder. However, these translational DOFs were introduced to take into account shoulder movements when estimating the joint torques. The mass of the upper arm was assigned to the 6th link, which had an offset equal to the length of the upper arm segment. The mass of the forearm, hand, and handle was assigned to the 7th link, associated with the elbow flexion. The inertial parameters for this link were computed from the inertial parameters estimated separately from the regression equations for the forearm and hand. As the estimated position of the center of mass of the hand and of the handle coincided, the mass of the handle (180 g) was summed to the mass of the hand. The moments of inertia were computed with respect to its center of mass. The model was implemented in Matlab using the Robotic Toolbox (Corke, 1996, 2011). The D-H parameters of the generic arm model are reported in Table 1 and the specific geometric and inertial parameters estimated for each subject are reported in Table 2.

TABLE 1

Table 1. D-H parameters for the 7 DOFs arm model.

TABLE 2

Table 2. Arm model parameters for individual subjects.

Joints kinematics

The arm model was used to estimate at each time sample the shoulder adduction angle, the shoulder flexion angle, the shoulder external rotation angle, the elbow flexion angle using the positions of the shoulder and elbow markers and the mean position between the two wrist markers. For each time sample and each joint angle, a vector between two markers aligned with the axis of the limb segment defining the rotation of that joint (i.e., shoulder and elbow markers for shoulder adduction and shoulder flexion, elbow and wrist markers for shoulder external rotation and elbow flexion) was computed first. Then, the segment vector was transformed into the reference frame associated to the joint according to the matrices defined by Equation (1) and the angle computed as

$\begin{array}{cc}{\vartheta }_i={\tan }^{-1}(y/x)& \left(2\right)\end{array}$

where x and y are the coordinate of the vector in the reference frame associated with the joint rotation axis (z axis). To compensate for potential misalignment between the tracker z axis and the vertical axis, the coordinates of the markers were first rotated into a Cartesian reference frame with the gravitational acceleration along the z axis. The direction of the gravitational acceleration was estimated by means of a calibration of based on tracking two markers attached to the fulcrum and the extremity of a pendulum.

Angular velocity and acceleration were computed by numerical differentiation. To validate the kinematic model, forward kinematics was used to compare estimated and measured end-point trajectories.

Inverse dynamics

Joint angles, joint velocities and joint accelerations were used to estimate the torque profiles via recursive Newton-Euler calculation (rne function of Matlab Robotics Toolbox). We computed the total torques τ

$\begin{array}{cc}\tau =M(q)\ddot{q}+C(q,\dot{q})\dot{q}+G(q)& \left(3\right)\end{array}$

where M is the manipulator inertia matrix, C is the Coriolis and centripetal torque, and G is the gravitational torques. To estimate non-gravitational (dynamic) torques we subtracted gravitational torques from the total torques.

To validate the inverse dynamics calculation we also performed a forward dynamics simulation (fdyn function of Matlab Robotics Toolbox) using the arm model and the estimated torque profiles to reproduce the original joint angle trajectories.

Data preprocessing

The EMGs for each trial were digitally full-wave rectified, low-pass filtered (FIR filter, 20 Hz cut-off, zero-phase distortion, Matlab fir1, and filtfilt functions), and integrated over 10 ms intervals. In a few cases muscle waveforms showed some artifacts, possibly due to a partial detachment of the electrode from the skin, or to an high correlation between two or more muscles, and those muscles were removed from further analysis (subject 1: PectLow; subject 2: TrapInf, PectLow).

EMGs and torques for all the trials in each experimental condition (2 planes × 8 targets × 2 directions) were aligned on the time of movement onset and averaged.

Finally, both torques and muscle waveforms were normalized in time to equal MT and resampled with 50 samples per MT. Samples from 0.5 MT before movement onset to 0.5 MT after movement end (total 100 samples) were considered for further analysis.

Dimensionality of motor commands

We consider a set of D command signals (joint torques or muscle patterns) delivered by a controller in a given time interval (sampled T times) to accomplish a task in one of K distinct task conditions (e.g., different reaching targets). We hypothesize that a modular controller generates these command signals by modulating and combining a small set of generators whose structure is invariant across all task conditions. The structure of such generators may be defined in the spatial (across signals, i.e., muscles or joints), temporal, and spatiotemporal domains. The dimensionality of the ensemble of command signals is then simply the number of generators necessary to accomplish all K tasks conditions.

Spatial dimensionality is the number of generators necessary to capture time-invariant relationships between the signals. For N generators:

$\begin{array}{cc}{x}^k(t)={\displaystyle \sum _{n=1}^N{c}_n^k}(t)w_n& \left(4\right)\end{array}$

where x^k(t) are the set of signals for condition k, i.e., a vector-valued (D-dimensional) function time (or a D × T matrix for discrete time samples), c^k_n(t) is a condition-dependent, time-varying combination coefficient for the n-th generator, w_n is the condition-independent, time-invariant n-th spatial generator, i.e., a D-dimensional vector capturing the relative activation weight of different signals.

Temporal dimensionality is the number of generators necessary to capture temporal components shared across all signals (i.e., space-invariant). For N generators:

$\begin{array}{cc}{x}^k(t)={\displaystyle \sum _{n=1}^N{c}_n}(t)w_n^k& \left(5\right)\end{array}$

where x^k(t) are again the set of signals for condition k, c_n(t) is the condition-independent, time-varying n-th generator (or temporal component), w^k_n is the condition-dependent, time-invariant n-th D-dimensional weight vector for the n-th component. Notice how the critical difference in the definition of spatial and temporal generators and dimensionality is in the dependence on the task condition (k). Indeed, generators are useful concepts only if they can be used for a variety of conditions thus allowing an effective reduction of the number of parameters to select for each condition.

Spatiotemporal dimensionality is the number of generators capturing simultaneously invariant spatial and temporal features in the signals. Thus, each generator includes a set of signal components that can be expressed as a time-varying vector. For N generators

$\begin{array}{cc}{x}^k(t)={\displaystyle \sum _{n=1}^N{a}_n^k}{v}_n(t)& \left(6\right)\end{array}$

where a^k is a condition-dependent combination coefficient for the n-th generator, v_n(t) is the n-th spatiotemporal generator, i.e., a condition-independent, time-varying D-dimensional vector (or a D × T matrix for discrete time samples). However, as different signals may be related synchronously or asynchronously, we can distinguish the case of synchronous spatiotemporal generators:

$\begin{array}{cc}{v}_n(t)=c_n(t)w_n& \left(7\right)\end{array}$

in which each generator v_n(t) can be expressed as the product of a scalar function of time c_n(t)times a time-invariant weight vector w_n. In contrast, asynchronous spatiotemporal generators cannot in general be factorized into separate spatial and temporal generators.

In addition to being scaled in amplitude, spatiotemporal generators may also be recruited at different times across task conditions, i.e., they may also show invariance for time shifts (d'Avella et al., 2003, 2006). If we assume that the duration of each spatiotemporal generator is smaller than the duration of the signals, we can incorporate condition-dependent onset times t^k_n into Equation (6):

$\begin{array}{cc}{x}^k(t)={\displaystyle \sum _{n=1}^N{a}_n^k}{v}_n(t-t_n^k)& \left(8\right)\end{array}$

Identification of generators and their dimensionality

To investigate the spatial, temporal, and spatiotemporal dimensionality of joint torques and muscle patterns we used multidimensional decomposition techniques to identify the different types of generators. We considered the dynamic component of the torques and the phasic component of the muscle activity waveforms. We then used PCA to identify torque generators and, because of the inherent non-negativity of muscle activity, we used NMF to identify muscle pattern generators. Finally, as discussed below, we selected the number of generators with three different criteria, two specific for each dataset and one for both datasets.

Dynamic torques and phasic muscle patterns. Reaching movements in vertical planes require torques and muscle activities to accelerate and decelerate the limb and to balance gravitational forces. In this work we focused on the former components, i.e., dynamic torques and phasic muscle patterns. Dynamic torques were computed as the total torques with the gravitational components [the last term of the right hand side of Equation (3)] removed (Gottlieb et al., 1997). Flanders and collaborators (Flanders and Herrmann, 1992) found that it is possible to distinguish the phasic component (related to the movement) from the tonic one (related to maintain a specific posture of the arm) of an EMG signal. As in d'Avella et al. (2006) we used a subtraction procedure to remove the tonic component, i.e., we subtracted a constant muscle activation level before and after the movement and a linear ramp between the two constant values during the movement. After the subtraction a small fraction of EMG samples assumed negative values, indicating that the phasic EMG activity was lower than the tonic activity. However, in order to use the NMF algorithm, we set to zero all negative values (ratio of negative area over total area of all muscles, 0.15 ± 0.04, mean ±SD over subjects). To assess the effect of this procedure on the number of generators and on their structure, we also identified generators from the original phasic muscle patterns without imposing a non-negativity constraint and using an iterative factorization algorithm based on gradient descent in place of NMF (see below).

Data matrices. To identify spatial, temporal, and spatiotemporal generators, joint torque and muscle patterns data, after pre-processing, were organized into three different data matrices that were factorized by either PCA (torques) or NMF (muscle patterns). For each subjects we identified generators from K task conditions (K = 32, except for subject 3 for which we had to exclude 2 conditions on the frontal plane and 2 on the sagittal plane because of missing data from the arm markers used to compute joint angles). To identify spatial generators, the data for each condition (D signals, EMG or torque, times T samples, with T = 100 after time normalization and resampling, see Figure 2) were arranged into a data matrix X with D row and T × K columns which was factorized, according to Equation (4) in matrix notation, as X = W C, where W is the condition-independent synergy matrix with D rows and N columns, N number of generators, and C is the matrix of condition- and time-dependent combination coefficients with N rows and T × K columns. For temporal generators, in contrast, the data matrix was constructed by arranging the waveforms from all signals in all conditions as columns, i.e., X had T rows and D × K columns, and it was factorized, according to Equation 5 in matrix notation, as X = C W, with C is the condition-independent matrix of temporal components, with T rows and N columns, and W is the condition- and signal-dependent matrix of weights, with N rows and D × K columns. Finally, for spatiotemporal generators, the data samples for all signals of each conditions were arranged in a column and the data matrix X, with D × T rows and K columns, was factorized, according to Equation 6 in matrix form, as X = V A, with V condition-independent matrix of time-varying synergies with D × T rows and N columns and A condition-dependent matrix of combination coefficients with N rows and K columns. For joint torque generators, the covariance of the data matrix was computed and, for each N, the first N principal components (extracted using MATLAB pcacov function) were considered. For muscle pattern generators, for each N, C and W matrices were initialized randomly and the best solution out of 20 runs of the NMF algorithm was selected. Each run of the iterative algorithm was terminated when the reconstruction R² increased in one iteration by less than 10⁻⁴ for five consecutive iterations. To assess the effect of clipping to zero the negative values of the phasic muscle patterns, we also identified generators without non-negativity constraints using an iterative gradient descent algorithm. We minimized the data reconstruction error by combination of spatial, temporal, and spatiotemporal generators iterating a step in which combination coefficients (C and A), in the spatial and spatiotemporal cases, and weights (W), in the temporal case, were updated and a step in which synergies (W and V), in the spatial and spatiotemporal cases, and temporal components (C), in the temporal case, were updated. Both updates were performed along the direction opposite to the error gradient with step size of 0.05 for the combination coefficients (spatial and spatiotemporal cases), 0.0018 for the weights (temporal case), 0.0005 for the synergies (spatial and spatiotemporal cases), 0.0019 for the temporal components. These gradient step sizes were selected within a range of values in order to achieve the highest reconstruction R². As in d'Avella et al. (2006), in the spatial and spatiotemporal cases, we added a term to the error function to penalize large negative values in the identified synergies. The same number of runs and termination conditions as for the NMF algorithm were used.

FIGURE 2

Figure 2. Types of multidimensional decomposition of data from different task conditions. Data collected from D channels (4 in this schematic illustration, represented by different patterns) over T time samples (6, represented by different color saturations) in a single task condition are represented by a grid of squares. Different task conditions are represented with different background colors. A spatial decomposition is obtained by factorizing the data matrix obtained by stacking the data from individual conditions horizontally (i.e., matching their spatial—channels—dimension) into a matrix of N (3) spatial generators (D rows and N columns) times a matrix of time-and condition-dependent coefficients. A temporal decomposition is obtained by factorizing the transpose of data matrix obtained by stacking the data from individual conditions vertically (i.e., matching their temporal dimension) into a matrix of N (3) temporal generators (T rows and N columns) times a matrix of channel- and condition-dependent coefficients. Finally, a spatiotemporal decomposition is obtained by arranging all the data samples of each condition into a column and factorizing the resulting matrix into a matrix of N (3) spatiotemporal generators (D × T rows and N columns) times a matrix of condition-dependent coefficients.

Selection of the number of generators. For torque generators, we selected their number as the minimum number which explained at least 90% of the data variation (VAF or R², defined as 1—SSE/SST, with SSE sum of square residuals of the data reconstruction by the generators, and SST sum of the squared residuals of the data with respect to the mean over the rows of the data matrix). Such criterion (“R² threshold”) has been frequently used in the muscle synergy literature (Tresch et al., 1999; Ting and Macpherson, 2005; Torres-Oviedo et al., 2006; Roh et al., 2012), even if sometimes with a different definition (i.e., with SST defined as the sum of the squared data, see Delis et al., 2013). Such criterion is based on the assumption that the fraction of data variation unexplained is due to noise and the threshold is supposed to separate structured variation due to the combination of generators and noise. However, if an independent estimation of the noise level is not available the choice of such threshold is necessarily ad-hoc. An alternative approach, also used in previous studies (d'Avella et al., 2003; Cheung et al., 2005; Tresch et al., 2006), that we used for selecting the number of muscle pattern generators is the detection of a “knee” in the curve of R² as a function of the number of generators. Such criterion (“R² knee”) relies on the assumption that the noise is isotropic, i.e., contributes equally to all dimensions, and does not depend on a specific assumption of the relative level of noise. To detect a change in slope in the R² curve, for each N, we performed a linear fit of the portion of the curve from N to the end (i.e., D) and we selected N for which the mean square error of the fit was <10⁻⁴, indicating that the “tail” of the curve after the “knee” was essentially straight. We could not use this second criterion for the torques as their maximum spatial dimension (4) was too low and it was impossible to identify a “knee” with such procedure. However, to compare torque and muscle pattern with the same criterion, we also determined their dimensionality with a criterion (“R² shuffle”) that took into account the different intrinsic noise levels of the two datasets. We then used a threshold on the slope of the R² curve according to slope of the curve obtained after a random shuffling the rows of the data matrix (Cheung et al., 2009). Data samples after shuffling were low-pass filtered to match the smoothness in the original data. By shuffling the data the multidimensional structure of the original data was lost but each dimension maintained the original variability. Thus, we selected the number of generators as the point on the original R² curve at which any further increase in the number of extracted generators yielded an R² increase smaller than 75% of that for the generators extracted from the shuffled data (mean over 50 extractions from reshuffled data).

Comparison of generators across subjects. To compare generators across subjects, we tested how well a set of generators identified in one subject could reconstruct the data of a different subject. We then computed a R²-value to assess the similarity of the subspaces spanned by the generators of different subjects. To assess the significance of these R²-values we performed a Monte Carlo simulation identifying generators from random data obtained by randomly shuffling the original data (50 runs for each subject and type of generator). We then computed the 95% percentile of the distribution of R²-values for the reconstruction of the original data with generators identified from random data.

Effect of potential contamination of EMG recordings by cross-talk. As surface EMG recordings can be affected by cross-talk due to volume conduction of the EMG signal from neighboring muscles, we performed a Monte Carlo simulation to assess the effect of such potential contamination on the dimensionality of the spatial muscle pattern generators. For each muscle, i.e., the i-th row of the normalized data matrix X, we simulated a cross-talk contamination from a second muscle (j-th row), randomly chosen from all other muscles and not limited to the neighboring ones, according to a cross-talk weight (α) randomly drawn from an exponential distribution with mean 0.1, i.e., X′_ik = X_ik + α X_jk. We then identified the generators from the contaminated data matrix X′ with NMF and we estimated their dimensionality. For each subject, we performed a total of 100 simulation runs, and, across all subjects, we found that the dimensionality estimated using the contaminated data matrix was different from the dimensionality estimated using the original matrix only in 6.7% of runs (0 runs for subjects 1, 2, 4; 27 runs for subject 3, mean dimensionality difference 0.27). Thus, while we cannot exclude that our EMG recordings were not affected by cross-talk, we are confident that such potential contamination did not significantly alter muscle pattern dimensionality estimation.

Results

Dynamic Torques

Joint torques were estimated by inverse dynamics from joint angle trajectories using a kinetic model of the arm parametrized by the height and weight of each subject. To validate the arm model and the results of the inverse dynamics computation we performed a forward dynamic simulation. The joint angle trajectories simulated using the torques estimated by inverse dynamics matched well the original joint angle trajectories of all subjects and conditions (R² = 0.99 ± 0.02, mean across subjects ±SD). Figure 3 shows an example of end point trajectories, end point speed profiles, joint angle trajectories, angular velocities, and gravitational and dynamic joint torque profiles for eight center-out movements on the frontal plane. As expected, end point trajectories were straight and velocity profiles bell-shaped. Joint angle trajectories and the corresponding angular velocities were modulated by movement direction. Dynamic torques, i.e., total torque with the gravitational torques removed, were bi-phasic, as observed before (Gottlieb et al., 1997). The time courses of the joint angle trajectories and angular velocities were different across joints and conditions but, because of the dynamic interaction between the different DOFs and they were generated by a synchronous biphasic pulse of torque distributed across joints with different balances depending on the movement direction. Such coordination patterns in the dynamic torque profiles is clearly visible in a scatter plot of a pair of joint torques. Figure 4 shows the six scatter plots of all pairs of joint torque profiles, during an interval of 250 ms around movement onset, approximately capturing the first phase of the profile, for the same 8 movements of Figure 3. If a pairs of torques were modulated synchronously, the corresponding trajectory in the scatter plot would appear as a straight line segment with a direction depending on the relative amplitude. Indeed, for most pairs and movement directions dynamic torques appeared to be modulated close to synchronously, especially in the initial (raising) portion of the profile. Finally, for two pairs of dynamic torques, shoulder external rotation-shoulder adduction and elbow flexion-shoulder flexion, the direction of the line segment in the scatter plot depended only weakly on the movement direction, suggesting that the dynamic torques were spanning a subspace of the four dimensional torque space orthogonal to those two directions. We then generalized these observations by identifying dynamic torque generators and estimating their dimensionality.

FIGURE 3

Figure 3. Example of endpoint speed, velocity, joint angles and torques. Example of endpoint trajectories, end-point speed profiles, joint angles, joint angular velocities, gravitational (light gray) and dynamic (dark gray) torques for eight center-out movements in the frontal plane of subject 1 (mean across repetitions of each movement). Vertical dashed lines represent the times of movement onset and movement end. Shaded areas around gravitational and dynamic torque profiles represent ±1 SD around the mean.

FIGURE 4

Figure 4. Example of coordination between pairs of dynamic torques. Each scatter plot illustrated the dynamic torques for a pair of joints recorded in an interval of 250 ms around the time of movement onset for 8 center-out movements in the frontal plane of subject 1.

Spatial dimensionality

We first assessed the spatial dimensionality of the dynamic torques by identifying spatial generators, i.e., vectors in the torque space capturing specific balances of torque magnitude which could reconstruct the data once multiplied by time- and condition-dependent coefficients (see Materials and Methods and Figure 2), using PCA. For each subject, the number of generators was selected as the minimum number for which the fraction of data variation explained exceeded 0.9 (“R² threshold” criterion) and as the number of generators for which adding an additional generator increased the R²-value less than 75% of the mean R²-values obtained identifying generators from shuffled data (“R² shuffle” criterion). The mean dimensionality across subjects was 2.25 according to the R² threshold criterion and 2.75 according to the R² shuffle criterion (see Table 3 for individual values). The maximum potential spatial dimensionality of the torques was 4, corresponding to the number of joints, i.e., the number of rows of the data matrix used for spatial decomposition (Figure 2).

TABLE 3

Table 3. Comparison of different types of dimensionality of dynamic torques and phasic muscle patterns estimated according to three criteria for the selection of the number of generators.

Figure 5A shows the R²-value as a function of the number of generators for subject 1 and Figure 5B the three spatial generators of the same subject selected according to the R² shuffle criterion. The first generator (w₁) is dominated by shoulder flexion torque. The second generator (w₂) combines a large shoulder adduction torque with a smaller shoulder internal rotation (i.e., negative external rotation) and elbow extension (i.e., negative elbow flexion). Finally, the third generator (w₃) represents a large elbow flexion torque and a smaller shoulder adduction torque. Notably, none of the generators or their combinations can generate coordinated shoulder adduction and shoulder external rotation torques, i.e., the direction orthogonal to the torques direction observed in the corresponding scatter plot of Figure 4. Thus, the structure of the spatial generators indicated that such torque coordination pattern was never used to perform reaching movements in the frontal and sagittal planes.

FIGURE 5

Figure 5. Spatial decomposition of dynamic torques. (A) R² curve for subject 1 obtained by spatial decomposition using PCA. (B) Three spatial generators selected for subject 1. (C) Example of the reconstructions of the dynamic torques for six movement conditions of subject 1 obtained with the generators illustrated in panel (B) (shaded area: original data, thick line: reconstructed data, bottom: time-varying combination coefficients).

Figure 5C illustrates an example of the reconstruction of the dynamic torque profiles of subject 1 in six different conditions by the combination of the three spatial generators of Figure 5B. The dynamic torques for the first two conditions, medial and lateral movements in the frontal plane, are generated by a comparable level of activation of all three generators with a bi-phasic activation of shoulder adduction and internal rotation followed by shoulder abduction and external rotation for the medial movement and the opposite order for the lateral movement captured mainly by the activation of the second generation with similar bi-phasic profiles but opposite signs of its combination coefficient (c₂). The last two conditions, backward and forward movements in the sagittal plane, require large shoulder flexion/extension torques that are generated by a bi-phasic activation of the first generator, captured by the first time-varying coefficient (c₁).

To assess the similarity between the subspaces spanned by the generators identified in each subject we reconstructed all dynamic torques of each subject with the generators of all subjects. Table 4 shows the R²-values obtained using the number of generators determined according to the R² shuffle criterion (see Table 3). The R²-values for the reconstruction of the data of each subject by the generators extracted from the other subjects (0.96 ± 0.04, mean ±SD, n = 12) were close to the R²-values of the reconstruction of the data of each subject by the generators extracted from the same data (0.98 ± 0.02, n = 4) and were significantly higher than the R²-values obtained with generators identified from randomly shuffled data, indicating that the dynamic torques of the different subjects shared a similar spatial organization.

TABLE 4

Table 4. R² for the reconstruction of the data of each subject with the torque generators identified in all subjects.

Temporal dimensionality

To identify generators of the temporal organization of dynamic torques we performed PCA on the collection of the torque profiles of all joints and conditions. The resulting temporal components were then waveforms with the same duration as the torque profiles and each profile was reconstructed by multiplying the component matrix by a weight specific for that joint and condition. The dimensionality was 1 for all subjects and for both criteria (see Table 3). In contrast, the maximum potential temporal dimensionality of the torques was 100, corresponding to the number of time samples after time-normalization and resampling from −0.5 MT before movement onset and 0.5 MT after movement end, i.e., the number of rows of the data matrix used for temporal decomposition (Figure 2).

Figure 6A illustrates the R² curve for the temporal decomposition up to 12 generators for subject 1 and Figure 6B the single temporal component identified in this subject and representative of all subjects, clearly showing a bi-phasic profile. Figure 6C illustrates the reconstruction of the joint torques for the same six conditions of Figure 5C by the temporal generator. The torque profiles for each condition are reconstructed multiplying the single temporal component (c₁) by a single condition-dependent weight vector (w₁). With respect to the reconstruction with spatial generators, the weight vector, which has the same dimensions of a spatial generator, is now modulated by the movement. For example the opposite signs in the bi-phasic profiles of shoulder adduction and external rotation for medial and lateral movements and for shoulder flexion for backward and forward movements are obtained by opposite signs of the components of the weight vector.

FIGURE 6

Figure 6. Temporal decomposition of dynamic torques. (A) R² curve for subject 1 obtained by temporal decomposition using PCA. (B) The single temporal generators selected for subject 1. (C) Example of the reconstructions of the dynamic torques for six movement conditions of subject 1 obtained with the generator illustrated in panel (B) (shaded area: original data, thick line: reconstructed data, bottom: joint-specific weights).

Finally, the temporal generators were also similar across all subjects. As for spatial generators, the reconstruction of the data of each subject by the generators of all other subjects had R²-values (0.93 ± 0.01, mean ±SD) which were comparable with the R²-values for the reconstruction of the data of each subject by the generator extracted from the same data (0.94 ± 0.01) and significantly higher than the R²-values obtained with generators identified from randomly shuffled data.

Spatiotemporal dimensionality

Spatiotemporal generators, which can be viewed as either time-varying vectors capturing a different spatial coordination among torques at each time or as collections of different waveforms for each torque, were identified by PCA on a data matrix obtained arranging all time samples from all joints in a single column for each movement condition. Thus, torque samples from different joints and times represented different dimensions and the possibility of generating the data with a number of generators smaller than the maximum potential dimension (400, corresponding to the number of joints times the number of samples) revealed a coordination in the activation of different joints at different times. Once a set of spatiotemporal generators are identified, the data are reconstructed by multiplying each generator by a single condition-dependent coefficient (see Figure 2). Thus the spatiotemporal decomposition provides a potentially very compact representation of the structure inherent in the data. The mean spatiotemporal dimensionality across subjects was 2.75 according to both criteria (see Table 3 for individual values). Notably, mean spatial and spatiotemporal dimensionalities were very close and even equal for each subject when considering the R² shuffle criterion. Moreover, as the temporal dimensionality was 1, the spatiotemporal dimensionality was essentially the product of the spatial and the temporal dimensionalities.

pt Figure 7A illustrates the R² curve for the spatiotemporal decomposition up to 12 generators and Figure 7B the three spatiotemporal component for subject 1. Comparing the structure of these generators with that of the spatial (Figure 5B) and temporal (Figure 6B) generators of the same subject, it is apparent how each spatiotemporal generator appears as the product of a spatial generator by a temporal one. Indeed, the activation waveforms of all spatiotemporal generators are approximately synchronous and similar to the waveform of the single temporal generator. Figure 7C illustrates the reconstruction of the joint torques for the same six conditions of Figures 5C, 6C. The torque profiles for each condition are reconstructed multiplying each spatiotemporal component by a single condition-dependent coefficient (c_i), represented by the height of the rectangle below the torque profiles. With respect to the reconstruction with spatial and temporal generators, movements requiring torque profiles with opposite signs are generated simply by changing the sign in a single combination coefficient, e.g., c₂ for medial and lateral movements and c₁ for backward and forward movements.

FIGURE 7

Figure 7. Spatiotemporal decomposition of dynamic torques. (A) R² curve for subject 1 obtained by spatiotemporal decomposition using PCA. (B) Three spatiotemporal generators selected for subject 1. (C) Example of the reconstructions of the dynamic torques for six movement conditions of subject 1 obtained with the generator illustrated in panel (B) (shaded area: original data, thick line: reconstructed data, bottom: combination coefficients represented by the height of the rectangle containing the temporal profile of each generators averaged over joints).

As in previous cases, spatiotemporal generators were similar across subjects. The R²-values for the reconstruction of the data of each subject by the generators of all other subjects (0.90 ± 0.04, mean ±SD) were close to the R²-values for the reconstruction of the data of each subject by the generator extracted from the same data (0.95 ± 0.02) and significantly higher than the R²-values obtained with generators identified from randomly shuffled data.

Potential misestimation of subject mass and height

The estimation of the joint torques from the recorded joint kinematics through the inverse dynamic calculation depends on geometric and inertial parameters which are estimated as a function of the subject mass and height according to anthropometric tables (see Arm Model section in Materials and Methods). We assessed the effect of a potential misestimation of such parameters on the estimated torque dimensionality by identifying joint torque generators after varying the mass and the height of each subject by ±5, ±10, ±15, ±20%. We recomputed joint torques for individual trials of all subjects and re-processed the torque data as with the original parameters. Across all subjects and types of generators, the dimensionality was affected by a change in mass in 8 out of 96 cases (4 subjects × 3 types of generators × 8 mass change levels) and by a change in height in nine cases. Thus, the estimation of torque dimensionality is robust to small errors in the estimation of anthropometric parameters.

Muscle Patterns

Phasic muscle patterns, obtained by subtracting the anti-gravity (tonic) components from the rectified, filtered, averaged, time-normalized, and resampled EMG waveforms, were decomposed with NMF to assess their dimensionality. Phasic muscle patterns for fast reaching movements in vertical planes have been described before (d'Avella et al., 2006). In contrast to our previous study, here we identified spatial generators, temporal generators, and spatiotemporal generators without onset delays and we compared their dimensionality with the dimensionality of the corresponding generators of dynamic torques.

Spatial dimensionality

The mean spatial dimensionality of the phasic muscle patterns across subjects was five according to the position of change in slope of the R² curve as a function of the number of generators (R² knee criterion) and five according to the R² shuffle criterion (see Table 3 for individual values). Thus, as expected, the dimensionality of the muscle pattern generator was larger than the number of spatial torque generators (2.75 according R² shuffle criterion) as muscle pattern generators could only be combined with non-negative time- and condition-dependent combination coefficients. However, the number of muscle pattern generators was larger than the minimum required for generating a space of the same number of linear dimensions as the torque generators (2.75) by non-negative combinations (3.75 = 2.75 + 1).

Figure 8A shows the R² curve for the spatial decomposition of the phasic muscle patterns of subject 1, in which a knee at four generators is clearly visible. The lower R²-value at the selected number of muscle patterns generators (0.80) with respect to the corresponding value for the torque generators (0.99) indicated that a much larger fraction of the muscle data variation was due to noise. The four spatial generators (or time-invariant muscle synergies) for the same subject illustrated in Figure 8B (w₁–w₄) show specific groupings of muscles spanning multiple joints and with the same muscle recruited by multiple generators. Finally, in Figure 8C the examples of the reconstruction of the phasic muscle patterns for six movement conditions by the combination of the spatial generators are presented. The temporal structure of muscle patterns and of combination coefficients is clearly more complex than that of the spatial generators with the tri-phasic organization of the muscle patterns generated by both the temporal structure of the combination coefficients and by the superposition of different generators.

FIGURE 8

Figure 8. Spatial decomposition of phasic muscle patterns. (A) R² curve for subject 1 obtained by spatial decomposition using NMF. (B) Four spatial generators selected for subject 1. (C) Example of the reconstructions of the muscle patterns for six movement conditions of subject 1 obtained with the generators illustrated in panel (B) (shaded area: original data, thick line: reconstructed data, bottom: time-varying combination coefficients).

Finally the spatial generators for the muscle patterns were less similar across subjects than the spatial generators for the torques (see Table 5). The reconstruction of the data of each subject by the generators of all other subjects had a R²-values (0.66 ± 0.06, mean ±SD) much lower than the mean R² for the reconstruction by the generators extracted from the same data (0.80 ± 0.04) but still significantly higher than the R²-values obtained with generators identified from randomly shuffled data.

TABLE 5

Table 5. R² for the reconstruction of the data of each subject with the muscle pattern generators identified in all subjects.

Temporal dimensionality

The mean number of temporal generators of the phasic muscle patterns was 4.5 according to the R² knee criterion and 3.75 according to the R² shuffle criterion (see Table 3 for individual values). As for the spatial generators, the temporal dimensionality of the muscle pattern generators was larger than the minimum number required to generate a space with the same linear dimensions as the number of temporal torque generators (1) by non-negative combinations (2 = 1 + 1).

Figure 9A shows the R² curve for the temporal decomposition of the phasic muscle patterns of subject 1 and Figure 9B the four temporal generators (or components) selected in that subject according to both criteria. The first three generators capture a single burst of muscle activity and the fourth component a small burst followed by a larger burst. The four components peak at different times and thus they appear to capture four distinct phases of the muscle patterns observed in different directions. However, the examples of muscle pattern reconstructions and combination weights for six movement directions (Figure 9C) show that in many cases the weight vectors loading the different components were similar within each movement condition (e.g., for the first two components of the medial movement and the last two components of downward movement), suggesting that such temporal decomposition was necessary to capture not only the major changes in the muscle patterns over the duration of the movement but also small asynchronous adjustments.

FIGURE 9

Figure 9. Temporal decomposition of phasic muscle patterns. (A) R² curve for subject 1 obtained by temporal decomposition using NMF. (B) The four temporal generators selected for subject 1. (C) Example of the reconstructions of muscle patterns for six movement conditions of subject 1 obtained with the generator illustrated in panel (B) (shaded area: original data, thick line: reconstructed data, bottom: muscle-specific weights of each generator).

In contrast to the spatial generators but similarly to the temporal generators for torques, muscle pattern temporal generators were similar across all subjects. The R²-values for the reconstruction of the data of each subject by generators of all other subjects (0.82 ± 0.04, mean ±SD) were close to the values for the reconstruction by the generators extracted from the same data (0.85 ± 0.03) and higher than the R²-values obtained with generators identified from randomly shuffled data.

Spatiotemporal dimensionality

The mean number of spatiotemporal generators of the phasic muscle patterns was 5.5 according to the R² knee criterion and 7.25 according to the R² shuffle criterion (see Table 3 for individual values). Both dimensionality estimates were larger than the minimum number of generators required for generating a space with the same linear dimensions as the number of torque spatiotemporal generators (2.75) by non-negative combinations (3.75 = 2.75 + 1). Moreover, differently from torques, the product of the spatial and temporal muscle pattern dimensionalities (22.5 according to the R² knee criterion and 18.7 according to the R² shuffle criterion) was much higher than the spatiotemporal dimensionality. Thus the spatiotemporal generators captured asynchronous muscle coordination patterns that were not simply the result of the synchronous combination of all possible spatial and temporal generators.

Figure 10A shows the R² curve for the spatiotemporal decomposition of the phasic muscle patterns of subject 1 and Figure 10B the seven spatiotemporal generators selected in that subject according to the R² shuffle criterion. The asynchronous nature of the muscle activation waveforms can be noticed in most of these spatiotemporal generators. For example, in TrLat and TrLong in w₁show two clearly delayed peaks. Finally, the examples of muscle pattern reconstructions and combination coefficients for six movement conditions illustrate how the organization of the muscle patterns is captured parsimoniously by the spatiotemporal generators as each movement is reconstructed specifying only 7 scalar combination coefficients (represented by the height of the rectangles depicting the mean generator waveform over all muscles).

FIGURE 10

Figure 10. Spatiotemporal decomposition of phasic muscle patterns. (A) R² curve for subject 1 obtained by spatiotemporal decomposition using NMF. (B) Seven spatiotemporal generators selected for subject 1. (C) Example of the reconstructions of the muscle patters for six movement conditions of subject 1 obtained with the generator illustrated in panel (B) (shaded area: original data, thick line: reconstructed data, bottom: combination coefficients represented by the height of the rectangle containing the temporal profile of each generators averaged over muscles).

Finally, muscle patterns of different subjects did not have similar spatiotemporal generators. The reconstruction of the data by generators of all other subjects had a much lower R²-values (0.37 ± 0.08, mean ±SD) than the R²-values for the reconstruction of the data by generators extracted from the same dataset (0.80 ± 0.03) but still higher than the R²-values obtained with generators identified from randomly shuffled data.

Effect of setting to zero negative values in the phasic muscle patterns

To identify muscle pattern generators from phasic muscle patterns using NMF we set to zero all negative values resulting from the subtraction of the tonic muscle activity from the filtered EMG waveforms. However, to assess the effect of such procedure we also extracted muscle pattern generators from the unclipped phasic muscle patterns using a gradient descent iterative algorithm (see Materials and Methods). In all cases the dimensionality of the generators identified with the gradient descent algorithm from the unclipped data was close to the dimensionality of the generators identified by NMF from the clipped data and the generators extracted in the two cases were very similar. The dimensionality of the spatial generators identified from the unclipped data was, on average across subjects, 4.75 (according to both R² knee and R² shuffle criteria) and thus differed only by 0.25 from the dimensionality of the generators identified from the clipped data. For the temporal generators the difference in dimensionality was 0.5 according to the R² knee criterion and 0.25 according to the R² shuffle criterion. Finally, for the spatiotemporal generators the difference was 0.25 according to both criteria. The similarity between the spatial generators identified from the clipped data and the same number of generators identified from the unclipped data, quantified by the mean normalized scalar product between matched pairs of generators, was 0.91 ± 0.09 (mean across subjects ± SD). For all subjects the similarity value was significantly higher that the value expected by chance, i.e., it was above the 95% percentile of the distribution of the similarity values between generators identified from the unclipped data and generators identified from the shuffled clipped data. The similarity between the temporal generators identified from the clipped data and the same number of generators identified from the unclipped data was 0.91 ± 0.05, also significantly higher than chance for all subjects. Finally, for spatiotemporal generators the similarity was 0.90 ± 0.02 and also significantly higher than chance for all subjects. We can thus conclude that clipping to zero the negative values of the phasic muscle patterns affected the dimensionality and the structure of the identified muscle pattern generators only minimally.

Discussion

We assessed the dimensionality of the dynamic joint torques responsible for accelerating and decelerating the arm during point-to-point reaching movements in different directions in the frontal and sagittal planes and the dimensionality of phasic muscle patterns underlying the production of those torques. We used multidimensional factorization techniques, PCA for the torques and NMF for the muscle patterns, to identify generators capturing the spatial, temporal, and spatiotemporal organization of the motor commands. The number of generators selected according to either a threshold in the total data variation explained, or a change in slope in the curve of the variance explained, or the increase in data variation explained adding an additional generator with respect to the increase obtained extracting generators from randomly shuffled data was taken as an estimate of the dimensionality. The spatial dimensionality of the dynamic torques was lower than the number of joints considered, indicating that some of the available spatial coordination patterns were never employed by the CNS when generating the joint torques for this task. A single temporal generator with a biphasic activation profile was identified in all subjects, in accordance and generalizing previous observations on the temporal organization of dynamic torques on a single vertical plane (Gottlieb et al., 1997). However, a higher number of temporal generators may be required to account for more complex changes in joint torques in other types of reaching movements (e.g., slow, egocentric etc., see Lacquaniti et al., 1986). The number of spatiotemporal generators was in most subjects equal to the product of the spatial and temporal dimensionality and their structure indicated that the spatiotemporal organization of the dynamic torques was essentially synchronous, obtained by the temporal modulation of the spatial generators by the biphasic profile of the single temporal generator. In contrast, the spatial, temporal, and spatiotemporal dimensionalities of the phasic muscle patterns were higher than the corresponding torque dimensionality, as expected because of the non-negativity constraints in the combination of muscle pattern generators, but also higher than the minimum number required according to this biomechanical constraint. Moreover, the spatiotemporal dimensionality of the muscle patterns was much lower than the product of their spatial and temporal dimensionality, suggesting that specific asynchronous coordination patterns were used in the generation of muscle patterns. In fact, most of the identified spatiotemporal generators showed peaks of activity in different muscles at different times, i.e., coordination patterns that cannot be captured by the synchronous modulation of one of the spatial generators by one of the temporal generators.

The CNS might generate motor commands by organizing a few generators, basic elements in a modular architecture capturing shared knowledge across tasks and conditions, to reduce the number of parameters required for control (Alessandro et al., 2013; Ruckert and d'Avella, 2013). Evidence for a modular organization of the motor commands has recently come from the observation of low-dimensionality in the muscle patterns recorded in many species, behaviors, and tasks (Tresch et al., 1999; d'Avella et al., 2003, 2006; Hart and Giszter, 2004; Ivanenko et al., 2004; Ting and Macpherson, 2005; Overduin et al., 2008; Muceli et al., 2010; Dominici et al., 2011; Berger et al., 2013) and from neural recordings and stimulation (Saltiel et al., 2001; Ethier et al., 2006; Gentner and Classen, 2006; Hart and Giszter, 2010; Overduin et al., 2012). Intramuscular recordings during isometric contractions have also revealed that the number of basic muscle activation patterns in complex movements is very limited (ter Haar Romeny et al., 1984; van Zuylen et al., 1988). However, motor task and behaviors are accomplished by the joint torques generated by the simultaneous and coordinated activation of many muscles and to understand how a small set of muscle pattern generators may accomplish a task it is necessary to understand the relationship between muscle patterns and joint torques. The transformation between muscle patterns and torques depends on several biomechanical characteristics and constraints. There are more muscles than joints, making motor commands at the level of muscle patterns redundant, i.e., the same torque pattern can be generated by different muscle patterns. Muscles can only pull and their activation can be expressed by non-negative values, thus introducing a fundamental non-negativity constraint in the generation of muscle patterns. These characteristics and constraints affect the potential dimensionality of the joint torques associated to the dimensionality of the underlying muscle patterns. For a linear mapping of muscle activity into force, an assumption that may be true only in specific conditions such as submaximal isometric contractions (Borzelli et al., 2013) but useful for illustrative purposes, because of the non-negativity constraint, at least D + 1 generators are required to span a D dimensional torque space (Davis, 1954; Valero-Cuevas, 2009). In fact, by linearity, the image of the pseudo-inverse transformation of the D dimensional torque space is a D dimensional subspace of muscle space but such subspace is not contained in the positive orthant of the muscle space, i.e., cannot be generated by non-negative activations and an additional dimension in the null space of the linear transformation needs to be used to achieve a non-negative muscle pattern for each torque. In the general non-linear case, the manifold in muscle space containing all the minimum patterns associated to a D dimensional torque space may be already higher than D even before considering the non-negativity constraint. Thus, the dimensionality of the muscle space must be at least D + 1 to generate a D dimensional torque space for the non-negativity constraint and possibly larger. However, because of redundancy, the CNS might use a number of muscle pattern generator larger than the minimum to optimize some other cost in addition to the number of control parameters, such as effort. Since the generation of muscle patterns with a smaller number of generation requires in general more effort, the dimensionality of the muscle pattern generators might result from a trade-off between computational complexity and effort. We found that the muscle pattern dimensionality is indeed larger than the minimum prescribed by non-negativity, likely an effect of non-linearity but possibly also due to a choice of generators capable of achieving the same motions with less effort. Future investigations comparing tasks with similar kinematics but different effort might help clarify this point.

As mentioned in the Introduction, muscle pattern generators can be associated to force-field primitives (Bizzi et al., 1991; Giszter et al., 1993; Kargo and Giszter, 2000a,b; Giszter and Hart, 2013), endpoint force or joint torque generators that depend on joint angles and velocities and are linearly combined using the activation coefficients of the muscle pattern generators. Competence of force-field primitives to generate observed force and kinematic behaviors has been demonstrated in the frog through forward biomechanical simulation (Kargo et al., 2010). Similarly, a forward dynamics simulation of a musculoskeletal model of the human leg has shown that the combination of a small number of muscle pattern generators is sufficient to perform the basic sub-tasks of walking (Neptune et al., 2009). Moreover, a recent simulation study using a musculoskeletal model of the human arm has indicated that spatiotemporal generators adequate to perform reaching movements can be learned through reinforcement (Ruckert and d'Avella, 2013). Here, in contrast, we did not perform forward dynamics simulation to assess the competence of muscle pattern generators and associated force-field primitive to control reaching movements. We focused, instead, onto the dimensionality of recorded muscle patterns and the dimensionality of joint torques estimated from recorded kinematics by inverse dynamics. Force-field primitives are associated one-to-one with muscle pattern generators, i.e., they have the same dimensionality. Because of the non-negativity of muscle pattern generator activation coefficients and the non-linearity in the muscle-to-force mapping, such dimensionality is necessarily higher than the dimensionality of the joint torques, i.e., the number of torque generators necessary to adequately reconstruct the observed torques by linear combinations. As different numbers of muscle pattern generators and force-field primitives can potentially generate the same set of observed torques and thus being equally competent to perform a given behavior, comparing the dimensionality of muscle patterns and joint torque provides additional information on the strategy that the CNS employs to organize a modular control architecture.

We assessed the dimensionality of joint torques and muscle patterns according to three different definitions of generators (spatial, temporal, and spatiotemporal) and we could then also compare, within each dataset, the different types of dimensionality. For the torques we found that the dimensionality of the spatiotemporal generators was equal to the product of the dimensionalities of the spatial and temporal generators, suggesting that such generators can be obtained as the product of spatial and temporal generators (Delis et al., 2014). The spatial dimensionality was two or three, i.e., less than the number of angular DoF involved in the task, indicating that the CNS selected specific coordination strategies already at the kinematic level. The temporal dimensionality was one, indicating, in accordance with previous observations (Gottlieb et al., 1997), that the temporal organization of the dynamic torques is very simple: a bi-phasic profile shared by all joints and movement conditions (but see Lacquaniti et al., 1986 for more complex torque profiles). The dimensionality of the spatiotemporal generators was in all subjects equal to the spatial dimensionality because, given a single temporal generator, each spatiotemporal generator was obtained by the temporal modulation of each spatial generator by the temporal generator. Consequently, the spatiotemporal organization of the torques was essentially synchronous. In contrast, the number of spatiotemporal generators for the muscle patterns was much less than the product of spatial and temporal dimensionalities. Indeed, spatiotemporal generators captured asynchronous activations across muscles that could not be obtained by the modulation of a single spatial generator by a single temporal generator, which necessarily produces a synchronous pattern. Thus, spatiotemporal generators appear to provide a very compact representation of the organization of the muscle patterns (Delis et al., 2014). However, differently from our previous analyses of muscle patterns during reaching (d'Avella et al., 2006, 2008, 2011; d'Avella and Lacquaniti, 2013), in the present spatiotemporal decomposition we did not take into account the possibility of shifting in time the onset of different generators (Equation 8) because in this way we could use the same NMF algorithm used for spatial and temporal decomposition. We found a larger number of generators without time-shifts than the number of time-varying muscle synergies (with time-shifts) reported before. Thus, additional structure in the muscle patterns can be captured by allowing the independent modulation of the time of recruitment of the generators thus allowing for an even more compact representation of the muscle pattern organization.

The validity of our observations depends on a number of assumptions made in the analysis of the torque and muscle activity data. Concerning the estimation of the joint torques from the recorded motions of markers positioned on the subjects' arm, we relied on a simplified model of the human arm. We assumed that the shoulder was a spherical joint, i.e., all three rotation axes intersect at a single point, we estimated the length (Winter, 1990) and inertial parameters (Zatsiorsky and Seluyanov, 1983) of each segment as a function of the height and weight of each subject. To assess the effect of potential inaccuracies in the parameters of our model, we performed inverse dynamics varying the mass and height of each subject up to ±20% and we found that the estimated joint torque dimensionality changed only in a small fraction of cases. Thus, we believe that the estimation of torque dimensionality is robust to small inaccuracies in the anthropometric parameters. Concerning the identification of the muscle pattern generators with NMF, in order to be able to run such algorithm, we set to zero all samples with a negative value after subtracting the tonic components. We assessed the effect of such procedure by identifying muscle pattern generators from unclipped data using a gradient descent algorithm and we found that both the dimensionality changed only minimally (less than 0.5 in all cases) and did not significantly affect the structure of the generators. Concerning the criteria for the selection of the number of generators, we used a threshold based criterion for torques, a criterion based on the detection of a change in slope (i.e., a “knee”) in the R² curve for muscle patterns, and a criterion based on the comparison of the slope of the R² curve for the generators extracted from the original data and those extracted after randomly shuffling the data (along the rows of the data matrix). All these criteria rely on the assumption that the data are generated by a number of generators smaller than the maximum dimensionality and that a fraction of the variation observed in the data is due to noise. Such assumption is shared by all previous studies using multidimensional decomposition to identify muscle synergies or temporal components but it is clear that it is not possible a-priori to exclude that, once a specific number of generators has been selected, the additional dimensions attributed to noise might be also necessary to capture the structure in the motor commands or, vice-versa, a generator might actually describe noise instead of structure in the motor commands. Moreover, unless an independent estimation of the level of noise in the data is available, the selection of the number of generators depends on ad-hoc choices of thresholds and parameters. However, as the determination of threshold on the MSE of the linear fit of the terminal portion of the R² curve, used for the R² knee criterion, is less dependent on the amount noise in the data than the threshold on the R²-value, used in the R² threshold criterion, we prefer, whenever possible, to use the former criterion, as we have done previously (d'Avella et al., 2003, 2006; Cheung et al., 2005), to the latter, also used in previous studies (Tresch et al., 1999; Ting and Macpherson, 2005; Torres-Oviedo et al., 2006; Roh et al., 2012). Unfortunately it was impossible to use the R² knee criterion for the spatial torque dimensionality, as the maximum dimension was 4 in that case, and we could only estimate the residual of a linear fit of the R² curve from 1 or 2 to 4 generators, i.e., evaluating the presence of a knee only at 1 or 2 generators. We then used the third criterion to compare more directly both datasets. The R² shuffle criterion (Cheung et al., 2009) is based on the assumption that the multidimensional structure but not the noise is affected by randomly shuffling the data. Thus, the selection of the number of generators in datasets with different amount of noise is based on the comparison of the slope of the R² curve for each data set with the slope of the R² curve of a random dataset with comparable level of noise. In most cases the number of generators selected by the two criteria either matched exactly or different by one.

In conclusion, whether spatial (time-invariant muscle synergies), temporal (temporal components or patterns), or spatiotemporal (time-varying muscle synergies) generators are fundamental building blocks in a modular control architecture and how are they implemented in the CNS remain open and debated questions. Our comparison of the dimensionality of muscle patterns and joint torques suggests that the larger dimensionalities and spatiotemporal complexity of the muscle patterns with respect to the joint torques may be required for the CNS to overcome the non-linearities of the musculoskeletal system and, exploiting its redundancy, to flexibly generate endpoint trajectories with simple kinematic features using a limited number of building blocks.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

Supported by the Italian Ministry of Health, the Italian Space Agency (DCMC, CRUSOE, and COREA grants), and the EU Seventh Framework Programme (FP7-ICT No 248311 AMARSi).

References

Alessandro, C., Carbajal, J. P., and d'Avella, A. (2013). A computational analysis of motor synergies by dynamic response decomposition. Front. Comput. Neurosci. 7:191. doi: 10.3389/fncom.2013.00191

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Berger, D. J., Gentner, R., Edmunds, T., Pai, D. K., and d'Avella, A. (2013). Differences in adaptation rates after virtual surgeries provide direct evidence for modularity. J. Neurosci. 33, 12384–12394. doi: 10.1523/JNEUROSCI.0122-13.2013

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bizzi, E., Cheung, V. C., d'Avella, A., Saltiel, P., and Tresch, M. (2008). Combining modules for movement. Brain Res. Rev. 57, 125–133. doi: 10.1016/j.brainresrev.2007.08.004

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bizzi, E., Mussa-Ivaldi, F. A., and Giszter, S. (1991). Computations underlying the execution of movement: a biological perspective. Science 253, 287–291. doi: 10.1126/science.1857964

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Borzelli, D., Berger, D. J., Pai, D. K., and d'Avella, A. (2013). Effort minimization and synergistic muscle recruitment for three-dimensional force generation. Front. Comput. Neurosci. 7:186. doi: 10.3389/fncom.2013.00186

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cheung, V. C., d'Avella, A., Tresch, M. C., and Bizzi, E. (2005). Central and sensory contributions to the activation and organization of muscle synergies during natural motor behaviors. J. Neurosci. 25, 6419–6434. doi: 10.1523/JNEUROSCI.4904-04.2005

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cheung, V. C., Piron, L., Agostini, M., Silvoni, S., Turolla, A., and Bizzi, E. (2009). Stability of muscle synergies for voluntary actions after cortical stroke in humans. Proc. Natl. Acad. Sci. U.S.A. 106, 19563–19568. doi: 10.1073/pnas.0910114106

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Chhabra, M., and Jacobs, R. A. (2006). Properties of synergies arising from a theory of optimal motor behavior. Neural Comput. 18, 2320–2342. doi: 10.1162/neco.2006.18.10.2320

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Corke, P. I. (1996). A robotics toolbox for MATLAB. IEEE Robot. Autom. Mag. 3, 24–32. doi: 10.1109/100.486658

CrossRef Full Text

Corke, P. I. (ed.). (2011). Robotics, Vision and Control. Berlin Heidelberg: Springer-Verlag. doi: 10.1007/978-3-642-20144-8

CrossRef Full Text

d'Avella, A., Fernandez, L., Portone, A., and Lacquaniti, F. (2008). Modulation of phasic and tonic muscle synergies with reaching direction and speed. J. Neurophysiol. 100, 1433–1454. doi: 10.1152/jn.01377.2007

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

d'Avella, A., and Lacquaniti, F. (2013). Control of reaching movements by muscle synergy combinations. Front. Comput. Neurosci. 7:42. doi: 10.3389/fncom.2013.00042

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

d'Avella, A., Portone, A., Fernandez, L., and Lacquaniti, F. (2006). Control of fast-reaching movements by muscle synergy combinations. J. Neurosci. 26, 7791–7810. doi: 10.1523/JNEUROSCI.0830-06.2006

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

d'Avella, A., Portone, A., and Lacquaniti, F. (2011). Superposition and modulation of muscle synergies for reaching in response to a change in target location. J. Neurophysiol. 106, 2796–2812. doi: 10.1152/jn.00675.2010

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

d'Avella, A., Saltiel, P., and Bizzi, E. (2003). Combinations of muscle synergies in the construction of a natural motor behavior. Nat. Neurosci. 6, 300–308. doi: 10.1038/nn1010

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Davis, C. (1954). Theory of positive linear dependence. Am. J. Math. 76, 733–746. doi: 10.2307/2372648

CrossRef Full Text

Delis, I., Berret, B., Pozzo, T., and Panzeri, S. (2013). Quantitative evaluation of muscle synergy models: a single-trial task decoding approach. Front. Comput. Neurosci. 7:8. doi: 10.3389/fncom.2013.00008

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Delis, I., Panzeri, S., Pozzo, T., and Berret, B. (2014). A unifying model of concurrent spatial and temporal modularity in muscle activity. J. Neurophysiol. 111, 675–693. doi: 10.1152/jn.00245.2013

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Dominici, N., Ivanenko, Y. P., Cappellini, G., d'Avella, A., Mondi, V., Cicchese, M., et al. (2011). Locomotor primitives in newborn babies and their development. Science 334, 997–999. doi: 10.1126/science.1210617

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ethier, C., Brizzi, L., Darling, W. G., and Capaday, C. (2006). Linear summation of cat motor cortex outputs. J. Neurosci. 26, 5574–5581. doi: 10.1523/JNEUROSCI.5332-05.2006

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Flanders, M., and Herrmann, U. (1992). Two components of muscle activation: scaling with the speed of arm movement. J. Neurophysiol. 67, 931–943.

Pubmed Abstract | Pubmed Full Text

Flanders, M., Pellegrini, J. J., and Geisler, S. D. (1996). Basic features of phasic activation for reaching in vertical planes. Exp. Brain Res. 110, 67–79. doi: 10.1007/BF00241376

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Flanders, M., Pellegrini, J. J., and Soechting, J. F. (1994). Spatial/temporal characteristics of a motor pattern for reaching. J. Neurophysiol. 71, 811–813.

Pubmed Abstract | Pubmed Full Text

Gentner, R., and Classen, J. (2006). Modular organization of finger movements by the human central nervous system. Neuron 52, 731–742. doi: 10.1016/j.neuron.2006.09.038

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Giszter, S., Patil, V., and Hart, C. (2007). Primitives, premotor drives, and pattern generation: a combined computational and neuroethological perspective. Prog. Brain Res. 165, 323–346. doi: 10.1016/S0079-6123(06)65020-6

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Giszter, S. F., and Hart, C. B. (2013). Motor primitives and synergies in the spinal cord and after injury—the current state of play. Ann. N.Y. Acad. Sci. 1279, 114–126. doi: 10.1111/nyas.12065

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Giszter, S. F., Mussa-Ivaldi, F. A., and Bizzi, E. (1993). Convergent force fields organized in the frog's spinal cord. J. Neurosci. 13, 467–491.

Pubmed Abstract | Pubmed Full Text

Gottlieb, G. L., Song, Q., Almeida, G. L., Hong, D., and Corcos, D. (1997). Directional control of planar human arm movement. J. Neurophysiol. 78, 2985–2998.

Pubmed Abstract | Pubmed Full Text

Hart, C. B., and Giszter, S. F. (2004). Modular premotor drives and unit bursts as primitives for frog motor behaviors. J. Neurosci. 24, 5269–5282. doi: 10.1523/JNEUROSCI.5626-03.2004

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hart, C. B., and Giszter, S. F. (2010). A neural basis for motor primitives in the spinal cord. J. Neurosci. 30, 1322–1336. doi: 10.1523/JNEUROSCI.5894-08.2010

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hartenberg, R. S., and Denavit, J. (1955). A kinematic notation for lower pair mechanisms based on matrices. J. Appl. Mech. 77, 215–221.

Ivanenko, Y. P., Poppele, R. E., and Lacquaniti, F. (2004). Five basic muscle activation patterns account for muscle activity during human locomotion. J. Physiol. 556, 267–282. doi: 10.1113/jphysiol.2003.057174

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kargo, W. J., and Giszter, S. F. (2000a). Afferent roles in hindlimb wipe-reflex trajectories: free-limb kinematics and motor patterns. J. Neurophysiol. 83, 1480–1501.

Pubmed Abstract | Pubmed Full Text

Kargo, W. J., and Giszter, S. F. (2000b). Rapid correction of aimed movements by summation of force-field primitives. J. Neurosci. 20, 409–426.

Pubmed Abstract | Pubmed Full Text

Kargo, W. J., Ramakrishnan, A., Hart, C. B., Rome, L. C., and Giszter, S. F. (2010). A simple experimentally based model using proprioceptive regulation of motor primitives captures adjusted trajectory formation in spinal frogs. J. Neurophysiol. 103, 573–590. doi: 10.1152/jn.01054.2007

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Krishnamoorthy, V., Goodman, S., Zatsiorsky, V., and Latash, M. L. (2003). Muscle synergies during shifts of the center of pressure by standing persons: identification of muscle modes. Biol. Cybern. 89, 152–161. doi: 10.1007/s00422-003-0419-5

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lacquaniti, F., Ivanenko, Y. P., and Zago, M. (2012). Patterned control of human locomotion. J. Physiol. 590, 2189–2199. doi: 10.1113/jphysiol.2011.215137

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lacquaniti, F., Soechting, J. F., and Terzuolo, C. A. (1982). Some factors pertinent to the organization and control of arm movements. Brain Res. 252, 394–397. doi: 10.1016/0006-8993(82)90410-3

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lacquaniti, F., Soechting, J. F., and Terzuolo, S. A. (1986). Path constraints on point-to-point arm movements in three-dimensional space. Neuroscience 17, 313–324. doi: 10.1016/0306-4522(86)90249-6

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Morasso, P. (1981). Spatial control of arm movements. Exp. Brain Res. 42, 223–227. doi: 10.1007/BF00236911

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Muceli, S., Boye, A. T., d'Avella, A., and Farina, D. (2010). Identifying representative synergy matrices for describing muscular activation patterns during multidirectional reaching in the horizontal plane. J. Neurophysiol. 103, 1532–1542. doi: 10.1152/jn.00559.2009

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Mussa-Ivaldi, F. A., Giszter, S. F., and Bizzi, E. (1994). Linear combinations of primitives in vertebrate motor control. Proc. Natl. Acad. Sci. U.S.A. 91, 7534–7538. doi: 10.1073/pnas.91.16.7534

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Neptune, R. R., Clark, D. J., and Kautz, S. A. (2009). Modular control of human walking: a simulation study. J. Biomech. 42, 1282–1287. doi: 10.1016/j.jbiomech.2009.03.009

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Overduin, S. A., d'Avella, A., Carmena, J. M., and Bizzi, E. (2012). Microstimulation activates a handful of muscle synergies. Neuron 76, 1071–1077. doi: 10.1016/j.neuron.2012.10.018

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Overduin, S. A., d'Avella, A., Roh, J., and Bizzi, E. (2008). Modulation of muscle synergy recruitment in primate grasping. J. Neurosci. 28, 880–892. doi: 10.1523/JNEUROSCI.2869-07.2008

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Roh, J., Rymer, W. Z., and Beer, R. F. (2012). Robustness of muscle synergies underlying three-dimensional force generation at the hand in healthy humans. J. Neurophysiol. 107, 2123–2142. doi: 10.1152/jn.00173.2011

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ruckert, E., and d'Avella, A. (2013). Learned parametrized dynamic movement primitives with shared synergies for controlling robotic and musculoskeletal systems. Front. Comput. Neurosci. 7:138. doi: 10.3389/fncom.2013.00138

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Saltiel, P., Wyler-Duda, K., d'Avella, A., Tresch, M. C., and Bizzi, E. (2001). Muscle synergies encoded within the spinal cord: evidence from focal intraspinal NMDA iontophoresis in the frog. J. Neurophysiol. 85, 605–619.

Pubmed Abstract | Pubmed Full Text

Soechting, J. F., and Lacquaniti, F. (1981). Invariant characteristics of a pointing movement in man. J. Neurosci. 1, 710–720.

Pubmed Abstract | Pubmed Full Text

ter Haar Romeny, B. M., van der Gon, J. J., and Gielen, C. C. (1984). Relation between location of a motor unit in the human biceps brachii and its critical firing levels for different tasks. Exp. Neurol. 85, 631–650. doi: 10.1016/0014-4886(84)90036-0

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ting, L. H., and Macpherson, J. M. (2005). A limited set of muscle synergies for force control during a postural task. J. Neurophysiol. 93, 609–613. doi: 10.1152/jn.00681.2004

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ting, L. H., and McKay, J. L. (2007). Neuromechanics of muscle synergies for posture and movement. Curr. Opin. Neurobiol. 17, 622–628. doi: 10.1016/j.conb.2008.01.002

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Torres-Oviedo, G., Macpherson, J. M., and Ting, L. H. (2006). Muscle synergy organization is robust across a variety of postural perturbations. J. Neurophysiol. 96, 1530–1546. doi: 10.1152/jn.00810.2005

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Tresch, M. C., Cheung, V. C., and d'Avella, A. (2006). Matrix factorization algorithms for the identification of muscle synergies: evaluation on simulated and experimental data sets. J. Neurophysiol. 95, 2199–2212. doi: 10.1152/jn.00222.2005

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Tresch, M. C., and Jarc, A. (2009). The case for and against muscle synergies. Curr. Opin. Neurobiol. 19, 601–607. doi: 10.1016/j.conb.2009.09.002

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Tresch, M. C., Saltiel, P., and Bizzi, E. (1999). The construction of movement by the spinal cord. Nat. Neurosci. 2, 162–167. doi: 10.1038/5721

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Valero-Cuevas, F. J. (2009). A mathematical approach to the mechanical capabilities of limbs and fingers. Adv. Exp. Med. Biol. 629, 619–633. doi: 10.1007/978-0-387-77064-2_33

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

van Zuylen, E. J., Gielen, C. C., and Denier van der Gon, J. J. (1988). Coordination and inhomogeneous activation of human arm muscles during isometric torques. J. Neurophysiol. 60, 1523–1548.

Pubmed Abstract | Pubmed Full Text

Winter, D. A. (ed.). (1990). Biomechanics and Motor Control of Human Movement. New York, NY: Wiley.

Zatsiorsky, V., and Seluyanov, V. (1983). “The mass and inertia characteristics of the main segments of the human body,” in Biomechanics, eds H. Matsui and K. Kobayashi (Champain, IL: Human Kinetics), 7791–7810.