Uncontrolled Manifold Reference Feedback Control of Multi-Joint Robot Arms

The brain must coordinate with redundant bodies to perform motion tasks. The aim of the present study is to propose a novel control model that predicts the characteristics of human joint coordination at a behavioral level. To evaluate the joint coordination, an uncontrolled manifold (UCM) analysis that focuses on the trial-to-trial variance of joints has been proposed. The UCM is a nonlinear manifold associated with redundant kinematics. In this study, we directly applied the notion of the UCM to our proposed control model called the “UCM reference feedback control.” To simplify the problem, the present study considered how the redundant joints were controlled to regulate a given target hand position. We considered a conventional method that pre-determined a unique target joint trajectory by inverse kinematics or any other optimization method. In contrast, our proposed control method generates a UCM as a control target at each time step. The target UCM is a subspace of joint angles whose variability does not affect the hand position. The joint combination in the target UCM is then selected so as to minimize the cost function, which consisted of the joint torque and torque change. To examine whether the proposed method could reproduce human-like joint coordination, we conducted simulation and measurement experiments. In the simulation experiments, a three-link arm with a shoulder, elbow, and wrist regulates a one-dimensional target of a hand through proposed method. In the measurement experiments, subjects performed a one-dimensional target-tracking task. The kinematics, dynamics, and joint coordination were quantitatively compared with the simulation data of the proposed method. As a result, the UCM reference feedback control could quantitatively reproduce the difference of the mean value for the end hand position between the initial postures, the peaks of the bell-shape tangential hand velocity, the sum of the squared torque, the mean value for the torque change, the variance components, and the index of synergy as well as the human subjects. We concluded that UCM reference feedback control can reproduce human-like joint coordination. The inference for motor control of the human central nervous system based on the proposed method was discussed.

Introduction

There are redundant relationships underlying the motor control of the human body. For instance, a simple voluntary reaching movement involves a redundant hand trajectory, redundant joint angles, and redundant muscles. We focused on the redundancy between the hand position as a single target and the greater number of degrees of freedom (DOFs) of the joints.

To approach the above redundancy problems, Bernstein (1967) pointed out that “the human controls his/her redundant DOF of body by using joint coordination.” Specifically, the joint coordination is a control strategy in which redundant elements (motor elements) are varied without affecting the variable that must be controlled to achieve the task (performance variable). In other words, the human central nervous system (CNS) facilitates families of solutions equally able to solve the task as described by Latash's Principle of Abundance (Latash, 2012). Therefore, joint coordination can flexibly stabilize the performance variables. To quantify the joint coordination of human movements, uncontrolled manifold (UCM) analysis has been proposed (Scholz and Schöner, 1999). UCM analysis divides the trial-to-trial variance of redundant motor elements into two orthogonal components: a UCM component that does not affect the performance variable and an ORT (orthogonal) component that directly affects the performance variable. If a particular performance variable is controlled by the coordination of motor elements, the UCM component is greater than the ORT component. UCM analysis has been used to investigate many types of human tasks including the reaching movement (Domkin et al., 2002, 2005; Tseng et al., 2002, 2003; Yang et al., 2007), and has suggested a coordinated structure for voluntary human movement (Latash et al., 2002, 2007; Latash, 2010).

Within the framework of the UCM, the task can be achieved by only controlling the ORT component. This concept in which only the task-relevant elements can be controlled was also mentioned as the “minimal intervention principle” by Todorov and Jordan (2002a). They introduced a traditional optimal feedback control framework to achieve the principle (Todorov and Jordan, 2002b). In the framework of optimal feedback control, however, there is neither a control model explicitly using the UCM nor an evaluation of the joint coordination. Moreover, to apply the framework to a nonlinear system, which is the nature of the human body, a complicated algorithm is needed (Todorov et al., 2005; Li and Todorov, 2007; Liu and Todorov, 2009).

In this paper, we propose a novel control method of human redundant joints, which achieves the “minimal intervention principle” using the UCM directly. The main purpose of the present study is to confirm that our proposed method can generate arm movements with human-like joint coordination at a behavioral level. Specifically, we consider a one-dimensional target-tracking task that is more redundant than the traditional two-point reaching movement. In this task, a target of the hand in the lateral direction is provided. As a conventional method, we considered the method that pre-determined a unique target joint trajectory by inverse kinematics or any other optimization method. The controller tracked the pre-determined unique target trajectory. Standard inverse-kinematics approaches cannot implement motor variability in the joint space without violating the task space constraint, which are given by the tracking of an end-effector pose. Therefore, the conventional method does not model appropriately the variability of joint angles. In contrast, our proposed method generated the target UCM from the given target position of the hand, and selected the optimal joint combination in the target UCM at each time step.

The proposed method does not generate optimal and unique target joint trajectory. Instead, the optimal joint combination in the target UCM is determined each time step. During arm movements, an arm posture is disturbed by a noise of motor command and an external perturbation. By selecting the optimal joint combination in the target UCM, the variability of joint due to the noise and perturbation is minimally corrected. Therefore, the variability of the joint along the UCM is allowed and results in joint coordination. Moreover, the controller can smoothly and effortlessly respond to the perturbation by allowing variability in the target UCM. We refer to the proposed method as “UCM reference feedback control.” To check whether our proposed method can generate human-like joint coordination, we conducted measurement experiments of human subjects performing a one-dimensional target-tracking task. We quantitatively compared the performances of the arm movements and joint coordination quantified using the UCM analysis in simulation experiments with those from the measurement experiments.

Materials and Methods

UCM Reference Feedback Control

One-Dimensional Target-Tracking Task

In this study, the common controlled object in the simulation and measurement experiments is a three-link arm consisting of an upper arm, forearm, and hand. The common task is a tracking task of the one-dimensional target of the hand in the horizontal plane (Figure 1A). Thus, the joint angles and hand position correspond to the motor elements and the performance variable in the UCM concept. Only the lateral position of the target is considered so that the relationship between the target and motor elements is considered redundant. The task coordinates (X- and Y-axes) and joint angles (θ_s, θ_e, and θ_w) are defined in Figure 1A.

FIGURE 1

Figure 1. Concept of UCM reference feedback control. (A) Schematic diagram of the task coordination and joint angles in the horizontal plane. For the one-dimensional target-tracking task, the hand is moved from the start position to the end position tracking the target hand position along the X-axis while variability along the Y-axis is allowed. (B) Schema of UCM reference feedback control. UCM_t denotes the target UCM at time t. The configuration of the joint angles is disturbed at time t, and deviates from the target UCM at time t + 1. Through the UCM reference feedback control method, the disturbed joint angles converge to the target UCM at time t + 2. (C) Block diagram of UCM reference feedback control. The UCM reference feedback controller is given the target hand trajectory x^d, and calculates the unique reference joint angles θ^d* within the target UCM so as to minimize the cost function. Then, the inverse dynamics of the arm generates the input torque τ using the current joint condition θ and $\stackrel{\cdot }{\theta }$ and reference joint angles θ^d*. The upper (pseudo) loop represents the visual feedback loop while the lower loop represents the somatosensory feedback.

The one-dimensional target of the hand is provided for in the medial–lateral direction (X-axis), and is calculated according to the minimum jerk criterion of the start and end velocities equaling zero (Flash and Hogan, 1985). Meanwhile, the hand is allowed to take any position in the anterior–posterior direction (Y-axis). The origin is the center of gyration of the shoulder. The arm movements are considered for two initial positions: a far position (θ_s = π/3 rad, θ_e = π/3, and θ_w = π/6 rad) and a near position (θ_s = π/4 rad, θ_e = π/2 rad, and θ_w = π/6 rad). These conditions indicate that initial hand position is either far or near from the shoulder. From these initial positions, the hand moves 0.3 m in 5 s. This movement is slower than an ordinary reaching movement because the present study considers accurate target-tracking movements of the human.

Brief Description of UCM Reference Feedback Control

UCM is a manifold that can be defined when the motor elements have redundant DOF(s) with respect to the performance variable. In the tracking task, any value of a joint angle within the UCM results in the same value of the hand position. Thus, the combination of the joint angles within the UCM accomplishes the task. To simplify the explanation of the UCM, we give an example of the UCM in the target-tracking task with a two-link arm consisting of a shoulder and elbow, and a one-dimensional target trajectory of the hand as shown in Figure 1B. The horizontal and vertical axes denote the shoulder and elbow angles (motor elements), and UCM|_t denotes the target UCM at time t. At each time step, our proposed method first generates and refers the target UCM, then, determines optimal joint combination in the target UCM, finally, input torque is calculated by inverse dynamics.

