Low-Dimensional Synergistic Representation of Bilateral Reaching Movements

Kinematic and neuromuscular synergies have been found in numerous aspects of human motion. This study aims to determine how effectively kinematic synergies in bilateral upper arm movements can be used to replicate complex activities of daily living (ADL) tasks using a sparse optimization algorithm. Ten right-handed subjects executed 18 rapid and 11 natural-paced ADL tasks requiring bimanual coordination while sitting at a table. A position tracking system was used to track the subjects’ arms in space, and angular velocities over time for shoulder abduction, shoulder flexion, shoulder internal rotation, and elbow flexion for each arm were computed. Principal component analysis (PCA) was used to generate kinematic synergies from the rapid-paced task set for each subject. The first three synergies accounted for 80.3 ± 3.8% of variance, while the first eight accounted for 94.8 ± 0.85%. The first and second synergies appeared to encode symmetric reaching motions which were highly correlated across subjects. The first three synergies were correlated between left and right arms within subjects, whereas synergies four through eight were not, indicating asymmetries between left and right arms in only the higher order synergies. The synergies were then used to reconstruct each natural-paced task using the l1-norm minimization algorithm. Temporal dilations of the synergies were introduced in order to model the temporal scaling of movement patterns achieved by the cerebellum and basal ganglia as reported previously in the literature. Reconstruction error was reduced by introducing synergy dilations, and cumulative recruitment of several synergies was significantly reduced in the first 10% of training task time by introducing temporal dilations. The outcomes of this work could open new scenarios for the applications of postural synergies to the control of robotic systems, with potential applications in rehabilitation. These synergies not only help in providing near-natural control but also provide simplified strategies for design and control of artificial limbs. Potential applications of these bilateral synergies were discussed and future directions were proposed.


INTRODUCTION
The human arm is a highly complex structure with an equally sophisticated control system. Each arm possesses 11 independent degrees of freedom (DoF) defined from the pectoral girdle to the wrist, which are actuated by approximately 32 muscles (Mackenzie and Iberall, 1994). The brain, therefore, has to coordinate over 60 different controls in order to operate both arms, yet accomplishes this task with apparent ease. How the brain handles real-time control of two redundant, high DoF manipulators during activities of daily living (ADL) is known as the degrees-of-freedom (DoF) problem (Bernstein, 1967;Latash et al., 2007) and is the subject of much research, including the present work. Progress in this field has applications in numerous areas including motor rehabilitation, assistive and prosthetic technology, and robotic control.
Evidence suggests that the brain may control the limbs by scaling, offsetting, and temporally dilating fundamental movements encoded in the sensorimotor system (Viviani and Terzuolo, 1980;Brooks, 1986). Previous research has shown that the brain executes tasks by using certain movement patterns while preserving their relative spatiotemporal proportions. Viviani and Terzuolo (1980) have shown in handwriting tasks that increased letter size still results in similar execution times by automatically increasing writing speed. Furthermore, slowing down a writing task results in temporal dilation of a common velocity pattern, preserving the relative occurrence of velocity profile features in time while reducing velocity amplitude (Brooks, 1986).
These patterns of motion have been developed into the concept of synergies, which can be defined as "a collection of relatively independent degrees of freedom that behave as a single functional unit" (Turvey, 2007). Synergies exist in either the joint angular velocity space, in the form of kinematic synergies, or neuromuscular activity space, in the form of neuromuscular synergies. Linear discriminant analysis, singular value decomposition (SVD), principal component analysis (PCA), non-negative matrix factorization (NMF), artificial neural networks, and many other algorithms have been used in the literature to derive synergies for hand grasps, gait patterns, and single-arm motion (Merckle et al., 1998;Santello et al., 1998;Vinjamuri et al., 2010;Roh et al., 2013;Alibeji et al., 2015). NMF is typically used to derive neuromuscular synergies (Tresch and Jarc, 2009), while PCA is frequently used to derive kinematic synergies as in Mason et al. (2001) and Vinjamuri et al. (2010). PCA-derived kinematic synergies have been demonstrated to perform favorably when directly compared to those from other linear and non-linear dimensionality reduction methods when applied to hand grasp reconstruction (Patel et al., 2015a).
Recent work has been aimed at integrating synergies into the control of robotic systems with the goal of producing a simplified control scheme for high DoF devices. The authors in Chen et al. (2015) have demonstrated an anthropomorphic robotic hand that has two mechanically implemented postural synergies which could successfully grasp various objects. Several groups have also proposed autonomous, control systems for high DoF robotic and virtual hands based on two postural synergies (Wimbock et al., 2011) and four postural synergies (Rosell et al., 2011;Segil and Weir, 2013), whereas Matrone et al. (2012) have demonstrated real-time myoelectric control of a robotic hand using two postural synergies with able-bodied subjects. An EMG-based control scheme was also introduced by Artemiadis and Kyriakopoulos (2010), which controls a 7-DoF robotic arm using kinematic and muscular synergies. The review recently published by Santello et al. (2016) gives a thorough description of the state of the art concerning dexterous hand control using synergies and highlights some future directions merging synergies with compliant design.
Synergies derived using NMF have also been applied to optimal movement generation for virtual arms (Fu et al., 2013) as well as myocontrol of a multi-DoF planar robotic arm using muscle synergies (Lunardini et al., 2015). So far, work has been focused on using time-invariant postural synergies in the kinematic domain and restricted to unimanual processes.
Bilateral spatiotemporal kinematic synergies such as those presented here may be used as the controlled variable in future robotic systems that can be manipulated using EMG, EEG, or some other biosignal input. Whereas postural/spatial synergies attempt to linearize joint motion relative to each other, a timevarying approach allows more flexibility to capture the non-linear behaviors inherent to complex systems. An open question for such a system is whether or not ADL tasks are within the "workspace" of a system that is only manipulated using time-varying kinematic synergies. In other words, is it possible to manipulate bilateral spatiotemporal kinematic synergies by scaling their amplitudes and temporal offsets in such a way as to replicate ADL-like tasks.
In this study, we derive spatiotemporal kinematic synergies from rapidly paced ADL tasks for 10 able-bodied subjects. Tasks that require coordination of both arms and can be classified as symmetric in-phase, symmetric out-of-phase, asymmetric, or coupled are chosen. PCA is used to derive time-varying kinematic synergies from eight joint angular velocity profiles across both arms recorded during these rapid tasks. A separate set of tasks performed at a natural pace are reconstructed using the l 1 -norm minimization algorithm to select optimal amplitudes and temporal offsets and dilations of these synergies. The derived synergies are characterized in terms of intersubject and interlimb correlations, accuracy of reconstruction, and trends in their recruitment levels throughout the task duration.

