The Nonlinear Relationship Between Sensory and Motor Primitives During Reaching Movements

The stabilizing role of sensory feedback in relation to movement dynamics remains to be poorly understood in realistic three-dimensional movements of limbs. The objective of this experimental and computational study was to classify the contribution of sensory feedback from muscle spindles to the control of assistive and resistive limb dynamics during human pointing movements. We integrated a human upper-limb musculoskeletal model with a model of Ia primary afferent discharge to analyze motion and muscle activation patterns during reaching movements in virtual reality (VR). The reaching target locations in VR were selected to define movements with varying roles of gravity and interaction torques that created diverse dynamical contexts. Nine healthy human subjects performed the VR task by pointing to the reaching targets with visual feedback of their arm location. Motion capture and electromyography (EMG) were recorded, and joint torques and Ia primary afferent discharge were estimated using the integrated model. The experimental and simulated data were analyzed with hierarchal clustering (Gritsenko et al., 2016). The clustering analysis of EMG and predicted proprioceptive signals showed a divergent relationship between muscle activation and sensory feedback patterns. Even though the Ia models had nonlinear dynamical components, their output was still most related to the anatomical grouping of muscles and less so to the dynamical contexts of each movement, reflected in muscle activations. Altogether, these results suggest that sensory feedback is nonlinearly related to muscle activation profiles and that it may contribute information necessary for coupling between proximal and distal muscle groups.