UCM reference feedback control first generates the target UCM (UCM|_t) for a given target position x^d(t). Next, a unique joint combination in the target UCM (the target joint combination) is selected so as to minimize a cost function for the current joint combination. Here we consider a cost function consisting of torque and the torque change terms that encourages small and smooth joint torques so as to represent a human's smooth and effortless movement. For instance, when the joint combination at time t (θ_t) is disturbed and the joint combination θ_t+1 turns away from the target (Figure 1B), the optimal target joint combination in the UCM at t + 2 (θ_t+2) is selected so as to minimize the torque and torque change.

Finally, the input torque is calculated to track the optimal target joint combination using a current state feedback and inverse dynamics of the arm. Since the optimal joint combination depends on the current state, the target joint combinations are different trial by trial. Therefore, the UCM reference feedback control can achieve the task while permitting variance of the joints in the UCM. Such a mechanism that allows task-irrelevant variability of joints generates the joint coordination. The block diagram of our proposed method is shown in Figure 1C. The upper (pseudo) loop and lower loop correspond to the visual feedback and somatosensory feedback, respectively. To simplify the problem, we assume that there is no time delay or noise in the somatosensory feedback system.

Algorithm for UCM Reference Feedback Control

An algorithm for the implementation of the UCM reference feedback control is explained. We now consider the one-dimensional target-tracking task in a horizontal plane with a three-link arm (Figure 1A). The equation of motion for the three-link arm is:

$\begin{array}{l}M\left(\theta \right)\ddot{\theta }+V\left(\dot{\theta },\theta \right)+D\dot{\theta }=\tau \left(1\right)\end{array}$

where $\mathbf{M}\ddot{\theta }$ is an inertial term, V is a centripetal and Coriolis term and $\mathbf{D}\dot{\theta }$ is a viscous term. The state variables are the joint angles and angular velocities.

A relationship between the joint angles (motor elements) and the hand position (performance variable) is represented by the following kinematics equation (UCM):

$\begin{array}{l}{x}^d=l_{1}cos({\theta }_s)+l_{2}cos({\theta }_s+{\theta }_e)+l_{3}cos({\theta }_s+{\theta }_e+{\theta }_w)\left(2\right)\end{array}$

where θ_s, θ_e, and θ_w are subject to Equation (2) when generating the UCM. A given target hand position is represented by x^d. Hence, the update of x^d changes the UCM in a step by step manner. The target joint combination subject to the UCM as represented by Equation (2) is determined so as to minimize a cost function consisting of the joint torque and torque change terms. Thus, we consider the cost function to select the target joint combination:

$\begin{array}{l}E={\tau }^T\text{W}\tau +{\frac{d\tau }{dt}}^T\text{Q}\frac{d\tau }{dt}\left(3\right)\end{array}$

where W = diag[W_s, W_e, W_w] and Q = diag[Q_s, Q_e, Q_w] are weighted matrices.

Next, we explain how to solve an optimization problem that minimizes Equation (3) with the constraint condition of Equations (1) and (2). The input torque and input torque change for the cost function (3) are calculated from the motion equation of the arm (1) as follows. To represent the torque and torque change by θ_s, θ_e, and θ_w in the UCM Equation (2) as variables, Equation (1) is approximated using the Euler method. To represent a relationship between the input torque (τ_t) at t and a joint combination θ_t+2 at t + 2, Equation (1) is approximated as:

$\begin{array}{l}{\theta }_{t+1}={\theta }_t+\Delta t{\dot{\theta }}_t\left(4\right)\end{array}$

$\begin{array}{l}{\dot{\theta }}_{t+1}={\dot{\theta }}_t+\Delta t{M}^{-1}({\tau }_t-D{\dot{\theta }}_t-V)\left(5\right)\end{array}$

$\begin{array}{l}{\theta }_{t+2}={\theta }_{t+1}+\Delta t{\dot{\theta }}_{t+1}\left(6\right)\end{array}$

where Δt indicates the sampling duration in the Euler method. Substituting Equations (4) and (5) into Equation (6), a relationship between the target joint combination two steps latter ${\theta }_{t+2}^d$ and the current input torque is obtained:

$\begin{array}{l}{\theta }_{t+2}^d={\theta }_t+2\Delta t{\dot{\theta }}_t+{\left(\Delta t\right)}^2{M}^{-1}\left({\tau }_t-D{\dot{\theta }}_t-V\right)\left(7\right)\end{array}$

To simplify a calculation, the joint combinations θ^†_t+2 and θ^‡_t+2, which are generated by setting τ_t = 0 and τ_t= τ_t−1, are calculated by using forward dynamics in Equation (7) according to:

$\begin{array}{l}{\theta }_{t+2}^{†}={\theta }_t+2\Delta t{\dot{\theta }}_t+{\left(\Delta t\right)}^2{M}^{-1}\left(0-D{\dot{\theta }}_t-V\right)\left(8\right)\end{array}$

$\begin{array}{l}{\theta }_{t+2}^{‡}={\theta }_t+2\Delta t{\dot{\theta }}_t+{\left(\Delta t\right)}^2{M}^{-1}\left({\tau }_{t-1}-D{\dot{\theta }}_t-V\right)\left(9\right)\end{array}$

Equations (8) and (9) are used for calculation of input torque and input torque change. Subtracting Equations (8) and (9) from Equation (7) and transforming the result, the input torque τ_t and input torque change $\dot{\tau }=({\tau }_t-{\tau }_{t-1})/\Delta t$ are obtained as:

$\begin{array}{l}{\tau }_t=M\frac{{\theta }_{t+2}^d-{\theta }_{t+2}^{†}}{{(\Delta t)}^2}\left(10\right)\end{array}$

$\begin{array}{l}\frac{{\tau }_t-{\tau }_{t-1}}{\Delta t}=M\frac{{\theta }_{t+2}^d-{\theta }_{t+2}^{‡}}{{(\Delta t)}^3}\left(11\right)\end{array}$

Therefore, the cost function in Equation (3) can be represented as a function of the target joint combination. By substituting Equation 2 for the UCM into the cost function, an optimal solution is obtained for the joint combination in the target UCM thereby guaranteeing the optimal joint combination that corresponds to the target hand position. To simplify the calculation, we use the UCM as linearized by the Jacobian:

$\begin{array}{l}{x}^d-x^{†}=J\left({\theta }^{†}\right)\left({\theta }_{t+2}^d-{\theta }^{†}\right)\left(12\right)\end{array}$

$\begin{array}{l}{x}^d-x^{‡}=J\left({\theta }^{‡}\right)\left({\theta }_{t+2}^d-{\theta }^{‡}\right)\left(13\right)\end{array}$

where J is the Jacobian matrix of Equation (2), and x^† and x^‡ are the hand positions corresponding to θ^† and θ^‡, respectively. Equations (12) and (13), which are linearized around θ^†_t+2 and θ^‡_t+2, are substituted into the expression for the torque in Equation (10) and the torque change in Equation (11), respectively. $\theta {}_s^d$ eliminated by substitution. Here, we assume that θ^†_t+2 and θ^‡_t+2 are near the target UCM two steps latter. To calculate the cost function (3), the other variables (i.e., the inertial term M, the target hand position x^d, θ^†_t+2, and θ^‡_t+2 obtained by the forward dynamics) are required. Combining them into a constant term A (details are shown in Appendix A), the cost function can be represented as a function of ${\theta }_e^d$ and ${\theta }_w^d$:

$\begin{array}{l}E=A_1{({\theta }_e^d)}^2+A_2{\theta }_e^d+A_3{\theta }_e^d{\theta }_w^d+A_4{\theta }_w^d+A_5{({\theta }_w^d)}^2+A_6\left(14\right)\end{array}$

To calculate the extreme values for ${\theta }_e^{d\ast }$ and ${\theta }_w^{d\ast }$, Equation (14) is partially differentiated with respect to ${\theta }_e^d$ and ${\theta }_w^d$. The value ${\theta }_s^{d\ast }$ is calculated by substituting ${\theta }_e^{d\ast }$ and ${\theta }_w^{d\ast }$ for which the cost function (3) takes a minimal value into Equation (2); therefore, that the target joint combination is in the target UCM. From the above calculations, the optimal solution can be analytically obtained, and the target joint combination can thus be calculated without any iteration. Then, substituting the target joint combination θ^d* into Equation (7) and solving the equation for the inverse dynamics, the actual input torque τ_t can be calculated as:

$\begin{array}{l}{\tau }_t=M\frac{{\theta }_{t+2}^{d*}-{\theta }_t-2\Delta t{\dot{\theta }}_t}{{(\Delta t)}^2}+D\dot{\theta }+V\left(15\right)\end{array}$

The above algorithm is implemented for the three-dimensional motor elements and the one-dimensional performance variable. In general, even if the motor elements are n-dimensional and the performance variable is m-dimensional (m < n), our proposed method can be directly applied. It should be noted that a more complex equational representation of the analytical solution is required for a higher dimensional system.

Simulation Experiment

A computer simulation of the UCM reference feedback control was performed for the one-dimensional target-tracking task. The physical parameters of a three-link arm were calculated using anthropometric data (Winter, 2004). The values given in Table 1 were calculated from the length of the body segments and the body weight of a typical subject in the measurement experiment. The value of the viscous matrix D is given in Table 2. The initial posture is given by the measurement data of the typical subject. The sampling duration was 0.0083 s (at 120 Hz, which corresponds to the sampling rate of the measurement experiments). The fourth-order Runge–Kutta method was used for the integration of the motion equation of the arm. Since the wrist torque should be smaller than shoulder and elbow torques in the three-link arm reaching, we specified W_s = 1/$I_1^2$, W_e = 1/$I_2^2$, and W_w = 1/$I_3^2$ as the weighted terms in the cost function in Equation (3). According to the minimum torque change criterion that can reproduce human reaching movements (Uno et al., 1989), the weight of the torque change term was a unit matrix (i.e., Q_s = Q_e = Q_w = 1).

TABLE 1

Table 1. Arm parameters.

TABLE 2

Table 2. Viscous matrix in the simulation.

To reproduce the variability of the human arm movement, noise was added to the three variables in Figure 1C. The added noise is considered reasonable for biological motor systems. Additionally, the simulation was run for 100 trials for both the far and near initial positions. The first of the three variables is the perceived target hand position x^d. We assumed uncertainty as caused by visual localization in visual perception:

$\begin{array}{l}{x}_{in}^d=x^d+e_x{v}^{d}N\left(0,1\right),\left(16\right)\end{array}$

where N(0, 1) is a normal distribution (mean 0, variance 1), x^d is the given target hand position, $x_{in}^d$ is the perceived target position used in UCM reference feedback control, v^d is the target velocity and e_x is the noise amplitude. We assumed that the perception of the moving target was more uncertain than that of a static target. Therefore, the visual perception noise was proportional to the target velocity. The second variable was the optimal joint angle θ^d*. Assuming that there was uncertainty in acquiring the target UCM, we used

$\begin{array}{l}{\theta }_{in}^{d*}={\theta }^{d*}+e_{\theta }N\left(0,1\right),\left(17\right)\end{array}$

where ${\theta }_{in}^{d\ast }$ is the target joint combination substituted into the inverse dynamics of the arm (Equation 15) and e_θ is the noise amplitude.

The third variable is the input torque τ_t. Harris and Wolpert (1998) suggested that the hand trajectory in a reaching movement varied owing to signal-dependent noise (SDN) in the motor command. In this study, we assumed SDN in joint torque:

$\begin{array}{l}{\tau }_{in}=\left(U+e_{\tau }N\left(0,1\right)\right){\tau }_t,\left(18\right)\end{array}$

where U is an identity matrix, τ_t is the joint torque generated by the inverse dynamics (Equation 15), τ_in is the actual joint torque inputted to the arm and e_τ is the noise amplitude. The noise amplitudes e_x, e_θ, and e_τ were heuristically determined so as to quantitatively fit the mean value data across all subjects, especially the total variance as shown in Figure S1 (Appendix B). The profile of the sum of the squared torque is shown in Figure 6. (specifically, e_x = 1.5 × 10⁻⁵, e_θ = 1.0 × 10⁻⁶, and e_τ = 0.22). We evaluated the initial and end positions, sum of the squared tracking error (i.e., the difference between the hand and target position), the tangential hand velocity, input joint torque and torque change in this simulation. Additionally, joint coordination was quantitatively evaluated through UCM analysis. We calculated the mean value data for the above kinematics and dynamics across 100 trials and compared them to mean value data across all subjects in the measurement experiments by one-sample t-test.

Measurement Experiment

Eight healthy right-handed male subjects participated in the measurement experiments. The experiments were approved by the Nagoya University Ethical Review Board. All subjects were provided with explanations regarding the experimental procedure and gave their written informed consent.

We used a 3D position measurement system (OPTOTRAK CERTUS, Northern Digital Inc.) to record the kinematics at 120 Hz. Infrared-ray markers with a diameter of 7 mm were placed on four anatomical landmarks of a subject's arm: the center of gyration of the shoulder, the elbow, the wrist, and the tip of the index finger (Figure 2). The index finger was fixed at an extended position so that the length of the hand (l₃) was sufficiently long. A plastic board that was easy to slide on a desk was placed under the hand. Subjects secured a head-mounted display (HMD) (HMZ-T1, SONY Inc.) to their heads so as to obtain their hand positions and the target hand position in the medial–lateral direction. The refresh rate of the screen of the HMD was 60 Hz.

FIGURE 2

Figure 2. Schematic diagram of the measurement experiment setup. Subjects perform a one-dimensional target-tracking task while sitting in a chair and wearing a seatbelt. A head-mounted display (HMD) is secured to the individual's head. Markers are placed on the arm segments identified in Figure 2. The one-dimensional hand position of the subject and the target hand position are shown by screen of the HMD.

Subjects performed the one-dimensional target-tracking task while sitting on a chair and wearing a seatbelt. The far position task and near position task were each performed for 100 trials. Thus, the subjects performed the one-dimensional target-tracking task for 200 trials. To confirm that there was no effect of the task order, four subjects performed the far position task first, and the other four subjects performed the near position task first. When the subjects felt fatigued, they took a short break arbitrarily. To adjust the initial positions for all trials, the subjects moved their elbow, wrist, and fingertip for the target of the initial posture with accuracy of ±2 mm using information from the HMD.

In the initial position control phase, the lengths of all links were calculated from the kinematics data, and the target of the initial posture was calculated with the kinematics equation so as to meet the task condition (far or near). After this phase, a countdown of 3 s was displayed on the HMD. The one-dimensional hand position at the movement end, target hand position and subject's hand position were then displayed as shown in Figure 2. Subjects performed the one-dimensional target-tracking task only using the information on the HMD. At the end of each trial, subjects were given feedback of the sum of the tracking error. They were instructed to track the target hand position so as to reduce the tracking error. In addition, they were asked to keep their elbow horizontal so as to avoid friction with the desk. To eliminate outlier or failure trials, we applied the following two screening processes. First, we eliminated the trials in which the value of the squared sum of the tracking error was larger than the mean value plus 2 standard deviations. Second, the trials in which the value of the mean sum of the squared torque change was larger than the mean value plus 1 standard deviation were rejected. Finally, 79.3 ± 13.3 trials were analyzed.

The position data were filtered with a second-order Butterworth low-pass filter with a 10-Hz cutoff frequency. We obtained the start and end positions, the tangential velocity of the hand and the joint torque from these position data. The tangential velocity of the hand was obtained by calculating the differential of the measured position and then filtered with a second-order Butterworth low-pass filter with a 10 Hz cutoff frequency. Using the same method, we obtained the joint angle, angular velocity, and angular acceleration from the measured position, and obtained the joint torque from the inverse dynamics of the arm in Equation (1) using the same viscous matrix as that in the simulation experiment (Table 2).

UCM Analysis

The joint coordination was evaluated through the following UCM analysis. The UCM is the nonlinear manifold associated with the redundant kinematics, and is locally linearized by the null space of the Jacobian J between the hand velocity in the X-direction and the joint velocity.

The variances parallel to the UCM and orthogonal to the UCM are called the UCM component (V_UCM) and ORT component (V_ORT), respectively. These variance components can efficiently be computed from the covariance matrix C of the observed joint angles (Yen and Chang, 2009):

$\begin{array}{l}{V}_{\text{U}\text{C}\text{M}}=\frac{\text{t}\text{r}\text{a}\text{c}\text{e}(\text{n}\text{u}\text{l}\text{l}{\left(J\right)}^{T}C\text{n}\text{u}\text{l}\text{l}(J))}{n-d}\left(19\right)\end{array}$

$\begin{array}{l}{V}_{\text{ORT}}=\frac{\text{trace}({(J{J}^T)}^{-1}JC{J}^T)}{d}\left(20\right)\end{array}$

$\begin{array}{l}{V}_{\text{T}\text{O}\text{T}}=\frac{\text{t}\text{r}\text{a}\text{c}\text{e}(C)}{n}\left(21\right)\end{array}$

where V_TOT is the total variance. These variance components are normalized by the dimensions of the joint combination n (n = 3), the dimensions of the hand position d (d = 1), and the dimensions of the null space n–d. To evaluate the degree of joint coordination, we define the index of synergy σ* (Verrel, 2010) as:

$\begin{array}{l}{\sigma }^{*}=log\left(\frac{\sigma +n∕d}{n∕(n-d)-\sigma }\right)\end{array}$

$\begin{array}{l}\sigma =\frac{V_{\text{U}\text{C}\text{M}}-V_{\text{O}\text{R}\text{T}}}{V_{\text{T}\text{O}\text{T}}}\left(22\right)\end{array}$

when V_UCM = V_ORT then the joints are not coordinated. Thus, σ* > log((n/d)/(n/(n–d))) = 0.69 means that the joints are varied parallel to the UCM, which implies a synergetic stabilization of the performance variable (Latash, 2010).

Results

We confirmed that there was no effect of task order on the performance and characteristics of the movements.

Kinematic and Dynamic Properties

Figure 3 shows the start and end arm postures. The upper (Figures 3A,B) and middle (Figures 3C,D) figures indicate the results of the simulation and the measurement experiments. In both experiments, hand positions at the movement end were more varied in the anterior–posterior direction (Y-direction), and the variance of joint angles did not affect task achievement. Hand positions in the near position task were nearer the trunk at the movement end than those in the far position task. Figure 3E shows the mean hand position at the movement end in the Y-direction for all experiments. Our proposed method generated quantitatively similar hand positions to those of all subjects. Statistically, a one-sample t-test between the simulation and measurement results showed no significant difference [the far position task: t₍₇₎ = −0.45, P = 0.67; the near position task: t₍₇₎ = −1.06, P = 0.33]. Moreover, a paired t-test demonstrated that the measured hand positions in the near position task were significantly nearer than those in the far position task [t₍₇₎ = 4.50, P = 0.0028 < 0.05], and the proposed method could generate same tendency (the far position task: 0.38 m; the near position task: 0.29 m).

FIGURE 3

Figure 3. Initial and end postures of the arm in the simulation and measurement experiments. The upper left (A,C) and upper right (B,D) graphs correspond to the far position task and near position task. The horizontal axis and vertical axis denote the medial–lateral direction and anterior–posterior directions. The blue lines denote the initial posture and the red lines denote the end postures. The origin is the center of gyration of the shoulder. (E) The mean value of the end hand position in the Y-direction for all experiments. The n. s. indicates not significant (one-sample t-test), and the asterisk denotes a significant difference (paired t-test, P < 0.05).

Figure 4 shows the sum of the squared tracking error. Our proposed method could more accurately track the target than the human subjects. A one-sample t-test between the simulation and measurement results showed significant difference between the tracking error of the simulation and measurement experiments [the far position task: t₍₇₎ = 9.39, P = 3.24 × 10⁻⁵ < 0.05; the near position task: t₍₇₎ = 4.85, P = 0.0019 < 0.05].

FIGURE 4

Figure 4. The sum of the squared tracking error. (A) The mean value for the tracking error for the simulation experiments. (B) The mean value for the tracking error for all subjects for the measurement experiments.

Figures 5A,B show the profiles of tangential velocity of the hand. The red and blue solid lines indicate the results of simulation and measurement experiments. The blue area denotes standard deviation across all subjects. The target hand trajectory in the one-dimensional target-tracking task had a bell-shaped and smooth velocity profile with a peak at the middle of the movement duration (2.5 s). In both the simulation and measurement experiments, the profiles of the tangential velocity of the hand were also bell-shaped. Figure 5C shows the peak of the tangential hand velocity for all experiments. A one-sample t-test between the simulation and measurement results indicated that our proposed method could generate a similar peak compared to the measurement experiments [the far position task: t₍₇₎ = 1.93, P = 0.09; the near position task: t₍₇₎ = 2.27, P = 0.06]. A paired t-test demonstrated that the measured peak of the tangential hand velocities was not significantly different between the far and near position tasks [t₍₇₎ = 1.10, P = 0.31]. However, in the measurement experiments the mean peak of the tangential hand velocity in the far position task (0.135 m/s) tended to be larger than that of the near position task (0.128 m/s). Our proposed method also showed the same tendency (the far position task: 0.124 m/s; the near position task: 0.119 m/s).

FIGURE 5

Figure 5. Tangential hand velocity for both the simulation and measurement experiments. The upper left (A) and right (B) graphs denote the far position and near position tasks. The red and blue lines indicate the mean value profiles for the tangential hand velocity across all trials for the simulation experiments, and across all subjects for the measurement experiments. The blue area denotes the standard deviation across all subjects. (C) Peak of tangential hand velocity of all experiments. The n. s. indicates not significant (one-sample t-test), and the asterisk denotes a significant difference (paired t-test, P < 0.05).

Figures 6A,B show the profiles of the sum of the squared torques for all joints. The red and blue solid lines indicate the results of the simulation and measurement experiments. The blue area denotes the standard deviation across all subjects. The sum of the squared torque was also bell-shaped while the profiles for the simulation and measurement experiments were similar. Figure 6C shows the peak of the sum of the squared joint torque for all experiments. A one-sample t-test between simulation and measurement results indicated that our proposed method could generate a similar peak to the measurement experiments [the far position task: t₍₇₎ = 0.84, P = 0.43; the near position task: t₍₇₎ = 0.38, P = 0.71]. A paired t-test demonstrated that the measured peak of the tangential hand velocities was not significantly different between the far and near position tasks [t₍₇₎ = 1.37, P = 0.21]. However, for the measurement experiments the mean peak of the tangential hand velocity in the far position task [0.200 (Nm)²] tended to be smaller than that for the near position task [0.246 (Nm)²]. Our proposed method also showed the same tendency [far: 0.174 (Nm)²; near: 0.239 (Nm)²]. Figure 6D shows the mean sum of the squared torque change for all experiments. Figure 6D shows the sum of the squared joint torque change of all experiments. A one-sample t-test between the simulation and measurement results indicated that our proposed method could generate a peak value similar to the measurement experiments [the far position task: t₍₇₎ = 1.30, P = 0.23; the near position task: t₍₇₎ = 0.80, P = 0.45].

FIGURE 6

Figure 6. Sum of the squared torques for all joints (shoulder, elbow, and wrist) for both the simulation and measurement experiments. The upper left (A) and right (B) graphs correspond to the far position and near position tasks. The red and blue lines indicate the mean value profiles of the sum of the squared torque across all trials for the simulation experiments, and across all subjects for the measurement experiments. The blue area denotes the standard deviation across all subjects. (C) Peak of the sum of the squared torque for all experiments. The n. s. indicates not significant (one-sample t-test). (D) The mean value of the sum of the squared torque change for all experiments.

Results of UCM Analysis

Figures 7A–D show the waveforms of the UCM and ORT components. The upper and middle figures show results for the UCM and ORT components, respectively. The red and blue solid lines indicate the results of simulation and measurement experiments. The blue area denotes the standard deviation across all subjects. In both the simulation and measurement experiments, the mean waveforms of the UCM components were larger than those of the ORT components throughout the duration of movement, which indicates that the variance of joint angles was more varied across the UCM. Moreover, the UCM and ORT components gradually increased from the movement initiation to the end and our proposed method could generate same tendency. Figures 7E,F show the mean values for the UCM and ORT components for all experiments. A one-sample t-test between the simulation and measurement results indicated that our proposed method could generate a similar mean value for the UCM and ORT components from the measurement experiments [the UCM component in the far position task: t₍₇₎ = −0.33, P = 0.75; in the near position task: t₍₇₎ = 0.54, P = 0.60; the ORT component in the far position task: t₍₇₎ = 0.28, P = 0.79; in the near position task: t₍₇₎ = 0.73, P = 0.49].

FIGURE 7

Figure 7. UCM and ORT components for both the simulation and measurement experiments. The upper left (A,C) and right (B,D) graphs correspond to the far position and near position tasks. The red and blue lines indicate the mean value profiles of the UCM and ORT components across all trials for the simulation experiments, and across all subjects for the measurement experiments. The blue area denotes the standard deviation across all subjects. (E) The mean value for the UCM components across the movement duration for all experiments. The n. s. indicates not significant (one-sample t-test). (F) The mean value for the ORT component.

A paired t-test demonstrated that the measured UCM and ORT components were not significantly different between the far and near position tasks [the UCM component: t₍₇₎ = −0.61, P = 0.56; the ORT component: t₍₇₎ = 1.66, P = 0.14]. However, for the measurement experiments the mean value for the UCM component in the far position task [0.0017 (rad)²] tended to be larger than that of the near position task [0.0014 (rad)²]. In contrast, the mean value for the ORT component in the far position task [1.52 × 10⁻⁴ (rad)²] tended to be smaller than that for the near position task [1.86 × 10⁻⁴ (rad)²]. Our proposed method also showed the same tendency, i.e., the UCM component in the far position task was 0.0018 (rad)²; the near position task was 0.0013 (rad)²; the ORT component in the far position task was 1.44 × 10⁻⁴ (rad)²; and the near position task was 1.69 × 10⁻⁴ (rad)².

Figures 8A,B show the profiles of index of synergy. The red and blue solid lines indicate the results of the simulation and measurement experiments. The blue area denotes the standard deviation across all subjects. The horizontal green line indicates the value of no coordination (V_UCM = V_ORT). In both the simulation and measurement experiments, the index of synergy was larger than the value of no coordination, which indicates that the joint angles were well-coordinated. These results indicate that the UCM reference feedback control could generate a high index of synergy as well as the human subject. Figure 8C shows the mean value for the index of synergy for all experiments. A one-sample t-test between simulation and measurement results indicated that our proposed method could generate a similar mean value for the index of synergy as the measurement experiments [the far position task: t₍₇₎ = −0.10, P = 0.92; the near position task: t₍₇₎ = 1.03, P = 0.34]. A paired t-test demonstrated that the measured indices of synergy were not significantly different between the far and near position tasks [t₍₇₎ = −0.94, P = 0.38]. However, for the measurement experiments the mean peak of the tangential hand velocity for the far position task (2.93) tended to be larger than that for the near position task (2.73). Our proposed method also showed the same tendency (the far position task was 2.96, and the near position task was 2.59).

FIGURE 8

Figure 8. Index of synergy σ^* for the simulation and measurement experiments. The upper left (A) and right (B) graphs correspond to the far position and near position tasks. The red and blue lines indicate the mean value profiles of the index of synergy across all trials for the simulation experiments and across all subjects for the measurement experiments. The blue area denotes the standard deviation across all subjects. (B) The mean value for the index of synergy across the movement duration for all experiments. The n. s. indicates not significant (one-sample t-test).

Discussion

The aim of this study was to propose a control method that directly used the UCM and generated smooth and effortless movements while quantitatively comparing against human arm movements. We proposed a UCM reference feedback control for coordinating a redundant joint arm. Our proposed algorithm generated the target UCM from the given target hand trajectory step by step. The target joint combination was then selected in that space so as to minimize the cost function with respect to the input joint torque and torque change. Results of the simulation and measurement experiments for the one-dimensional target-tracking task were quantitatively compared.

According to the statistical results, our proposed method quantitatively reproduced kinematic and dynamic properties such as end postures, the tangential velocity of the hand and the joint torque. In addition, good agreement of the joint coordination was confirmed by the UCM analysis. For the one-dimensional target-tracking task, subjects showed a larger UCM component than an ORT component (Figure 7) indicating that multi-joints were coordinated to control the hand. This result suggests that there is a coordinated control system in the human visuomotor system allowing for task-irrelevant variability. Our proposed method could quantitatively reproduce the mean value for the UCM and ORT components, the index of synergy and the tendency depending on the posture difference. These results indicate that the UCM reference feedback control can generate human-like joint coordination.

In addition, these results could not be reproduced by a control framework in which a unique pattern of joint trajectories was determined. In robotics, a pseudo-inverse matrix is often used to control a redundant joint to generate unique target joint trajectories throughout the movement duration beforehand, and each joint tracks the target (Whitney, 1969). If a human used such a control strategy in the one-dimensional target-tracking task, the person could repetitively generate unique target joint trajectories since the same target hand trajectory and the same initial position were provided. In such a case, if the joints were disturbed by some noise or perturbation, they would converge to unique target joint trajectories. Using the above strategy, however, the UCM component does not increase with a small ORT component, as shown in Figures 7A,B. Therefore, at least for the one-dimensional target-tracking task, it is suggested that the human uses a strategy such as the UCM reference feedback control rather than a strategy that preliminarily generates unique joint angle trajectories for the whole movement duration.

Our UCM reference feedback control references online whether the variance of joint angles affects task achievement. In human behavior studies, it was reported that the response to the task-relevant perturbation was enhanced and the task-irrelevant perturbation was inhibited (Diedrichsen, 2007; Franklin and Wolpert, 2008). Moreover, such a modulation of the response can be achieved only using somatosensory information without visual information (Dimitriou et al., 2012). Thus, the human visuomotor system may have a neural representation of the UCM, especially in the somatosensory feedback loop (lower loop in Figure 1C). For the target-tracking task, the target performance variable representing the one-dimensional hand position was given by visual information. Thus, the UCM was generated by using visual information and could be utilized in the somatosensory feedback loop. Examining whether the UCM is represented in the brain and how to represent the UCM are considered important future works.

As shown in Figure 4, the error in the workspace (X–Y) coordinates indicates that the UCM reference feedback control could more accurately track the target than human subjects. In contrast, the ORT components, which were properties in the joint space and would be affected by the variability of the initial postures, showed quantitatively similar results between the simulation and measurement experiments (Figure 7F). Since the delay of human visual feedback system is relatively long, the tracking error of the human subjects was large.

Model Assumptions

In this paper, our UCM reference feedback control employed two assumptions. First, the three variables in Figure 1C (i.e., x^d, ${\theta }_{*}^d$, and τ) were exposed to noise. The noise of the input torque τrepresents the signal-dependent noise of a motor command (Harris and Wolpert, 1998). Additionally, the noise of the target hand position x^d represents the variability caused by target perception and localization with the human visual feedback system. Tseng et al. (2003) reported the relationship between the target size and the variability of joint angles in point-to-point reaching movements. They found that the variability of the joint angles increased with the size of the target. This suggests that the uncertainty of the perceived target affects the variability of joint angles. In addition, we assumed that the perception of the moving target would be more uncertain than the static target. Therefore, in our study such a relationship was represented by additive noise (Equation 16). The noise of the target joint combination ${\theta }_{*}^d$ represents the acquisition error of the UCM reference controller. As discussed above, we believe that the human CNS has an internal representation of information corresponding to the UCM, and we assume the acquisition uncertainty to be additive noise (Equation 17). The amplitudes of the noise e_τ, e_x, and e_θ affect the variability of the arm posture shown in Figure 3. The total variance of the joint is shown in Figure S1, and the profile of the input torque is shown in Figure 6.

Second, we used a cost function consisting of the input torque and torque change to select the unique joint combination for the target UCM. Alternatively, if we had used a cost function consisting of the minimum distance from the current joint combination to the target UCM in joint space it would be equivalent to using the pseudo-inverse matrix. In such a case, the required input torque would be much larger than the input torque as shown in Figure 6. Additionally, it would not be possible to reproduce the results of the measurement experiments. Therefore, a cost function based on the arm dynamics (e.g., torque and torque change) rather than a cost function related to the kinematics such as the pseudo-inverse matrix is required. Our cost function includes the torque change as a criterion for the smoothness of the human movement (Uno et al., 1989). Only the torque change term does not remain stationary since τ = 0 at the movement end; therefore, our cost function also includes the torque term (Equation 3) to generate effortless movements.

Comparison with the Other Control Method

When we used a cost function consisting of the sum of the squared norm of the joint angles instead of Equation (3), our proposed method was equivalent to a method in which the target joint angles were calculated for each time step using a pseudo-inverse matrix. Since the pseudo-inverse matrix does not consider the arm dynamics i.e., torque and torque change, the controller cannot smoothly, and effortlessly respond to noise and perturbation. Moreover, the pseudo-inverse matrix cannot reproduce a smooth bell-shape torque profile with motor noise as shown in Figure 6.

Figures 7, 8 show that our proposed method generates joint coordination while allowing variance that does not affect the task achievement. Such a “minimal intervention principle,” mentioned in the introduction, can also be realized within a framework of optimal feedback control (Todorov and Jordan, 2002b). The remarkable difference between a traditional optimal feedback control and our UCM reference feedback control is that the UCM is explicitly used in our method. In other words, our control scheme determines the optimal joint angle in UCM. This leads to the following advantages of the UCM reference feedback control. First, our proposed method can modularize an optimization process. In the traditional optimal feedback control framework, an input and output are simultaneously optimized using one cost function; e.g., E = ∫(x^TWx+τ^TQτ) dt. Therefore, the generated input and output depend on the weights used in the cost function.

For the one-dimensional target-tracking task, the accuracy of target tracking and the amplitude of the input torque are a weight-dependent trade-off. It would be difficult for the traditional framework to optimally solve this weight-dependent trade-off, especially for a nonlinear system. Moreover, even if the trade-off problem is solved, new weights are required for the situation where the task conditions are changed (e.g., a large input torque may or may not be required). On the other hand, our UCM reference feedback control can divide the optimization process into a target joint combination decision (Figure 1C: UCM reference controller) and an input torque decision (Figure 1C: inverse dynamics). By this modularization, the weight-dependent trade-off can be avoided and new weights are not required even if the task condition is changed.