SUBJECTS AND METHODS
The present study was conducted under IRB Approved Protocol # 2014-026/2015-022 at the Stevens Institute of Technology. Ten subjects were recruited in the study after obtaining written informed consent. Subjects performed ADL-like tasks while their movements were recorded using an electromagnetic motion tracking system (Polhemus LIBERTY). Positional data from each sensor were converted into joint angles, synergies were derived using PCA in the joint angle velocity domain, and a separate set of tasks were reconstructed from the derived synergies using the l 1 -norm minimization algorithm.

Data Capture
An electromagnetic tracking system (Polhemus LIBERTY, TX4 source) was used to record positional data of the subject during each task using their proprietary software (PiMgr). The study was executed in a minimal-metal environment with a compensation map calibration executed monthly to account for disturbances due to metal in the construction of the room. The workspace was calibrated such that the origin was on the edge of the table, centered in front of the subject. Positive Z extended upwards toward the ceiling, positive X extended forward away from the subject, and positive Y extended to the subject's left. Data were captured at 240 Hz and filtered using a 3 Hz fourth-order low-pass Butterworth filter.

Kinematic Model
Several groups have developed refined anatomical models with the intent of capturing kinematic data from subjects as they perform tasks. Most of this work has addressed hurdles using optical systems such as 3D interpolation of one or multiple 2D viewpoints (Sidenbladh et al., 2000;Chen et al., 2010), and soft tissue deformation (Gabiccini et al., 2013). Under guidance from Wu et al. (2005), this work utilized an electromagnetic tracking system to capture the positions of several convenient landmarks on the torso, left, and right arms with a positional accuracy of approximately 2.5 mm in x, y, and z. The tracking system lacks line of sight issues and readily supplies Cartesian positions of these landmarks.
Seven sensors were placed on the body as shown in Figure 1A. Three sensors defined the trunk of the subject, while two additional sensors per arm tracked elbow and wrist movements in space. S 1 and S 2 were placed at the lateral head of the clavicle on the subject's right and left shoulder, respectively, while S 3 was placed on the subject's right side near the middle of the rib cage on the midaxillary line. S 4 and S 6 were placed on the lateral side of the subject's elbows over the joint's center of rotation. S 5 and S 7 were placed on the dorsal side of the subject's wrists and were centered between the distal ulnar and radial heads. The filtered X, Y, and Z trajectories captured by the tracking system were converted to joint angles as follows.
Three shoulder angles and one elbow angle were calculated for each side of the subject: shoulder abduction/adduction, flexion/extension, and internal/external and elbow flexion/extension. These angles are calculated using six vectors: V should , V side , V ae , and V aw , where a = L, R to indicate the left or right arm. V should is a vector from the subject's left shoulder to their right shoulder sensors, V side is a vector from the right shoulder sensor to the sensor on the right side of the torso, V ae is a vector from the shoulder to the elbow on each respective side, and V aw is the vector from the elbow to the wrist on each respective side. These vectors were calculated as: Abduction/adduction was found by first projecting the elbow vector, V ae , onto the coronal plane: where n c is the vector normal to the coronal plane and is defined as the cross product between V should and V side . The shoulder abduction angle, θ sa , between V ae,proj and V should was found using The flexion/extension angle was found in the same manner as in Eqs 1 and 2 except by projecting the vectors onto the sagittal plane using n s , the normal vector to the sagittal plane, which was found by normalizing V should . The internal/external rotation of the upper FIGURE 1 | (A) Sensor placement on body. S 1 and S 2 are positioned at the lateral head of the clavicles, S3 is placed on the right side on the body's midaxillary line, S 4 and S 6 are placed on the outer side of the elbow, and S 5 and S 7 are placed on the wrists between the distal heads of the radius and ulna. P 1 indicates the coronal plane, P 2 indicates the sagittal plane, P 3 indicates the plane normal to Vae, and P 4 indicates the plane containing Vae and Vaw. P 3 and P 4 exist for both the left and right arms. (B) Tasks executed during study. Each panel is labeled with a number corresponding to the task in Table 2 which is shown. Subjects start each task with their hands in the rest positions marked by the visible red boxes.
Frontiers in Bioengineering and Biotechnology | www.frontiersin.org February 2017 | Volume 5 | Article 2 arm was found by projecting the wrist vector, V aw , and V should onto the plane normal to V ae and calculating the angle between the two projections as in Eq. 2. Elbow flexion/extension was calculated using V aw and V ae in Eq. 2. Sign changes were determined by comparing the vector cross product to the corresponding normal vector on the plane. Positive angles correspond to shoulder abduction, shoulder flexion, shoulder external rotation, and elbow flexion. These joint angle calculations were performed offline using a MATLAB function to get joint angle trajectories over time, and the resulting motion profiles were differentiated to get joint angular velocities. Table 1 shows the accuracy of the joint angle calculations reported here, compared to goniometer measurements. An iGaging 7 ′′ goniometer (Anytime Inc., Los Angeles, CA, USA) with a precision of 0.1°and resolution of 0.05°was used to enforce the measured (ground truth) readings presented in Table 1. Each of the eight joints (two angles per joint, indicating minimum and maximum angles) were measured independently in each trial. Sensor positions were recorded and filtered using PiMgr (see Data Capture) in 2 s recordings. Three repetitions were captured for each trial. Joint angles were calculated from the time-series position data using our model to yield 2-s long joint angle postures, which were then averaged across the 2 s to get a joint angle for each repetition of each trial. Mean and SD for the calculated angles shown in Table 1 are computed across repetitions. The difference between the model-estimated joint angles and the ground truth measured by goniometers for minimum and maximum angles is 8.8°and 8.6°, respectively. Error normalized to ground truth measurements reveals a mean of 36% error for the minimum measurement and 10% error for the maximum measurement, in large part because the minimum angle is smaller in magnitude than the maximum. Since reported results are in the velocity domain a large part of this error is negated in differentiation.