Introduction
Movements are the product of interactions between neural control signals and the musculoskeletal dynamics that depend on limb anatomy (Yakovenko, 2011;Gritsenko et al., 2016). The control of this complex system is accomplished by the hierarchically organized central nervous system (CNS). Musculoskeletal morphology has traditionally been viewed as an additional complexity with redundant characteristics that the CNS is required to solve (Bernstein, 1967). The supraspinal systems are thought to solve this complexity through the dynamical cortical mechanisms (Churchland et al., 2012) that generate a variety of movements through combinations of relatively few neural signals (Tresch et al., 2006;Cheung et al., 2009;d'Avella and Lacquaniti, 2013). The precision of such control is particularly important to compensate for passive moments around joints that arise from the inertial properties of the limb and external forces, such as gravity. These passive moments play two main roles during movement, resistive or assistive. An example of the resistive role of passive moments is during simultaneous extension or flexion of both shoulder and elbow joints in a horizontal plane. During such movement, interaction torques at shoulder and elbow joints are in the opposite direction to net torques, that correspond to the direction of the desired movement. These resistive passive moments need to be overcome by additional muscle contraction, and they have been shown to be predictively counteracted by the CNS (Almeida et al., 1995;Gribble and Ostry, 1999;Koshland et al., 2000;Pigeon et al., 2003;Debicki and Gribble, 2005;Kurtzer et al., 2008;Gritsenko et al., 2011). An example of the assistive passive moments is akin to a whiplash action between adjacent joints. Skilled baseball players increase the efficiency of their throws by using the whiplash action of interaction torques at the elbow that are initiated by torques applied at the shoulder (Hirashima et al., 2007). To allow for the assistive passive moments to have effect, the magnitude of active moments due to muscle contractions need to be reduced. However, the assistive action of passive moments also needs to be kept in check, so that it does not cause instability. This could be accomplished by the viscoelastic properties of muscles and spinal reflexes, which can resist perturbations and ensure stable control with minimal supraspinal input (Asatryan and Feldman, 1965;Brown and Loeb, 2000;Yakovenko et al., 2004;Prochazka and Yakovenko, 2007;Valero-Cuevas et al., 2015). Other evidence points toward the stabilizing action of sensory feedback through inter-joint coordination (Sainburg et al., 1993;. Altogether, evidence suggests that the supraspinal structures need to embed the complex mechanical dynamics of the limb in order to appropriately shape the control signals and incorporate sensory feedback (Kluzik et al., 2008;Kurtzer et al., 2008;Wagner and Smith, 2008).
The objective of this study is to classify the contribution of Ia primary afferent discharge to the control of assistive and resistive limb dynamics during human pointing movements. To achieve this goal, we selected a small set of whole arm movements that were accompanied by either assistive or resistive passive joint torques as defined by an inverse dynamic model of the arm. Then, we designed a virtual reality task that guided healthy human subjects through these movements. The recorded muscle activity and limb kinematics were then used to compare experimental motor commands to sensory feedback computed using models.

Experimental design and human subjects
We have recruited 9 healthy adults (5 males, 4 females, 24.3 ± 1.8 years old, 76.3 ± 14.5 kg) to perform upper extremity reaching movements in a virtual reality (VR) environment. All procedures were approved by the West Virginia University Institutional Review Board (Protocol #1309092800). The movement task was created with a wearable VR helmet (Oculus Rift) integrated with motion capture. Subjects moved to virtual targets on cue with the visual feedback of their arm position (Fig. 1A). To minimize the inter-subject differences in motion, the locations of all targets were calculated mathematically based on subject's segment lengths and shoulder and elbow joint angles, so that in movement #1 subjects' shoulder extended while elbow flexed, in Movement #2 both shoulder and elbow flexed, and in Movement #3 shoulder extended while elbow flexed to a different end position. These movements defined diverse dynamical contexts where the movement was largely passive (#1), interaction torques were resistive with increasing gravitational load (#2), and interaction torques were assistive with decreasing gravitational load (#3, Fig. 1B). During movement, we recorded kinematics of shoulder, elbow, and wrist joints and electromyography (EMG) of 12 muscles that span those joints. Muscles recorded were anterior and posterior deltoids (ADelt and PDelt respectively), pectoralis major (Pec), teres major (TMaj), biceps brachii long and short heads (BicL and BicS respectively), triceps brachii lateral and long heads (TriLa and TriLo respectively), brachioradialis (BR), extensor carpi radialis (ECR), flexor carpi radialis (FCR), and flexor carpi ulnaris (FCU). These muscle abbreviations are used consistently throughout the manuscript and figures. Motion capture data were recorded at 480 Hz using the Impulse system (PhaseSpace) and EMG was recorded at 2000 Hz with MA400-28 (MotionLab Systems). The onset and offset of each movement was defined using a velocity threshold method using wrist and elbow LED markers.

Models
To calculate joint torques, an inverse dynamic model of the subject's arm was constructed in Simulink (MathWorks). The model comprised five degrees of freedom (DOFs) including shoulder (flexion/extension, abduction/adduction, pronation/supination), elbow (flexion/extension), and wrist (flexion/extension) and three segments approximating inertial properties of the arm, forearm, and hand ( Fig. 1B). Inertia of the segments was approximated with a cylinder with 3-cm radius and the length equal to that of the corresponding segment. The masses and centers of mass for each segment were determined by their anthropometric ratios to the subjects' segment lengths and weights (Winter, 2009). Angular kinematics averaged per movement and per subject was used in the subject-specific inverse model to calculate joint torques (Fig. 2). These computed torques are proportional to the sum of all moments generated by muscles spanning the joints, so these torques are referred to as muscle torques in the rest of the manuscript. The muscle torques were then differentiated to obtain rotatums, whose sign is indicative of the directionality of the active forces produced to make a given movement. Rotatums were then separated into two separate positive signals based on the sign of the rotatum signals for comparison with EMG profiles (Fig. 2). Subjects moved primarily in a sagittal plane, therefore shoulder abduction/adduction and pronation/supination DOFs were excluded from all following analyses.  To estimate the sensory contribution from muscle spindles during movement, we have used a published model of Ia primary afferent discharge (Prochazka, 1999). The model was derived from cat afferent recordings and validated for human afferent recording by Malik et al. (Malik et al., 2016). The offset of the model was adjusted from 80 to 10 impulses per second to reflect human microneurography data as in Malik et al. (2016). The spindle model relates afferent firing rate (Ia) to the rate of change of muscle length (v), muscle length (l), and normalized EMG (a) as follows: Muscle length was normalized to vary between 0 and 1, which corresponded to minimal and maximal possible muscle lengths respectively, calculated as described below. Velocity was in the units of rest length per second; a given muscle rest length was defined as described below. EMG was high-pass filtered at 5 Hz, rectified, low-pass filtered at 20 Hz, and normalized to the maximum across all movements.
Muscle lengths were calculated using a modified musculoskeletal model of the human arm in OpenSim (Saul et al., 2015;Gritsenko et al., 2016) (Fig. 1C). The model was scaled for each participant's segment lengths by adjusting proportionally the muscle origin and insertion points and muscle paths to individual segment dimensions. We used the model to estimate muscle length changes during movement from joint kinematics. The recorded motion capture data were used to calculate averaged joint kinematics for each movement and each subject as described above. The DOFs of the OpenSim model were set using the angular kinematics and the corresponding muscle lengths were calculated for each sample. Muscle rest length was generalized as the average muscle length within the physiological range of motion. These data and EMG were used in equation (1).

Analysis
We used regression analysis to explore the relationships between muscle activations and Ia primary afferent discharge. All signals were aligned on movement start and normalized to the duration of movement using onset and offset events defined from motion capture as described above ( Fig. 2 and 3). We created normalized EMG and Ia profiles for each muscle and each movement per subject by averaging signals across 20 repetitions. Then, we calculated correlation coefficients (r) between all pairs of EMG and Ia profiles. These r values were then converted into the heterogeneous variance explained (HVE) as follows: This transformed the large positive r values that were characteristic of agonistic relationships into short distances close to 0 and the large negative r values corresponding to antagonistic relationships into long distances close to 2. Lastly, we used hierarchical clustering across participants on unbiased HVE distance matrix using the linkage function with un-weighted average distance method (Gritsenko et al., 2016).  Consistent muscles were determined based on belonging within the same clusters across the three movement types. Cluster labels were first matched across subjects and subsequently across movements based on the majority of signal types inside each cluster. Then the total number of signals that fall within clusters with the same label across all movements was calculated for hierarchical clustering of EMG and Ia and separately for hierarchical clustering of EMG and rotatums. The cluster assignments were obtained for 2 to 12 clusters for EMG/Ia clustering and 2 to 9 clusters for EMG/rotatum clustering. The upper limit of the total number of clusters was set at half the total number of signals, e.

Statistics
Consistency of hierarchal clustering was evaluated using analysis of variance with repeated measures (rANOVA). Five rANOVAs were applied in MATLAB to the outcomes of hierarchical clustering for individual subjects (Table 1). rANOVA1, rANOVA2, and rANOVA3 will refer to the statistical analysis of cluster assignments of EMG and Ia signals, while rANOVA4 and rANOVA5 will refer to the statistical analysis of cluster assignments of EMG and rotatum signals. All rANOVAs compared cluster assignments of recorded data with noise estimated with a Bootstrap procedure (Efron and Tibshirani, 1993). The Bootstrap procedure included random permutation of cluster assignments 900 times, 100 times per subject. The rANOVAs were designed to test five null hypotheses. The null hypotheses of rANOVA1 and rANOVA4 were that the number of agonist signals in each cluster is equal to that occurring by chance (rANOVA1 on EMG/Ia and rANOVA3 on EMG/rotatum). EMG and Ia signals were considered agonistic if they originated from the same muscles. EMG and rotatum signals were considered agonistic if they span the same joint and acted in the same direction, e.g. BicL and elbow flexion rotatum. The other three null hypotheses were that the number of experimental clusters containing only signals of one type, e.g. only EMG, or only Ia, or only rotatum, is equal to that occurring by chance (rANOVA2, rANOVA3, and rANOVA5 respectively). All rANOVAs contained three factors. The first was a between-factor Group defining subject cluster assignments vs. randomized cluster assignments. The randomized cluster assignments were averaged across 100 permutations to match the number of subjects for a balanced ANOVA design. The second was a within-factor Movement defining data belong to the three movements. The third was a within-factor Cluster defining the level of hierarchy, i.e. 2-cluster breakdown, or 3-cluster breakdown, and so on until the number of clusters equaled half of the number of signals. Main and interaction effects were investigated using MATLAB ranova function, and post-hoc analysis between levels of the within factors was done using MATLAB multcompare function.

Results
Subjects performed movements in a highly-consistent manner. The angular excursions of each joint were very similar across subjects for the three movements (Fig. 2, top row). This resulted in highly consistent muscle torques and their components, the flexion and extension rotatums (see Methods, Fig.  2). Movement #1 was initiated with active extension shoulder torque accompanied by largely passive motion of elbow and wrist. This movement was stopped through primarily active elbow flexion and wrist extension torques. Movement #2 was initiated with active flexion torques at the shoulder and elbow and extension torque at the wrist. This movement was stopped through simultaneous active shoulder and elbow extension torques and wrist flexion torque. Movement #3 was initiated with active flexion torques at the elbow and wrist and accompanied by largely passive shoulder motion (Fig. 2, Shoulder). This movement was stopped through simultaneous active elbow extension torque and wrist flexion torque. The movements were accompanied by somewhat more variable muscle activation and Ia profiles (Fig. 3). The profiles of EMG were less variable during Movement #2 compared to the other movements. During Movement #2 shoulder was flexing, while during the other movements it was extending, thus the simulated profiles of Ia were markedly different for shoulder muscles during Movement #2 compared to the other movements. In contrast, the profiles of EMG during Movement #3 were largely similar across most muscles, which reflects co-contraction. The Ia profiles that most closely followed the pattern of co-contraction were from wrist flexor muscles (Fig. 3, FCR and FCU).
Overall, we found that EMG profiles were highly correlated with each other and with Ia profiles in all three recorded movements (Fig 4A). The clustering was very robust across subjects, so that most signals ended up in the same clusters at multiple levels of the hierarchy (Fig 5A). Surprisingly, EMG and Ia profiles from agonist muscles did not group together into the same clusters. Statistical analysis has shown that there were fewer agonistic EMG and Ia profiles grouped into common clusters than expected by chance ( Table 2, significant effects of Group factor). This was true for movements #2 and #3, but not #1 (Table 2, post-hoc tests per movement) at all hierarchal levels except for 2-and 3-cluster breakdown. We also found that EMG and Ia profiles tended to cluster into separate groups. EMG clustering changed between movements and EMG clusters included muscles spanning different joints that are not traditional agonists (Fig. 5, red lines). The number of homogeneous EMG clusters was higher than expected by chance ( Fig 6A; Table 3, main effects and insignificant interactions). This was significant for movements 1 and 3 (Table 3, post-hoc tests) at all hierarchal levels except for 3-and 4-cluster breakdown.        The clustering of Ia profiles tended to change less across movements. Furthermore, Ia profiles from muscles that are traditional agonists tended to cluster together. For example, FCU and FCR were always closely correlated in all movements. This was also true for PDelt and TMaj, and for BicL, BicS, and Br ( Fig 5A, blue lines). The number of homogeneous Ia clusters was higher than expected by chance, but to a lesser degree than the number of homogeneous EMG clusters (Fig. 6B, Table 3). This was significant only for the movement #3 (Table 4, post-hoc tests) and only at 4-, 6-and 7-cluster breakdown. During the same movements, fewer strong correlations were observed between EMG profiles and rotatums, although the latter were correlated among themselves (Fig. 4B). Statistical analysis has shown that there were fewer agonistic EMG and rotatum signals grouped into common clusters than expected by chance ( Table 5, significant effects of Group factor). This was true for movements #2 and #3, but not #1 (Table 5, post-hoc tests per movement) at all hierarchal levels. Hierarchal clustering further showed that rotatums were mostly grouped into two clusters (Fig. 5B, green lines). One cluster consisted of Shoulder flexion, Elbow flexion, and Wrist extension, and the other cluster consisted of Shoulder extension, Elbow extension, and Wrist flexion. The only exception was movement #3, in which shoulder flexion rotatum was largely zero (Fig. 2) and both shoulder rotatums changed their cluster allocations. Statistical analysis has shown that the number of homogeneous rotatum clusters was higher than expected by chance ( Fig 6C; Table 6, main effects and insignificant interactions). This was significant for movements #2 and #3 (Table 6, post-hoc tests) at all hierarchal levels except for 2-cluster breakdown.

Discussion
The objective of this study is to quantify the contribution of Ia primary afferent discharge to the control of assistive and resistive limb dynamics during human pointing movements. The main finding of this study is that clustering of Ia profiles is distinct from that of EMG profiles during movements with distinct dynamic conditions. These results suggest that the primary afferent signals are providing information that is of a different modality than the outgoing activity of motor neurons. Indeed, the ensemble firing rate of motor neurons has been closely linked to the forces produced by the muscles. In contrast, the primary afferent firing rate is most related to the muscle length and its rate of change (Prochazka, 1999). There is a known monosynaptic relationship between primary afferents and motor neurons innervating the same and synergistic muscles that underlies stretch reflexes, which compensate for perturbations. This anatomical arrangement would result in common signals that would be observed between profiles of the activity of primary afferents and the profiles of the activity of homologous motor neurons. However, our results have shown that there is very little common variance between the Ia profiles and EMG profiles from agonistic muscles. This indicates that the simple reflex relationship is not a major contributor to common muscle activations, or synergies, observed during reaching movements. This supports prior studies showing that short-latency stretch reflexes do not scale together with mechanical interactions between arm joints when perturbations are applied prior to movement (Kurtzer et al., 2009;Weiler et al., 2015). Instead, these studies have shown that long-latency reflexes, that presumably involve descending neural signals, do scale with mechanical interactions between arm joints in a goal-dependent manner. Our results add the evidence that during movement execution primary afferent feedback needs to be processed by the CNS prior to being incorporated into the ongoing motor commands.
The movements selected for this study represent different dynamical contexts, where the movement #1 was largely passive, movement #2 was against gravity and accompanied by resistive interaction torques between all joints, and movement #3 was with gravity and accompanied by assistive interaction torques between shoulder and elbow. To compensate for these dynamic conditions, the forces generated by the muscles, as reflected in active muscle torques, were different in the three movements, even when the joint motion resulting from these torques was the same (Fig. 2, wrist). Despite the varied dynamical conditions, we observed consistent clustering of rotatum signals from different joints across all movements. This occurred even though the movement trajectories were unconstrained. This suggests that there is a robust inter-joint coordination that is independent of dynamic condition. This is consistent with known invariant characteristics of motion kinematics and dynamics (Soechting and Lacquaniti, 1981;Hollerbach and Flash, 1982). However, unlike joint torques, the common signals in muscle activation were not consistent with the observed invariances in inter-joint coordination. Indeed, EMG clustering changed between movements and did not appear to include consistent agonists or muscles spanning the same joints within a given cluster ( Fig. 4 and 5). Instead, the common signals in EMG may be more likely related to the different dynamic conditions experienced by the limb during these movements. For example, in the movement #3 the clustering of antagonistic muscles with assistive interaction torques may serve to increase joint stiffness, which would help stabilize the movement. In movement #2, the clustering of predominantly antigravity flexors may serve to overcome resistive interaction torques and increasing gravity load experienced by the limb during that movement (Fig. 5). Thus, common activation of muscle groups may reflect the neural compensation for complex limb dynamics (Gritsenko et al., 2011).
Most muscles span at least two joints and even more DOFs, because the anatomy of most joints support motion along several rotational DOFs. For example, shoulder can move along three rotational DOFs represented by functional motions of flexion/extension, abduction/adduction, and internal/external rotation. Although each DOF is independent by definition, muscle contractions would naturally couple the DOFs that belong to the joints that those muscles span. Thus, we would expect coupling between shoulder and elbow DOFs, because multiple muscles span both of those joints. In a prior computational study with the musculoskeletal model of the human arm, we have observed such coupling between muscle lengths across all joint excursions even without muscle contractions (Gritsenko et al., 2016). We termed it mechanical coupling, as it arose purely from the anatomical arrangement of muscles on the skeleton. Interestingly, the Ia clusters observed in the current study are similar to those identified through mechanical coupling, i.e. common signals in Ia primary afferent profiles occurred among agonist muscles that span the same joints ( Fig. 5 red lines). Notably, mechanical coupling separated muscles into proximal and distal groups, where lengths of muscles that spanned shoulder and elbow joints did not correlate with the lengths of muscles that spanned wrist and finger joints. This separation between proximal and distal muscles is also evident in the Ia primary afferent clusters of the present study. In contrast, we have observed clustering of rotatum signals across all three modeled joints, the wrist rotatums were not clustering separately from the shoulder and elbow rotatums (Fig. 5 green lines). Similarly, the clustering of EMG in different movements was not divided into proximal and distal groups. Altogether, this suggests that one of the functions of the neural motor command for the observed pointing movements may be to couple proximal and distal muscle groups defined by anatomy using the proprioceptive information about the state of each of these groups.

Funding
This study is supported by a scholarship from NIH/NIGMS U54GM104942 (RLH), Ruby Fellowship (MTB), and NIH CoBRE P20GM109098 (SY, VG). The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH.