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.

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₁ and S₂ were placed at the lateral head of the clavicle on the subject’s right and left shoulder, respectively, while S₃ was placed on the subject’s right side near the middle of the rib cage on the midaxillary line. S₄ and S₆ were placed on the lateral side of the subject’s elbows over the joint’s center of rotation. S₅ and S₇ 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: Si=xiyizi Vshould=S1−S2 Vside=S1−S3 VRe=S1−S4,VLe=S2−S6 VRw=S4−S5,VLw=S6−S7

Abduction/adduction was found by first projecting the elbow vector, V _ae, onto the coronal plane: (1)Vae,proj=Vae−(Vae⋅nc)nc∥nc∥ 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 (2)θsa=cos−1Vshould⋅Vae,proj∥Vshould∥∥Vae,proj∥

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 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.

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.

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.

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 non-dominant 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.

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 hand-over-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. Vn=v1n(1)…v1n(Tmax)⋮⋱⋮vJn(1)…vJn(Tmax)=v1n(1)…v1n(Tmax)⋮⋱⋮v8n(1)…v8n(Tmax)

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: V=V¯1V¯2⋮V¯n=v11(1)⋯v81(1)⋯v11(Tmax)⋯v81Tmaxv12(1)⋯v82(1)⋯v12(Tmax)⋯v82Tmax⋮v1N(1)⋯v8N(1)⋯v1N(Tmax)⋯v8NTmax

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 instant of a J⋅T_max-joint system. The task matrix is resolved into three-component matrices using SVD: V=U∑S where U is an N × N square matrix with orthogonal columns, Σ is an N × J⋅T_max diagonal matrix, and S is a J⋅T_max × J⋅T_max square matrix with orthogonal rows. The diagonal elements of Σ correspond to the singular values, λ, of V. ∑=diag{λ1,λ2, …, λN}

In this form, the first m rows of S correspond to the first m principal components, or synergies, where m < N: S=s11(1)⋯sJ1(1)⋯s11(Tmax)⋯sJ1(Tmax)s12(1)⋯sJ2(1)⋯s12(Tmax)⋯sJ2(Tmax)⋮s1M(1)⋯sJM(1)⋯s1M(Tmax)⋯sJM(Tmax)⋮s1N(1)⋯sJN(1)⋯s1N(Tmax)⋯sJN(Tmax)

U, Σ, and S are computed from V using the SVD function in MATLAB. The original V matrix can be approximated as  V˜ by isolating the first M columns of U, M × M elements of Σ, and M rows of S. The matrix U_Mdiag{λ₁, λ₂,… λ_M} is now denoted as the weight matrix for the nth task and first M synergies for S_M: (3)V≈ V˜=UMdiagλ1,λ2,…,λMSM=WMnSM

It follows from Eq. 3 that the time-series components of  V˜ over time can be expressed as a weighted sum of M principal components: (4)v˜jn(t)=∑m=1Mwmnsjm(t)

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. (5)λ12+λ22+...+λM2λ12+λ22+...+λN2

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).

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: R=r11(1)⋮rj1(1)⋮r11(Td)⋮rj1(Td),r12(1)⋮rj2(1)⋮r12(Td)⋮rj2(Td),⋯r1L(1)⋮rjL(1)⋮r1L(Td)⋮rjL(Td)

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: (6)R=BC

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: BDm=S(m)[0][0][0]S(m)[0][0][0]…[0]⋮⋮⋮[0][0]S(m) 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: B=B01…B31…B0M…B3M l₁-norm minimization was used to sparsely select values for C which satisfy the following optimization problem: Minimize ∥C∥1+1λ∥BC−R∥22 where ||⋅||₁ and ||.||₂ represent the l₁- and l₂-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: R≈ R˜=BC

The error between the measured and reconstructed profiles across all joints J = 8 is computed for each task, l using (7)el=∑  j=1J∑  t=0Tl(rjjl(t)−rjl(t))2∑  j=1J∑  t=0Tlrjl(t)2

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%.

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² 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² 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² 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 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 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², 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 significantly more correlated between left/right than Synergies 4–8, and Synergy 3 was more correlated between left/right than Synergies 4 and 8.

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 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).