Subjects
Subjects were recruited, with written and informed consent, based on the criteria specified in the approved IRB. Healthy subjects with no history of right or left upper limb injury or weakness and no cognitive or motor impairments were allowed in the study after signing a consent form and filling out a basic medical questionnaire. Ten subjects of age 18-25 years were recruited (mean 20); of which 4 were female and 10 were self-reported right-hand dominant.

Experiment Procedure
Upon arrival subjects were fitted with hook-and-loop harnesses for the position sensors. Straps were adjusted so the sensors were held firmly in place without impeding the subject's motion. Basic range of motion exercises were performed to ensure that all straps and wires were settled and that the sensors remained in the proper locations. Subjects first executed 18 rapid training tasks, during which they were instructed to execute each task as quickly as they could successfully be completed. Three repetitions of each task were completed before performing the next task in the same fixed order for every subject. Eleven testing tasks were then performed in a fixed order at the subject's natural pace with three repetitions each. Each session lasted approximately 90 min with a short break offered between task phases. Figure 1B and Table 2 show each of the tasks performed in the study. These tasks were selected based on Barreca et al. (2005) and Foti and Koketsu (2013) as a cross-section of what is performed during ADL while requiring coordinated motion of both arms and being executable in the testing environment. These tasks were grouped into four categories to describe their type of motion: symmetric in-phase, symmetric out-of-phase, asymmetric, and coupled. Symmetric in-phase motions involve mirror symmetry between the two arms. This symmetry is typically about the midline but can be present in any direction. Symmetric out-of-phase motions have the same mirror symmetry as the in-phase category, except the motions of the left and right arms are offset in time. Asymmetric motions involve no symmetry between the two arms: each arm executes a different motion trajectory from the other such as washing a dish or using a fork and knife. Coupled motions involve both arms manipulating one object, such as when moving a box or tray, which results in a fixed relationship between the endpoints of each arm. Table 3 shows the number of each type of task present in the training and testing phase. Subjects were given instruction on what to do in each task with special care taken not to coach how to execute the motions. Subjects began each task with their hands in a flat resting position marked with red rectangles. Subjects were instructed to stop at the end of the task without returning to the rest position.

Minimum angle (°)
Maximum angle (°)   Tasks 1 and 22, the knife and fork task, involved the subject picking up a knife and fork, positioning them on a plate as if to cut food, and executing a cutting motion with the knife. Task 1 involved one cutting motion and task 22 involved three. Tasks 2 and 23, picking object off plate, involved the subject picking up a plate holding a small wooden object with their non-dominant hand and using their dominant hand to pick the object off the plate and place it on a target marked on the table. Tasks 3 and 4, scrub dish, involved the subject picking up a plate in their nondominant hand and using a sponge with their dominant hand to wipe the plate in a circular clockwise (task 3) or counterclockwise (task 4) motion for one complete cycle. Tasks 5 and 6, scrub table, consist of the subject reaching for a sponge and wiping it in a shallow upright (task 5) or upside down (task 6) V pattern along the surface of the table. The open box task, number 7, 8, and 24, involved the subject simultaneously opening the left and right flaps (task 7), the front and back flaps (task 8), or both sets of flaps sequentially (task 27) of a medium shipping box. Tasks 9 and 25, the fold clothes task, consist of identical motions between the training and testing set. The subject grasped a large piece of cloth along the outer edges and folded the left and right thirds over the center.

Joints
Task 10, drink from cup, consisted of the subject grasping a weighted cup placed in front of them with both hands and raising it up to their mouth as if to drink. Tasks 11-14, place cup on shelf, involved the subject picking up a weighted cup with both hands and placing it on one of four locations on a small set of shelves (top left, top right, bottom left, and bottom right sections of shelf for tasks 11-14, respectively). Tasks 15-18, pick cup off shelf, has the same sequence as 11-14 except the subject moves the cup from the shelf to the table.
Tasks 19-21, manipulate tray, involved the subject picking up a tray while tilting it to the left (task 19), right (task 21), or not tilting (task 20). Tasks 26 and 27, imagined steer wheel, involved the subject pantomiming a steering motion. Subjects were instructed to pretend to turn a steering wheel in a handover-hand fashion in either the clockwise direction (task 26) or the counterclockwise direction (task 27) for three complete cycles. Tasks 28 and 29, imagined ladder climb, involved the subject pantomiming climbing up (task 28) and down (task 29) a ladder while seated.
The measured velocity trajectories were windowed using a threshold of 5% maximum repetition velocity to identify task start/end. Windowed training tasks were then averaged across repetitions for each subject to produce 18 windowed, filtered, averaged joint angular velocity profiles for synergy derivation. Testing task data were converted to joint angular velocity as above and were windowed to task onset without averaging across repetitions. Testing tasks were only windowed to task start to ensure that task duration always exceeded synergy duration. Table 4 documents all symbols used in this section for reference. The measured velocity profiles for the training tasks were formatted into a matrix for each subject during data processing with J rows, T max columns, and N pages. For the training tasks, the study used J = 8 joints, T max = number of samples in the longest windowed training task, and n training tasks with N = 18 total.

Synergy Derivation
Tasks were padded with zeros on the end to reach T max length. T max varies for each subject and ranged from 400 to 608 samples. In order to extract spatiotemporal synergies using PCA, the 8 × T max dimensional velocity matrix V n for a given task is manipulated into a row vector, V n , with J·T max columns where each instant of time is represented as an additional set of J joints. Each task's row vector is concatenated into an N × J·T max task matrix, V: Principal component analysis operates on this matrix as though it were a J·T max -joint system observed at N instants of time, effectively treating the time-series velocity profiles as a single Frontiers in Bioengineering and Biotechnology | www.frontiersin.org February 2017 | Volume 5 | Article 2 instant of a J·T max -joint system. The task matrix is resolved into three-component matrices using SVD: In this form, the first m rows of S correspond to the first m principal components, or synergies, where m < N: U, Σ, and S are computed from V using the SVD function in MATLAB. The original V matrix can be approximated asṼ by isolating the first M columns of U, M × M elements of Σ, and M rows of S. The matrix U M diag{λ 1 , λ 2 ,. . . λ M } is now denoted as the weight matrix for the nth task and first M synergies for S M : It follows from Eq. 3 that the time-series components ofṼ over time can be expressed as a weighted sum of M principal components:ṽ Since the singular values found in Σ are related to the spread of data along each principle axis, i.e., variance in that direction, an index known as fraction of sum-squared variance can be calculated from the diagonal elements of Σ. This index describes the fraction of total variance accounted for by the first m synergies and is useful as an indicator of how many principal components are needed to adequately represent the data.
An index threshold of 95% variance is used to determine how many synergies, M, are required to represent the training task data.
Principal component analysis assumes stationary input variables, but the input data contain time-varying joint angular velocities. Here, "time-varying" refers to a motion consisting of a sequence of postures that change with time. However, the statistical properties of the synergies are quite stationary. PCA also assumes independent and identically distributed variables, but the input data that contain bilateral arm postures are not strictly independent (as there are biomechanical constraints that lead to joint correlations) as is the case with many real-world variables. PCA is a non-parametric method, i.e., it does not require any prior knowledge. Although this makes the application of this method simple, the method itself assumes linearity, which could be a weakness in many applications. We have compared the performance of PCA with other non-linear methods in Patel et al. (2015a) and found that PCA outperformed other methods. Using this exploratory analysis has previously led us to anatomically informing and meaningful synergies as principal components. These synergies could represent 100 postural movements with as low as six synergies with accuracy greater than 90% (Vinjamuri et al., 2010). These synergies also showed the effect on visual and tactile feedback in reaching and grasping movements (Patel et al., 2015b).

Reconstruction
Before performing reconstruction, the synergies and testing data were downsampled from 240 to 60 Hz due to the relatively long duration of synergies and testing tasks leading to excessive computation time. The testing task matrix, R, was reformatted from a J × T d × L matrix for J joints, T d samples in the downsampled task, and L = 11 tasks into a J·T d row, L column 2D matrix with each column defined as a separate task: Each task shorter than the longest task was padded with zeros to equal the same length. The objective of reconstruction is closely analogous to finding a representation of R in a new basis B. This is accomplished as an optimization problem to find the elements of column vector C which most closely satisfies: In this case, B will be made up of temporal offsets and dilations of the synergies defined by S m . B is formed by first transposing the downsampled version of S m . For notation, let S(m) be a J·T s element column vector containing synergy m and let [0] be a null column vector of J elements long. The first T d − T max columns of B are: where m = 1 is the synergy number, D = 0 identifies an undilated synergy (1, 2, or 3 indicating the first, second, or third synergy dilation), and T max is the length of the undilated synergy. Next, the synergy is dilated by interpolating S(m) to be some proportion of the difference between testing task length and synergy length. In this study, we dilate each synergy to be longer than the original synergy by 25, 50, and 75% of the difference in task and synergy sample length. Each dilation, S D (m), is used to form the next T d − T m,D columns of B, where T m,D is the number of samples in the Dth dilations of the mth synergy. This is done by forming each column as a sequential temporal offset of S D (m) as done above. This process is repeated for each dilation of each synergy, resulting in a B matrix, which contains every temporal offset of every dilation examined of the first M synergies: l 1 -norm minimization was used to sparsely select values for C which satisfy the following optimization problem: where ||·|| 1 and ||.|| 2 represent the l 1 -and l 2 -norms, respectively, and λ is a regulation parameter. Since the columns of B represent each synergy in different temporal offsets and dilations, the elements of C serve as recruitment weights for a synergy at a particular instant of time. The reconstructed task profiles,R, can be generated by multiplying B·C: The error between the measured and reconstructed profiles across all joints J = 8 is computed for each task, l using

Extraction of Synergies Using PCA
Spatiotemporal kinematic synergies were computed for each subject using PCA. Figure 2 shows the squared variance for each synergy along with the fraction of the sum of squared variance averaged across subjects. Derived synergies were 2.14 ± 0.29 s long, ranging from 1.67 to 2.53 s. The first synergy accounts for 57.98 ± 6.4% of total variance, while the first six synergies account for 91.05 ± 1.7% and the first eight account for 94.82 ± 0.85%.

FIGURE 2 | Fraction of variance accounted for by each synergy (bars)
and total from 1 to n synergies (line). The first synergy accounts for 57.98 ± 6.35% of variance, the first six synergies account for 91.05 ± 1.69%, and the first eight synergies account for 94.82 ± 0.85%. Dotted line shows 0.95 threshold.
The first three of each subject's synergies were compared to each other using Pearson's correlation coefficient. Synergies were put into the column form discussed in the Section "Reconstruction" and compared. Figure 3A shows the Pearson's r 2 averaged across subject comparisons for each combination of the first three synergies, leading to 90 unique comparisons for each synergy pair. Statistically significant differences were found in each of the three groups by one-way ANOVA tables, α = 0.05. All comparisons between synergy 1, 2, and 3 were statistically significant (p ≪ 0.005 for all). Tukey post hoc tests were computed to determine specific differences. As expected, the correlation between synergy 1 and synergy 1 was greater than the correlation between synergies 1 and 2 and between synergies 1 and 3. Correlations between synergies 1 and 2 and between synergies 1 and 3 were not found to be different from each other. All r 2 values are statistically different for the synergy 2 and 3 comparisons. A closer examination of the synergy correlations was conducted using Pearson's correlation coefficient, r. Figure 3B is a color-scale grid displaying the comparison between synergy one, two, and three between each pair of subjects. Synergy 1 appears positively correlated across all subjects except subject 6, who appears to have a strong negative correlation. Synergy 2 appears more mixed: there exist some positive/negative correlations that contribute to the overall mean r 2 of 0.3593 ± 0.2652, although subjects 1 and 10 have statistically insignificant mean r values at α = 0.05 of 0.0504 ± 0.265 (p = 0.5841) and 0.0594 ± 0.4582 (p = 0.7074), respectively. Figure 4 shows the synergy velocity profiles for the first eight synergies of subject 6. Each row corresponds to a joint labeled with a letter indicating the side ("R" for right or "L" for left), the joint ("S" for shoulder or "E" for elbow), and the rotation ("A" for abduction, "F" for flexion, "I" for internal rotation). Since the synergy is unscaled, the y axis is a unit-less velocity amplitude, while the x axis is time in seconds. For subject 6, the first synergy Frontiers in Bioengineering and Biotechnology | www.frontiersin.org February 2017 | Volume 5 | Article 2

FIGURE 3 | Correlation analysis of the first three synergies. (A)
Pearson coefficient of determination, r 2 , averaged across the 45 unique combinations between the specified synergies of each of 10 subjects. Statistical differences found using one-way ANOVA tables with α = 0.05 with Tukey post hoc tests. Comparison of synergy 1 to synergy 1 and synergy 1 to synergies 2 and 3, among all pairings of synergy 2, and among all pairings of synergy 3 yielded significant differences. (B) Correlation coefficient, r, between each subject for synergies 1, 2, and 3. Only unique pairings of subjects are shown using the upper triangle matrices. Synergy 1 appears highly positively correlated except for subject 6, who is highly negatively correlated.

FIGURE 4 | Angular velocity profile for first eight synergies of subject 6.
Vertical axis is unitless velocity since synergy is unscaled. Rows are labeled with letters indicating the side ("R" for right or "L" for left), the joint ("S" for shoulder or "E" for elbow), and the rotation ("A" for abduction in positive direction, "F" for flexion in positive direction, "I" for internal rotation in positive direction) for each DoF. Synergy 1 and 2 involve flexion and internal rotation of shoulder along with extension of elbow, implying a reaching motion, whereas synergy 3 involves extensions at the shoulder and flexion at the elbow.
consists of bilateral shoulder flexion, abduction, and internal rotation, which suggests a forward reaching motion. The elbows both show a small initial flexion, perhaps as the subject raises their arms from the rest posture, followed by an extension motion as they complete the reach. Synergy 2 behaves similarly, except shoulder abduction is delayed relative to synergy 1 and shoulder flexion and internal rotation are executed for a longer period of time.
As discussed above, synergy 3 is relatively uncorrelated among subjects. For subject 6, synergy 3 involves slight shoulder flexion followed by extension and adduction with an external rotation. The elbows appear to flex slightly, extend, then flex again.
In order to better visualize the other synergies for subject 6, each profile was integrated, multiplied by a gain, and added to the average starting joint angles of a particular subject. The gain was chosen such that the resulting angular profile remains within natural range of motion throughout the whole path. Figure 5 shows a FIGURE 5 | Posture visualization for first eight synergies of subject 6 (columns) over six normalized time instances (rows). Position at T = 0% is the subject's position averaged across the first 50 samples of all tasks. Synergies were integrated up to each time point and multiplied by a gain such that normal joint range of motion is not violated. Mirror symmetry between left/right arms can be seen in the first three synergies whereas synergies 4-8 have asymmetric motions. Subject's non-dominant hand tends to go to a single position and hold steady while the dominant hand appears to move continuously.
virtual mannequin posed at the resulting postures for six normalized time points. Subject 6, with a negatively correlated synergy 1, has a distinct reach-and-grasp motion. The downward dip of both hands observed at t = 0.4, or nearly halfway through the motion, could roughly correspond to the end of the reach phase and the beginning of object grasp and manipulation involving picking up an object. Note that subject 6's first synergy is strongly negatively correlated with the rest of the subjects, so synergy 1 for most of the subjects involves outward extension of the arms, similar to synergy 3 in Figure 4. Synergy 2 was similar, except it involved a reach-and-grasp motion instead of reach and manipulate/pick up. The shoulder motions observed in the synergy 3 profiles clearly result in an overall bilateral extension movement, with the arms spanning outward behind the back. Synergy 4 appears similar to picking up and holding a box at chest level.
The higher order synergies shown in Figure 5 appear to show a level of handedness, with the left hand tending to go to a certain position and holding while the right hand moved in various profiles through the duration of each synergy. Figure 6 shows the Pearson coefficient of determination, r 2 , averaged across subjects comparing joints 1-4 (right arm) to joints 5-8 (left arm) of all eight synergies. A one-way ANOVA (α = 0.05) was performed with Tukey post hoc tests to establish significant differences (ANOVA p ≪ 0.005). Synergies 1 and 2 were each FIGURE 6 | Pearson's coefficient of determination comparing right arm joints to left arm joints of the each synergy. One-way ANOVA and Tukey post hoc show that synergy 1 and synergy 2 each had significantly higher coefficients of determination than synergies 4-8, while the coefficient for synergy 3 was significantly higher than synergy 4 and 8. Statistical difference found in first three synergies implies asymmetric motion between left and right arms in higher order synergies.
Frontiers in Bioengineering and Biotechnology | www.frontiersin.org February 2017 | Volume 5 | Article 2 significantly more correlated between left/right than Synergies 4-8, and Synergy 3 was more correlated between left/right than Synergies 4 and 8.

