Edited by: Martin Giese, University Clinic Tuebingen, Germany
Reviewed by: Simon Giszter, Drexel Med School, USA; Stan Gielen, Radboud University Nijmegen, Netherlands
*Correspondence: Andrea d'Avella, Laboratory of Neuromotor Physiology, Santa Lucia Foundation, Via Ardeatina 306, 00179 Rome, Italy e-mail:
This article was submitted to the journal Frontiers in Computational Neuroscience.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
Muscle activities underlying many motor behaviors can be generated by a small number of basic activation patterns with specific features shared across movement conditions. Such low-dimensionality suggests that the central nervous system (CNS) relies on a modular organization to simplify control. However, the relationship between the dimensionality of muscle patterns and that of joint torques is not fixed, because of redundancy and non-linearity in mapping the former into the latter, and needs to be investigated. We compared the torques acting at four arm joints during fast reaching movements in different directions in the frontal and sagittal planes and the underlying muscle patterns. The dimensionality of the non-gravitational components of torques and muscle patterns in the spatial, temporal, and spatiotemporal domains was estimated by multidimensional decomposition techniques. The spatial organization of torques was captured by two or three generators, indicating that not all the available coordination patterns are employed by the CNS. A single temporal generator with a biphasic profile was identified, generalizing previous observations on a single plane. The number of spatiotemporal generators was equal to the product of the spatial and temporal dimensionalities and their organization was essentially synchronous. Muscle pattern dimensionalities were higher than torques dimensionalities but also higher than the minimum imposed by the inherent non-negativity of muscle activations. The spatiotemporal dimensionality of the muscle patterns was lower than the product of their spatial and temporal dimensionality, indicating the existence of specific asynchronous coordination patterns. Thus, the larger dimensionalities of the muscle patterns may be required for CNS to overcome the non-linearities of the musculoskeletal system and to flexibly generate endpoint trajectories with simple kinematic features using a limited number of building blocks.
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.,
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,
To generate the joint torques
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
The torque does not depend in general linearly on the muscle activations and, consequently, on the combination coefficients
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
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,
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.,
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.
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
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,
The rotation axis of each joint coincides with the
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,
1 | Sh X | π/2 | 0 | π/2 | 0 | P | 0 |
2 | Sh Y | π/2 | 0 | π/2 | 0 | P | 0 |
3 | Sh Z | π/2 | 0 | π/2 | 0 | P | 0 |
4 | Sh adduction | π/2 | 0 | 0 | 0 | R | −π/2 |
5 | Sh flexion | π/2 | 0 | 0 | 0 | R | π/2 |
6 | Sh external rotation | π/2 | 0 | 0 | LU | R | π |
7 | El flexion | 0 | LF | 0 | 0 | R | π/2 |
Height (cm) | 180 | 162 | 177 | 181 |
Weight (Kg) | 84 | 58 | 75 | 78 |
LU (cm) | 33.48 | 30.13 | 32.92 | 33.67 |
LF (cm) | 45.36 | 40.82 | 44.60 | 45.61 |
rU (cm) | 13.91 | 12.16 | 13.48 | 13.78 |
rF (cm) | 26.38 | 23.57 | 25.83 | 26.39 |
MU (kg) | 2.29 | 1.56 | 2.03 | 2.11 |
MF (kg) | 2.00 | 1.52 | 1.84 | 1.90 |
I(lo) U (kg cm2 s−2) | 46.54 | 28.54 | 40.45 | 42.61 |
I(ap) U (kg cm2 s−2) | 137.84 | 74.02 | 120.09 | 130.03 |
I(tr) U (kg cm2 s−2) | 152.50 | 84.74 | 133.92 | 144.65 |
I(lo) F (kg cm2 s−2) | 21.91 | 12.93 | 18.63 | 19.56 |
I(ap) F (kg cm2 s−2) | 465.00 | 295.87 | 416.54 | 445.74 |
I(tr) F (kg cm2 s−2) | 475.79 | 302.27 | 425.77 | 455.45 |
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
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.
Joint angles, joint velocities and joint accelerations were used to estimate the torque profiles via recursive Newton-Euler calculation (
To validate the inverse dynamics calculation we also performed a forward dynamics simulation (
The EMGs for each trial were digitally full-wave rectified, low-pass filtered (FIR filter, 20 Hz cut-off, zero-phase distortion, Matlab
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.
We consider a set of
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.,
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.
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 (
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
S 1 | R2 threshold | 2 | 1 | 3 | – | – | – |
R2 knee | – | – | – | 4 | 5 | 5 | |
R2 shuffle | 3 | 1 | 3 | 4 | 4 | 7 | |
S 2 | R2 threshold | 3 | 1 | 3 | – | – | – |
R2 knee | – | – | – | 5 | 5 | 6 | |
R2 shuffle | 3 | 1 | 3 | 5 | 4 | 7 | |
S 3 | R2 threshold | 2 | 1 | 3 | – | – | – |
R2 knee | – | – | – | 6 | 4 | 6 | |
R2 shuffle | 3 | 1 | 3 | 6 | 3 | 7 | |
S 4 | R2 threshold | 2 | 1 | 2 | – | – | – |
R2 knee | – | – | – | 5 | 4 | 5 | |
R2 shuffle | 2 | 1 | 2 | 5 | 4 | 8 |
Figure
Figure
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
Spatial | S 1 | 0.99 | 0.99 | 0.98 | 0.99 |
S 2 | 0.99 | 0.99 | 0.97 | 0.99 | |
S 3 | 0.99 | 0.98 | 0.99 | 0.99 | |
S 4 | 0.92 | 0.90 | 0.88 | 0.95 | |
Temporal | S 1 | 0.95 | 0.92 | 0.92 | 0.93 |
S 2 | 0.94 | 0.93 | 0.94 | 0.93 | |
S 3 | 0.91 | 0.92 | 0.94 | 0.92 | |
S 4 | 0.94 | 0.92 | 0.93 | 0.93 | |
Spatiotemporal | S 1 | 0.97 | 0.92 | 0.90 | 0.94 |
S 2 | 0.94 | 0.95 | 0.92 | 0.93 | |
S 3 | 0.90 | 0.89 | 0.96 | 0.92 | |
S 4 | 0.88 | 0.84 | 0.82 | 0.91 |
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
Figure
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
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
pt Figure
As in previous cases, spatiotemporal generators were similar across subjects. The
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.
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.,
The mean spatial dimensionality of the phasic muscle patterns across subjects was five according to the position of change in slope of the R2 curve as a function of the number of generators (R2 knee criterion) and five according to the R2 shuffle criterion (see Table
Figure
Finally the spatial generators for the muscle patterns were less similar across subjects than the spatial generators for the torques (see Table
Spatial | S 1 | 0.80 | 0.60 | 0.54 | 0.62 |
S 2 | 0.64 | 0.74 | 0.61 | 0.70 | |
S 3 | 0.67 | 0.71 | 0.80 | 0.69 | |
S 4 | 0.73 | 0.75 | 0.60 | 0.85 | |
Temporal | S 1 | 0.88 | 0.83 | 0.88 | 0.83 |
S 2 | 0.85 | 0.86 | 0.88 | 0.85 | |
S 3 | 0.81 | 0.79 | 0.86 | 0.79 | |
S 4 | 0.77 | 0.78 | 0.83 | 0.81 | |
Spatiotemporal | S 1 | 0.84 | 0.30 | 0.27 | 0.32 |
S 2 | 0.37 | 0.77 | 0.41 | 0.51 | |
S 3 | 0.25 | 0.41 | 0.79 | 0.36 | |
S 4 | 0.38 | 0.49 | 0.39 | 0.80 |
The mean number of temporal generators of the phasic muscle patterns was 4.5 according to the R2 knee criterion and 3.75 according to the R2 shuffle criterion (see Table
Figure
In contrast to the spatial generators but similarly to the temporal generators for torques, muscle pattern temporal generators were similar across all subjects. The
The mean number of spatiotemporal generators of the phasic muscle patterns was 5.5 according to the R2 knee criterion and 7.25 according to the R2 shuffle criterion (see Table
Figure
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
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 R2 knee and R2 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 R2 knee criterion and 0.25 according to the R2 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 ±
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.,
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.,
As mentioned in the Introduction, muscle pattern generators can be associated to force-field primitives (Bizzi et al.,
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.,
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,
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.
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.
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).