Second, our method can be applied to a nonlinear system such as the human arm while the tradition optimal feedback control cannot be directly applied. Todorov et al. applied the framework of optimal feedback control to a nonlinear system using a hierarchical control framework in which the controlled object was linearized around the virtual target trajectory resulting in the convergence of the cost function to an optimal solution through iterative computation (Todorov et al., 2005; Li and Todorov, 2007; Liu and Todorov, 2009). However, we demonstrated that UCM reference feedback control can be directly applied a nonlinear system without any iteration. Finally, our method can be applied to many types of redundant control reported in the previous studies of UCM analysis.

For example, previous studies examined the relationships between the joint angles and the center of mass (Freitas et al., 2010), individual finger forces and the total force (Kang et al., 2004), electromyograms and the center of pressure (Wang et al., 2006), electromyograms and the center of mass (Robert et al., 2008), and angular momenta of different components and whole-body angular momentum (Robert et al., 2009). In addition, tasks that were examined include postural tasks (Krishnamoorthy et al., 2003; dos Santos et al., 2007; Scholz et al., 2007; Robert et al., 2008; Wu et al., 2009), hopping (Auyang et al., 2009; Yen et al., 2009), pistol shooting (Scholz et al., 2000), walking (Black et al., 2007; Robert et al., 2009), writing by hand (Shim et al., 2010), and carrying a cup of water (Togo et al., 2012). Our UCM reference feedback control can be applied to all the above tasks by only providing the appropriate target of the performance variable and the cost function. We believe that UCM reference feedback control can be developed from joint coordination to comprehensive “synergy.”

The most remarkable difference between UCM reference feedback control and optimal feedback control is the process of optimization. In optimal feedback control, the optimal gain is calculated so as to minimize some cost function over the whole duration of the task. In contrast, our UCM reference feedback control model selects an optimal point in the target UCM at each time step. However, at the present stage of research it is difficult to conclude which model is more appropriate for human control model. UCM reference feedback control could generate human-like coordinated arm movements, but our results do not establish that the human CNS solves an optimization problem step by step. Instead, we conjecture that the human CNS generates an output close to the solution of an optimization problem based on the properties of the musculoskeletal system and a neural structure.

Extending the framework of UCM reference feedback control and improving physiological appropriateness, we would like to discuss the validity of the control model for humans in future work. As an example of an earlier study of human control model with redundancy, Martin et al. (2009) proposed a redundant control model with the human muscle model. They solved the redundancy problem by utilizing a pseudo-inverse matrix with a kinematic cost function, while our control model can deal with a cost function related to the joint torque for the redundancy. We would like to extend the UCM reference feedback control as a control model for humans by considering the muscle model, adding biological delay and including neural representation from a neural network. In that case, the method for controlling redundant muscles is also an important and interesting problem. Recently, it was reported that human CNS controls redundant muscle through some hierarchical motor primitives (d'Avella et al., 2003; Bizzi and Cheung, 2013; d'Avella and Lacquaniti, 2013) rather than optimal control (de Rugy et al., 2012; Berger et al., 2013). In future work, we need to combine the “muscle synergy,” which is the motor primitive for redundant muscle; and the “motor synergy,” which is a neural mechanism for generating joint coordination.

Future Works

In this study, we proposed the UCM reference feedback control to reproduce human multi-joint coordination. There are some other control techniques for optimally determining unique joint combination at each time step, e.g., inverse kinematics (Buss, 2004; Peters et al., 2007; Liu et al., 2015) and optimal feedback control (Toussaint, 2009; Rückert et al., 2013). We consider statistical comparisons of the performance of the related control methods as important future work. Moreover, our proposed method requires an internal model, specifically an inverse dynamics model. In future work, we would embed the learning schema of an internal model, e.g., the feedback error learning (Kawato et al., 1987).

Author Contributions

Conceived and designed the model and experiments: ST. Performed the experiments: ST. Analyzed the data: ST, TK, and YU. Wrote the paper: ST, TK, and YU. All authors gave final approval for publication.

Funding

This research was supported by a Japan Society for the Promotion of Science (JSPS) Grant-in-Aid for JSPS Fellows 267576, and Grants-in-Aid for Scientific Research (B) No. 26289129 and 26280101.

Conflict of Interest Statement

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

Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fncom.2016.00069

Figure S1. Total variance for the simulation and measurement experiments. The upper left (A) and right (B) graphs correspond to the far position and near position tasks. The red and blue lines indicate the mean values profiles for the total variance across all trials for the simulation experiments and across all subjects for the measurement experiments. The blue area denotes the standard deviation across all subjects. (C) The mean value for the total variance across the movement duration for all experiments. The n. s. indicates not significant (one-sample t-test).

References

Auyang, A. G., Yen, J. T., and Chang, Y. H. (2009). Neuromechanical stabilization of leg length and orientation through interjoint compensation during human hopping. Exp. Brain Res. 192, 253–264. doi: 10.1007/s00221-008-1582-7

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Bernstein, N. A. (1967). The Coordination and Regulation of Movements. Oxford: Pergamon.

Bizzi, E., and Cheung, V. C. K. (2013). The neural origin of muscle synergies. Front. Comput. Neurosci. 7:51. doi: 10.3389/fncom.2013.00051

PubMed Abstract | CrossRef Full Text | Google Scholar

Black, D. P., Smith, B. A., Wu, J., and Ulrich, B. D. (2007). Uncontrolled manifold analysis of segmental angle variability during walking: preadolescents with and without down syndrome. Exp. Brain Res. 183, 511–521. doi: 10.1007/s00221-007-1066-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Buss, S. R. (2004). Introduction to inverse kinematics with Jacobian transpose, pseudoinverse and damped least squares methods. IEEE J. Robotics Automation 17, 1–19.

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

de Rugy, A., Loeb, G. E., and Carroll, T. J. (2012). Muscle coordination is habitual rather than optimal. J. Neurosci. 32, 7384–7391. doi: 10.1523/JNEUROSCI.5792-11.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Diedrichsen, J. (2007). Optimal task-dependent changes of bimanual feedback control and adaptation. Curr. Biol. 17, 1675–1679. doi: 10.1016/j.cub.2007.08.051

PubMed Abstract | CrossRef Full Text | Google Scholar

Dimitriou, M., Franklin, D. W., and Wolpert, D. M. (2012). Task-dependent coordination of rapid bimanual motor responses. J. Neurophysiol. 107, 890–901. doi: 10.1152/jn.00787.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Domkin, D., Laczko, J., Djupsjöbacka, M., Jaric, S., and Latash, M. L. (2005). Joint angle variability in 3D bimanual pointing: uncontrolled manifold analysis. Exp. Brain Res. 163, 44–57. doi: 10.1007/s00221-004-2137-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Domkin, D., Laczko, J., Jaric, S., Johanson, H., and Latash, M. L. (2002). Structure of joint variability in bimanual pointing tasks. Exp. Brain Res. 143, 11–23. doi: 10.1007/s00221-001-0944-1

PubMed Abstract | CrossRef Full Text | Google Scholar

dos Santos, A. D., Slomka, K., Zatsiorsky, V. M., and Latash, M. L. (2007). Muscle modes and synergies during voluntary body sway. Exp. Brain Res. 179, 533–550. doi: 10.1007/s00221-006-0812-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Flash, T., and Hogan, N. (1985). The coordination of arm movements: an experimentally confirmed mathematical model. J. Neurosci. 5, 1688–1703.

PubMed Abstract | Google Scholar

Franklin, D. W., and Wolpert, D. M. (2008). Specificity of reflex adaptation for task-relevant variability. J. Neurosci. 28, 14165–14175. doi: 10.1523/JNEUROSCI.4406-08.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Freitas, S. M. S. F., Scholz, J. P., and Latash, M. L. (2010). Analysis of joint variance related to voluntary whole-body movements performed in standing. J. Neurosci. Methods 188, 89–96. doi: 10.1016/j.jneumeth.2010.01.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Harris, C. M., and Wolpert, D. M. (1998). Signal-dependent noise determines motor planning. Nature 394, 780–784. doi: 10.1038/29528

PubMed Abstract | CrossRef Full Text | Google Scholar

Kang, N., Shinohara, M., Zatsiorsky, V. M., and Latash, M. L. (2004). Learning multi-finger synergies: an uncontrolled manifold analysis. Exp. Brain Res. 157, 336–350. doi: 10.1007/s00221-004-1850-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Kawato, M., Furukawa, K., and Suzuki, R. (1987). A hierarchical neural-network model for control and learning of voluntary movement. Biol. Cybern. 57, 169–185. doi: 10.1007/BF00364149

PubMed Abstract | CrossRef Full Text | Google Scholar