Task Reconstruction
Reconstruction was carried out using up to eight synergies given the PCA results above and the 8-DoF nature of the bilateral arm model presented. Tasks were reconstructed with and without dilations. Figure 7 shows the mean reconstruction error calculated using Eq. 7 and averaged across degrees of freedom, tasks, and subjects using no dilations (blue) and dilations at 25, 50, and 75% of task/synergy length difference. Synergies 1-4 show a nearly linear decrease in normalized error in both cases; each subsequent synergy begins showing less of an improvement. Recruiting the first six synergies yielded a normalized reconstruction error of 0.1757 ± 0.0347 and 0.104 ± 0.0161 for no dilations and three dilations, respectively. Ultimately the error reduces down to 0.062 ± 0.0098 by synergy 8.
FIGURE 7 | Normalized reconstruction error when recruiting from synergies 1-8 with and without dilations. The reconstruction error for each synergy was averaged over degrees of freedom, subjects, and tasks. Dilated synergies were longer than undilated synergies by 25, 50, and 75% of the difference between minimum reconstruction task length and the synergy length. SDs are across subjects and tasks. Figure 8 shows examples of the reconstruction progression for left shoulder abduction using dilations for several tasks, demonstrating a clear progression from two available synergies (blue dotted line) to 8 synergies (red line). Figure 8A shows that the best-performing task that was repetition 2 of task 27, pretend steering wheel counterclockwise for subject 8. Reconstruction error went from 74.2% using only the first synergy to 2.22% using the first eight. Figure 8B shows repetition 2 of task 24, open box, for subject 5. Reconstruction error for this task was 27.9% using the first synergy and 2.41% using the first eight. Figure 8C shows repetition 1 of task 22, knife and fork, for subject 10. Reconstruction error was 43.8% using the first synergy and 7.32% using the first eight.
The optimization algorithm appeared to handle cyclic profiles relatively well, whereas the initial and ending phases of the tasks would often be less accurate. Figures 8A,C show this quality within the first 1 s of the tasks. Reconstruction was able to capture the overall shape of the task shown in Figure 8B, but some of the finer, irregular motions were not captured.
Tables 5 and 6 show the integrated recruitment gain for subject 4 averaged across tasks during normalized time bins for each synergy without and with dilations, respectively. Synergies and their dilations were not allowed to be recruited beyond the time at which the last synergy and task sample would coincide. Time bins which would include "missing" recruitment gains were therefore omitted. Significant differences found using one-way ANOVA with α = 0.00089 (0.05 over 56 comparisons) between synergy recruitments with and without dilations are indicated in Table 5 as bolded and underlined with an asterisk. For subject 4, synergies 3 and 5 in time bin 2 were significantly more recruited in the no dilations reconstruction than in the dilations reconstruction (p ≈ 0 for both). No other synergies had significantly different recruitments. Table 6 shows significant differences in recruitment between dilations of each synergy in each time bin, found using one-way ANOVA tables with α = 0.00125 (0.05 over 40 comparisons) and Tukey post hoc tests. Differences are bolded and marked with an asterisk. The undilated synergy 1 was significantly more recruited than all dilations in time bin 1 (p ≈ 0). The undilated synergy 8 was significantly more recruited than dilation 2 and 3 in time bin 1 (p ≪ 0.00125). The first dilation of synergy 7 was significantly more recruited than the un-dilated synergy in time bin 3 (p ≈ 0). The undilated synergy 8 was significantly more recruited than the    second dilation in time bin 3 (p ≪ 0.00125), although the absolute value of their recruitments likely would not be different. A similar analysis was performed after pooling all subjects. Significant differences in recruitment of an undilated synergy including and excluding dilations were found across tasks and subjects using ANOVA with α = 0.0013 (α = 0.05 for 40 comparisons). Synergy 1 was significantly less recruited in time bins 1 and 2 when dilations were included in reconstruction (p ≈ 0). Synergy 8 was significantly more recruited in time bin 1 with dilations compared to without dilations (p ≪ 0.0013). Differences in recruitment among dilations of each synergy in each time bin were found across tasks and subjects using ANOVA tables with α = 0.00156 (α = 0.05 for 32 comparisons). The undilated synergy 1 in time bin 1 is significantly more recruited than the dilated versions (p ≈ 0). The undilated and first dilation of synergy 5 are more recruited than dilation 3, but this difference was not statistically significant (p = 0.0034). Dilation 4 of synergy 7 was recruited significantly more than the undilated, first, and second dilation in time bin 1 (p ≪ 0.001). The undilated synergy 8 was statistically different from the first and second dilation in time bin 1 (p = 0.0015).

DISCUSSION
Previous work in implementing postural synergies in robotic control systems (see Introduction) aims to create a simplified control scheme, which can control high-dimensional systems with a reduced number of control inputs or actuators. The present study expands on this work by evaluating time-varying spatiotemporal synergies defined for both shoulder and elbow joints. Furthermore, we believe a key question for robotic systems (such as prosthetics and other assistive devices) meant to perform ADL, is how capable such systems are of reproducing ADL-like motion. Most experiments evaluate performance via endpoint variables such as task success rate, task completion time, or endpoint accuracy. This work attempts to evaluate such a system based on its ability to replicate ADL-like kinematics measured from reach-and-grasp experiments. This representation gives a more complete representation of the state of the system through the entire motion being performed. We present the optimal kinematic performance of time-varying synergies by sparsely optimizing the recruitment of synergies in time and amplitude and also introduce temporal dilations as discussed later.
Our results indicate that although the first three synergies can account for up to 80% of variance in training data, replicating the actual kinematics of ADL tasks could require the use of higher order synergies. Using two, three, or even four synergies as in previous literature could still be placing an absolute optimal kinematic accuracy limit of over 30% error in angular velocity, as shown in Figure 7. In the context of robotic control and assistive robotic interfaces, our results are reminiscent of the tradeoff between the optimal performance and complexity of control. An autonomous system may be able to incorporate higher order synergies in order to finely coordinating two arms, whereas an assistive system meant to be used by a person would quickly become too complex to use. These results may be improved on by experimenting with training task composition, dividing joints into different subgroups, or several other approaches which were not examined here.
Beyond simplifying control schemes, synergies have also been proposed as an avenue to make robotic systems exhibit more humanoid behaviors with biomimetic control. Robots designed to interact and cooperate with humans are deemed more esthetically pleasing when exhibiting anthropomorphic motion and are believed to be safer since human-like motions are more readily predictable than robotic ones (Duffy, 2003;Hegel et al., 2008). Control schemes aimed at anthropomorphizing robotic motions through the use of synergies derived from reach-and-grasp tasks have been carried out by Liarokapis et al. (2015) and reviewed by Santello et al. (2016). Whereas this work leveraged postural synergies, our results indicate that time-varying synergies could be used to produce more accurate replications of human motion. It is possible that the reconstruction method applied here could be performed on non-anthropomorphic trajectories to generate a near-natural approximation of a synthetictrajectory.
The recent trend of mechanically embedding synergies into robotic systems may also be extended to time-varying principal components. The SoftHand, for example, uses a compliant structure actuated by a single principal component to achieve various grasps (Catalano et al., 2014). Implementing the complex motions described in Figure 4 may be achieved by software-generated temporal postural synergies coupled with a mechanical design which together introduce non-linearities. These time-varying synergies in bilateral arm movements can be easily integrated in upper limb prosthetics and rehabilitation devices to render natural movements.
Our results could also have implications outside of strict robotic control. The symmetric nature of the first two synergies between the left and right arms led to an examination of the similarities between the left and right joint angular velocities of synergies 1-8 for each subject. These results indicate that lower order synergies tended to have similar profiles between the left and right arms, whereas higher order synergies were dissimilar. The lower order synergies seem to provide broad reaching motions, while higher order synergies appear to fine-tune the motion to suit particular manipulations required by a task.
The asymmetries in higher order synergies could contain information related to the handedness of the subject. Observationally, the postural visualization of the higher order synergies in Figure 3 appears to show the left arm maintaining stable positions and motions, while the right arm follows a more complex path. This result is in line with previous work done on bimanual upper limb control, and a more detailed kinematic analysis of bilateral kinematic synergies may also be conducted to explore this link. Sainburg and Kalakanis (2000) and Bagesteiro and Sainburg (2002) found significant differences in hand path curvature and torque efficiency between dominant and non-dominant hands; the same group found advantages in the dominant arm for speed and directional control, whereas the non-dominant arm was specialized for accurate position control through increased limb impedance (Wang and Sainburg, 2007). Recently, Yokoi et al. (2014) conducted bimanual lever tasks in which the non-dominant arm operated in an artificial force field. They were able to demonstrate that the non-dominant limb was more adaptable to the force field when trained with a specific motion direction of the dominant arm, but that altering the dominant-side target location led to a decrease in performance of the non-dominant arm. In this paper, a model using movement primitives to approximate motor learning of the non-dominant arm relative to dominant arm dynamics was able to replicate this effect, supporting the idea that bilateral synergies may be used by the central nervous system to produce motion.
The synergies are hypothesized to be abstract representations of motion encoded in the sensorimotor system of the brain, which can be combined to form the complex motions required to execute ADL tasks. The actual generation and execution of complex upper limb movements is a neural process that still has much to be discovered. The original drive or goal which the sensorimotor system acts on has long been believed to be the limbic system (Brooks, 1986). The abstract goals produced by the limbic system are processed by the association cortex that then produces a plan for motion. This information is passed to the projection system, composed of the sensorimotor cortex, cerebellum, and basal ganglia, which converts the "trajectory" to a series of commands which can be passed to the spinal cord. This process is dictated by sensory feedback, developmental history, and a sense of the dynamics of the body (Brooks, 1983). The internal dynamic model has been hypothesized to be composed of both feedforward and feedback components, aimed at minimizing error while still allowing for quick responses. The existence of either of these models or the combination thereof is still much debated, with evidence pointed either way (Wolpert et al., 1998). As mentioned earlier, the studies conducted by Viviani and Terzuolo (1980) imply the projection system's ability to dilate motion building blocks in time. Single impulses at specific interneuronal sites result in coactivation of multiple muscles, as shown by Tresch et al. (1999) and Saltiel et al. (2001), and each muscle generates torque around multiple degrees of freedom in the hand and arm (Santello et al., 2013). This torque is influenced by numerous anatomical parameters which are unique to each individual, as is the resulting motion of the limb (Buchanan et al., 2004). Visual, proprioceptive, haptic, and numerous other sensory signals are then returned to the brain and integrated into motion planning.
Task reconstruction as performed in this study attempts to model the neural process described above computationally. The measured training task can be viewed as the motion plan output by the associative cortex. The synergies derived from the rapid tasks serve as the "history" programmed into the projective system which dictates how the brain converts planned motion to spinal cord commands. This would be analogous to a feedforward component of cerebellar processing discussed by Wolpert et al. (1998): given a certain output from the brain, there is expected to be some dynamic response observed at the arms that can take on an arbitrary pattern across joints. Dilations of the synergies are also provided which attempts to mimic the time-scaling ability of the basal ganglia and cerebellum. l 1 -norm minimization serves the role of the projection system by using synergies and their dilations to convert from a desired motion to a time-series of gains, which can be viewed as the time-series of signals being sent to the spinal cord. The neuromuscular and physiological construction of the limbs is approximated by multiplying the synergies by this time-series of gains, producing the final motion.
The time-series of gains generated by l 1 -norm minimization could give some insights into how useful temporal dilations of synergies could be in motion production. Table 5 shows that recruitment of synergy 1 is reduced when including dilations in reconstruction for the first two time bins. Synergy 8 had a larger recruitment in the first time bin with dilations compared to without. Table 6 shows that the only statistical difference between recruitment of dilations is in time bin 1 for synergies 1, 7, and 8. The undilated synergy 1 was significantly more recruited than the dilations, dilation 3 was the most recruited for synergy 7, and the undilated synergy was the most recruited for synergy 8. The most significant effect of including dilations in the recruitment process appears to be the reduction in recruitment of synergy 1 during task onset. Recruitment of the dilations of synergy 1 also appears to be relatively large compared to recruitment of the other synergy dilations.
The present experiment could be refined to further study dominant/non-dominant performance in bimanual ADL. The subject pool can be expanded to include left-and right-handed individuals, and tasks can be modified to selectively engage the left, right, or bilateral arms. The relative influence of dominant/non-dominant hemispheres on ipsilateral limb motion could be examined with a more focused data analysis by examining whether left or right side joints are better reconstructed for left-and right-handed subjects. Devising a way to "coax" reconstruction to replicate the structural asymmetry and increased influence of the dominant over the non-dominant hemisphere as opposed to vice versa (Snyder et al., 1995;Amunts et al., 1996;Ziemann and Hallett, 2001;Hayashi et al., 2008) could allow us to study the strength of this influence. The upper limb kinematic data recorded in this study may also be augmented with EMG and/or EEG signals.