Krishnamoorthy, V., Latash, M. L., Scholz, J. P., and Zatsiorsky, V. M. (2003). Muscle synergies during shifts of the center of pressure by standing persons. Exp. Brain Res. 152, 281–292. doi: 10.1007/s00221-003-1574-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Latash, M. L. (2010). Stages in learning motor synergies: a view based on the equilibrium-point hypothesis. Hum. Mov. Sci. 29, 642–654. doi: 10.1016/j.humov.2009.11.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Latash, M. L. (2012). The bliss (not the problem) of motor abundance (not redundancy). Exp. Brain Res. 217, 1–5. doi: 10.1007/s00221-012-3000-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Latash, M. L., Scholz, J. P., and Schöner, G. (2002). Motor control strategies revealed in the structure of motor variability. Exerc. Sport Sci. Rev. 30, 26–31. doi: 10.1097/00003677-200201000-00006

PubMed Abstract | CrossRef Full Text | Google Scholar

Latash, M. L., Scholz, J. P., and Schöner, G. (2007). Toward a new theory of motor synergies. Motor Control 11, 276–308.

PubMed Abstract | Google Scholar

Li, W., and Todorov, E. (2007). Iterative linearization methods for approximately optimal control and estimation of non-linear stochastic system. Int. J. Control 80, 1439–1453. doi: 10.1080/00207170701364913

CrossRef Full Text | Google Scholar

Liu, D., and Todorov, E. (2009). “Hierarchical optimal control of a 7-dof arm model,” in IEEE Symposium on Adaptive Dynamic Programming and Reinforcement Learning (Nashville, TN), 50–57.

PubMed Abstract | Google Scholar

Liu, M., Hak, S., and Padois, V. (2015). “Generalized projector for task priority transitions during hierarchical control,” in 2015 IEEE International Conference on Robotics and Automation (ICRA) (Seattle, WA: IEEE), 768–773.

Martin, V., Scholz, J. P., and Schöner, G. (2009). Redundancy, self-motion, and motor control. Neural Comput. 21, 1371–1414. doi: 10.1162/neco.2008.01-08-698

PubMed Abstract | CrossRef Full Text | Google Scholar

Peters, J., Mistry, M., Udwadia, F., Nakanishi, J., and Schaal, S. (2007). A unifying framework for robot control with redundant DOFs. Auton. Robots 24, 1–12. doi: 10.1007/s10514-007-9051-x

CrossRef Full Text | Google Scholar

Robert, T., Bennett, B. C., Russell, S. D., Zirker, C. A., and Abel, M. F. (2009). Angular momentum synergies during walking. Exp. Brain Res. 197, 185–197. doi: 10.1007/s00221-009-1904-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Robert, T., Zatsiorsky, V. M., and Latash, M. L. (2008). Multimuscle synergies in an unusual postural task: quick shear force production. Exp. Brain Res. 187, 237–253. doi: 10.1007/s00221-008-1299-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Rückert, E. A., Neumann, G., Toussaint, M., and Maass, W. (2013). Learned graphical models for probabilistic planning provide a new class of movement primitives. Front. Comput. Neurosci. 6:97. doi: 10.3389/fncom.2012.00097

PubMed Abstract | CrossRef Full Text

Scholz, J. P., and Schöner, G. (1999). The uncontrolled manifold concept: identifying control variables for a functional task. Exp. Brain Res. 126, 289–306. doi: 10.1007/s002210050738

PubMed Abstract | CrossRef Full Text | Google Scholar

Scholz, J. P., Schöner, G., Hsu, W. L., Jeka, J. J., Horak, F., and Martin, V. (2007). Motor equivalent control of the center of mass in response to support surface perturbations. Exp. Brain Res. 180, 163–179. doi: 10.1007/s00221-006-0848-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Scholz, J. P., Schöner, G., and Latash, M. L. (2000). Identifying the control structure of multijoint coordination during pistol shooting. Exp. Brain Res. 135, 382–404. doi: 10.1007/s002210000540

PubMed Abstract | CrossRef Full Text | Google Scholar

Shim, J. K., Hooke, A. W., Kim, Y. S., Park, J., Karol, S., and Kim, Y. H. (2010). Handwriting: hand-pen contact force synergies in circle drawing tasks. J. Biomech. 43, 2249–2253. doi: 10.1016/j.jbiomech.2010.04.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Todorov, E., and Jordan, M. I. (2002a). “A minimal intervention principle for coordinated movement,” in Advances in Neural Information Processing Systems 15, eds S. Becker, S. Thrun, S. and K. Obermayer (Cambridge, MA: MIT), 27–34.

Todorov, E., and Jordan, M. I. (2002b). Optimal feedback control as a theory of motor coordination. Nat. Neurosci. 5, 1226–1235. doi: 10.1038/nn963

PubMed Abstract | CrossRef Full Text | Google Scholar

Todorov, E., Li, W., and Pan, X. (2005). From task parameters to motor synergies: a hierarchical framework for approximately optimal control of redundant manipulators. J. Robot. Syst. 22, 691–710. doi: 10.1002/rob.20093

PubMed Abstract | CrossRef Full Text | Google Scholar

Togo, S., Kagawa, T., and Uno, Y. (2012). Motor synergies for dampening hand vibration during human walking. Exp. Brain Res. 216, 81–90. doi: 10.1007/s00221-011-2909-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Toussaint, M. (2009). “Robot trajectory optimization using approximate inference,” in Proceedings of the 26th Annual International Conference on Machine Learning (New York, NY), 1049–1056.

Tseng, Y. W., Scholz, J. P., and Schöner, G. (2002). Goal equivalent joint coordination in pointing: affect of vision and arm dominance. Motor Control 6, 183–204.

PubMed Abstract | Google Scholar

Tseng, Y. W., Scholz, J. P., Schöner, G., and Hotchkiss, L. (2003). Effect of accuracy constraint on joint coordination during pointing movements. Exp. Brain Res. 149, 276–288. doi: 10.1007/s00221-002-1357-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Uno, Y., Kawato, M., and Suzuki, R. (1989). Formation and control of optimal trajectory in human multi-joint arm movement - minimum torque-change model. Biol. Cybern. 61, 89–101. doi: 10.1007/BF00204593

PubMed Abstract | CrossRef Full Text | Google Scholar

Verrel, J. (2010). Distributional properties and variance-stabilizing transformations for measures of uncontrolled manifold effects. J. Neurosci. Methods 191, 166–170. doi: 10.1016/j.jneumeth.2010.06.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Asaka, T., Zatsiorsky, V. M., and Latash, M. L. (2006). Muscle synergies during voluntary body sway: combining across-trials and within-a-trial analysis. Exp. Brain Res. 174, 679–693. doi: 10.1007/s00221-006-0513-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Whitney, D. (1969). Resolved motion rate control of manipulators and human prostheses. IEEE Trans. Man-Machine Syst. 10, 47–53. doi: 10.1109/TMMS.1969.299896

CrossRef Full Text | Google Scholar

Winter, D. A. (2004). Biomechanics and Motor Control of Human Movement, 3rd Edn. New York, NY: John Wiley & Sons, Inc.

Google Scholar

Wu, J., McKay, S., and Angulo-Barroso, R. (2009). Center of mass control and multi-segment coordination in children during quiet stance. Exp. Brain Res. 196, 329–339. doi: 10.1007/s00221-009-1852-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, J. F., Scholz, J. P., and Latash, M. L. (2007). The role of kinematic redundancy in adaptation of reaching. Exp. Brain Res. 176, 54–69. doi: 10.1007/s00221-006-0602-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Yen, J. T., Auyang, A. G., and Chang, Y. H. (2009). Joint-level kinetic redundancy is exploited to control limb level forces during human hopping. Exp. Brain Res. 196, 439–451. doi: 10.1007/s00221-009-1868-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Yen, J. T., and Chang, Y. H. (2009). Rate-dependent control strategies stabilize limb forces during human locomotion. J. R. Soc. Interface 7, 801–810. doi: 10.1098/rsif.2009.0296

PubMed Abstract | CrossRef Full Text | Google Scholar

Appendix

A. Calculation of Cost Function

Using constant terms B and F which are associated with the torque and torque change terms, the specific constant term A in Eq. (14) is as follows:

$\begin{array}{l}{A}_1=B_1^2+B_4^2+B_7^2+F_1^2+F_4^2+F_7^2\\ A_2=2(B_1{B}_3+B_4{B}_6+B_7{B}_9+F_1{F}_3+F_4{F}_6+F_7{F}_9)\\ A_3=2(B_1{B}_2+B_4{B}_5+B_7{B}_8+F_1{F}_2+F_4{F}_5+F_7{F}_8)\\ A_4=2(B_2{B}_3+B_5{B}_6+B_8{B}_9+F_2{F}_3+F_5{F}_6+F_8{F}_9)\\ A_5=B_2^2+B_5^2+B_8^2+F_2^2+F_5^2+F_8^2\\ A_6=B_3^2+B_6^2+B_9^2+F_3^2+F_6^2+F_9^2.\left(\text{A}1\right)\end{array}$