CONCLUSION
This paper presents a study of spatiotemporal kinematic synergies in bilateral arm movements during ADL tasks. These tasks were selected as a representation of the fundamental categories of bilateral motion (symmetric, asymmetric, coupled), with further work aimed at tailoring the task list to target specific aspects of sensorimotor processing. Derived synergies accounted for up to 94.82 ± 0.85% of variance and were demonstrated to reconstruct ADL tasks to within 6.2 ± 0.98% using the l 1 -norm minimization algorithm. The concept of temporal synergy dilations were incorporated to replicate movement processing found in the basal ganglia and cerebellum. Recruitment patterns were examined throughout task duration and found that the dilated version of a synergy was used equally as much as the undilated version of the same synergy in most time bins, with the beginning of the motion having the most difference. Potential uses of time-varying synergies to anthropomorphize robotic motion and inform mechanical construction were discussed. Interesting features of the synergies in the context of handedness corroborating others' results were discussed. We believe that synergies will be instrumental in building next-generation biomimetic prosthetics and orthotics in the near future.

ETHICS STATEMENT
The present study was conducted under IRB Approved Protocol # 2014-026/2015-022 at the Stevens Institute of Technology. All the subjects gave written informed consent in accordance with the Declaration of Helsinki.

AUTHOR CONTRIBUTIONS
MB designed the experiment, ran instrument calibrations, designed the kinematic model and data processing scripts (with inputs from IF and KP), executed the experiment, and wrote the manuscript. VP contributed to the data processing scripts, assisted with executing the experiment, and provided valuable feedback in writing the manuscript. RV set up the hypothesis and direction/scope of the study, provided invaluable feedback in experiment design, task selection, and data processing, coordinated subject recruitment, and provided key revisions and editing on the manuscript.