The specific constant terms B and C are as follows:

$\begin{array}{l}{B}_1=\frac{W_s}{{(\Delta t)}^2}\left(-\frac{J_2^{†}}{J_1^{†}}{M}_{11}+M_{12}\right)\\ B_2=\frac{W_s}{{(\Delta t)}^2}\left(-\frac{J_3^{†}}{J_1^{†}}{M}_{11}+M_{13}\right)\\ B_3=\frac{W_s}{{(\Delta t)}^2}\left(-M_{12}{\theta }_e^{†}-M_{13}{\theta }_w^{†}-\frac{x_d-x^{†}+J_2^{†}{\theta }_e^{†}+J_3^{†}{\theta }_w^{†}}{J_1^{†}}{M}_{11}\right)\\ B_4=\frac{W_e}{{(\Delta t)}^2}\left(-\frac{J_2^{†}}{J_1^{†}}{M}_{21}+M_{22}\right)\\ B_5=\frac{W_e}{{(\Delta t)}^2}\left(-\frac{J_3^{†}}{J_1^{†}}{M}_{21}+M_{23}\right)\\ B_6=\frac{W_e}{{(\Delta t)}^2}\left(-M_{22}{\theta }_e^{†}-M_{23}{\theta }_w^{†}-\frac{x_d-x^{†}+J_2^{†}{\theta }_e^{†}+J_3^{†}{\theta }_w^{†}}{J_1^{†}}{M}_{21}\right)\\ B_7=\frac{W_w}{{(\Delta t)}^2}\left(-\frac{J_2^{†}}{J_1^{†}}{M}_{31}+M_{32}\right)\\ B_8=\frac{W_w}{{(\Delta t)}^2}\left(-\frac{J_3^{†}}{J_1^{†}}{M}_{31}+M_{33}\right)\\ B_9=\frac{W_w}{{(\Delta t)}^2}\left(-M_{32}{\theta }_e^{†}-M_{33}{\theta }_w^{†}-\frac{x_d-x^{†}+J_2^{†}{\theta }_e^{†}+J_3^{†}{\theta }_w^{†}}{J_1^{†}}{M}_{31}\right)\left(\text{A}2\right)\end{array}$

$\begin{array}{l}{F}_1=\frac{Q_s}{{(\Delta t)}^3}\left(-\frac{J_2^{‡}}{J_1^{‡}}{M}_{11}+M_{12}\right)\\ F_2=\frac{Q_s}{{(\Delta t)}^3}\left(-\frac{J_3^{‡}}{J_1^{‡}}{M}_{11}+M_{13}\right)\\ F_3=\frac{Q_s}{{(\Delta t)}^3}\left(-M_{12}{\theta }_e^{‡}-M_{13}{\theta }_w^{‡}-\frac{x_d-x^{‡}+J_2^{‡}{\theta }_e^{‡}+J_3^{‡}{\theta }_w^{‡}}{J_1^{‡}}{M}_{11}\right)\\ F_4=\frac{Q_e}{{(\Delta t)}^3}\left(-\frac{J_2^{‡}}{J_1^{‡}}{M}_{21}+M_{22}\right)\\ F_5=\frac{Q_e}{{(\Delta t)}^3}\left(-\frac{J_3^{‡}}{J_1^{‡}}{M}_{21}+M_{23}\right)\\ F_6=\frac{Q_e}{{(\Delta t)}^3}\left(-M_{22}{\theta }_e^{‡}-M_{23}{\theta }_w^{‡}-\frac{x_d-x^{‡}+J_2^{‡}{\theta }_e^{‡}+J_3^{‡}{\theta }_w^{‡}}{J_1^{‡}}{M}_{21}\right)\\ F_7=\frac{Q_w}{{(\Delta t)}^3}\left(-\frac{J_2^{‡}}{J_1^{‡}}{M}_{31}+M_{32}\right)\\ F_8=\frac{Q_w}{{(\Delta t)}^3}\left(-\frac{J_3^{‡}}{J_1^{‡}}{M}_{31}+M_{33}\right)\\ F_9=\frac{Q_w}{{(\Delta t)}^3}\left(-M_{32}{\theta }_e^{‡}-M_{33}{\theta }_w^{‡}-\frac{x_d-x^{‡}+J_2^{‡}{\theta }_e^{‡}+J_3^{‡}{\theta }_w^{‡}}{J_1^{‡}}{M}_{31}\right).\left(\text{A}3\right)\end{array}$

The specific components of the inertia matrix M and Jacobian J are as follows:

$\begin{array}{l}{M}_{11}=m_1{l}_{g1}^2+I_1+m_2\left(l_1^2+l_{g2}^2+2l_1{l}_{g2}\cos {\theta }_e\right)+I_2\\ +m_3(l_1^2+l_2^2+l_{g3}^2+2l_1{l}_2\cos {\theta }_e+2l_1{l}_{g3}\cos ({\theta }_e+{\theta }_w)\\ +2l_2{l}_{g3}\cos {\theta }_w)+I_3\\ M_{12}=M_{21}=m_2\left(l_{g2}^2+l_1{l}_{g2}\cos {\theta }_e\right)+I_2\\ +m_3(l_2^2+l_{g3}^2+l_1{l}_2\cos {\theta }_e+l_1{l}_{g3}\cos ({\theta }_e+{\theta }_w)\\ +2l_2{l}_{g3}\cos {\theta }_w)+I_3\\ M_{13}=M_{31}=m_3\left(l_{g3}^2+l_1{l}_{g3}\text{cos}({\theta }_e+{\theta }_w)+l_2{l}_{g3}\cos {\theta }_w\right)+I_3\\ M_{22}=m_2{l}_{g2}^2+I_2+m_3\left(l_2^2+l_{g3}^2+2l_2{l}_{g3}\cos {\theta }_w\right)+I_3\\ M_{23}=M_{32}=m_3\left(l_{g3}^2+l_2{l}_{g3}\cos {\theta }_w\right)+I_3\left(\text{A}4\right)\end{array}$

$\begin{array}{l}{M}_{33}=m_3{l}_{g3}^2+I_3,\\ J_1^{†}=J_2^{†}-l_{1}sin{\theta }_s^{†}\left(\text{A}5\right)\end{array}$

$\begin{array}{l}{J}_2^{†}=J_3^{†}-l_{2}sin({\theta }_s^{†}+{\theta }_e^{†})\\ J_3^{†}=-l_{3}sin({\theta }_s^{†}+{\theta }_e^{†}+{\theta }_w^{†}),\\ J_1^{‡}=J_2^{‡}-l_{1}sin{\theta }_s^{‡}\\ J_2^{‡}=J_3^{‡}-l_{2}sin({\theta }_s^{‡}+{\theta }_e^{‡})\\ J_3^{‡}=-l_{3}sin({\theta }_s^{‡}+{\theta }_e^{‡}+{\theta }_w^{‡}).\left(\text{A}6\right)\end{array}$

Using the constant term A, the extreme values ${\theta }_s^{d*}$ and ${\theta }_s^{d*}$ are represented as follows:

$\begin{array}{l}{\theta }_e^{d*}=\frac{A_3{A}_4-{2A}_2{A}_5}{4A_1{A}_5-A_3^2}{\theta }_e^{d*}=\frac{A_2{A}_3-{2A}_1{A}_4}{4A_1{A}_5-A_3^2}.\left(\text{A}7\right)\end{array}$

B. Supplementary Figures

Figures S1A and S1B show profiles for the total variance. The red and blue solid lines indicate the results of the simulation and measurement experiments. The blue area denotes the standard deviation across all subjects. In both the simulation and measurement experiments, the total variance gradually increased from the movement start to end. Our proposed method could generate same tendency. Figure S1C shows the mean total variance for all experiments. A one-sample t-test indicated that our proposed method could generate a similar mean value for the total variance to the measurement experiments [the far position task: t₍₇₎ = −0.32, P = 0.76; the near position task: t₍₇₎ = 0.59, P = 0.57]. A paired t-test demonstrated that the mean total variances were not significantly different between the far and near position tasks [t₍₇₎ = −0.56, P = 0.60]. However, for the measurement experiments the mean total variance for the far position task [0.0012 (rad)²] tended to be larger than that in the near position task [9.95 × 10⁻⁴ (rad)²]. Our proposed method also showed the same tendency [the far position task: 0.0012 (rad)²; the near position task: 9.32 × 10⁻⁴ (rad)²].