## Abstract

Here we computationally investigate how encumbering the hand could alter predictions made by the minimum torque change (MTC) and minimum endpoint variance hypotheses (MEPV) of movement planning. After minutes of training, people have made arm trajectories in a robot-generated viscous force field that were similar to previous baseline trajectories without the force field. We simulate the human arm interacting with this viscous load. We found that the viscous forces clearly differentiated MTC and MEPV predictions from both minimum-jerk predictions and from human behavior. We conclude that learned behavior in the viscous environment could arise from minimizing kinematic costs but could not arise from a minimization of either torque change or endpoint variance.

## INTRODUCTION

Despite decades of basic scientific research attempting to explain normal human motor control, the fundamental neural and computational bases of even the simplest motor behaviors remain poorly understood (Flash and Sejnowski 2001; Wolpert and Ghahramani 2000). The experimental and computational investigations of motor control have primarily concerned unencumbered horizontal movements of the supported arm (Desmurget et al. 1998). Even in these limited explorations, several theories of motor planning currently compete. Different theories posit that humans aim to optimize kinematic, dynamic or variance costs in the planning of movement. Although each theory explains some facets of horizontal movement trajectories, none can fully explain real human behavior. A more complete exploration of human motor behaviors across a variety of tasks, and simulations that consider a broader set of possible factors, are needed to understand to what degree each of these theories explain normal human behavior (Englebrecht 2001).

Optimal control strategies have been implemented to overcome the ill-posed nature of motor planning by minimizing cost functions to identify potential motor planning criteria. Optimal controllers minimizing jerk (Flash and Hogan 1985) and torque change (Uno et al. 1989) consider desired behavioral performance, which can be experimentally ascertained by measuring mean behavior. A more recent theory (Todorov and Jordan 2002) posits that people do not optimize their mean behaviors but rather aim to optimize their variance around the metrics of behavior most relevant to the task. In the case of point-to-point movement, this suggests that people focus on minimizing the error of their hands at the endpoint of movement (Harris and Wolpert 1998). These movement strategies interrogate the influence of mean extrinsic kinematics (jerk), mean intrinsic dynamics (torque change), and variance of extrinsic kinematics (endpoint variance) on motor planning and replicate the straight trajectories and Gaussian-shaped speed profiles of horizontal plane, point-to-point reaching movements quite well. Exploration of other motor behaviors reveal that the minimization of torque change fits longer movements that feature joint reversals (Uno et al. 1989); the minimum endpoint error hypothesis fits well the avoidance of obstacles (Hamilton and Wolpert 2002) and the differing variance observed in different horizontal movement directions (van Beers et al. 2004). Altering the positional characteristics (length, path, or direction) of movement have therefore revealed subtleties of motor behavior that are better mimicked by a particular motor planning hypothesis.

Here we examine whether altering the dynamic demands of movement differentiate between the predictions of minimum jerk (MJ), minimum torque change (MTC), or minimum endpoint variance (MEPV) hypotheses. People have, after hundreds of movements of training, made movements in viscous forces that closely resemble their movements before the application of forces (Shadmehr and Brashers-Krug 1997; Thoroughman and Shadmehr 2000). We therefore simulated these additional dynamics and asked whether viscous forces clearly differentiated predictions of the MTC or MEPV hypotheses from the straight-line trajectories predicted by the MJ hypothesis or from previously recorded human behavior.

## METHODS

### Canonical simulation

We simulated the human arm as a two-link system moving in the horizontal plane, supported against gravity. The arm moved either unencumbered or while holding a manipulandum (Fig. 1) that generated a viscous force field (as modeled in Bhushan and Shadmehr 1999). All simulations (Fig. 2) and analyses were conducted in Matlab and Simulink (Mathworks, Natick, MA).

The simulations generated 10-cm movements toward eight targets away from a central starting location. Movements were 500 ms in duration.

The first part of the model (Fig. 2) generated a desired neural activation off-line. Desired trajectories are generated by the optimization of costs (see *Spline parameterization*). Given a desired trajectory, an inverse dynamic transformation computes desired joint torques. This transformation accounted for loads generated or by the manipulandum. The equation of motion used for inverse dynamics was (1) where τ was the 2 × 1 joint torque (shoulder and elbow), **A** was the inertia matrix, **B** was the Coriolis/centripetal matrix, **J**_{h} was the Jacobian of the arm, and **F** was the robot-generated force. We changed the **A** and **B** to add the mechanics of the passive dynamics of the robot (see appendix).

### Torque change

We calculated torque change, throughout the desired movement duration, using the deterministic desired torque from the inverse dynamic transformation (Uno et al. 1989): (2) where subscripts s and e correspond to shoulder and elbow joints.

### Simulating noisy movements

To consider the effect of signal-dependent noise on endpoint variance, we simulated the transformation of desired muscle torques into actual movement. The desired muscle torques were each passed through derivative filters that served to invert the muscle model (detailed in the following text). The inverse muscle transformation generated a desired neural activation, which was used to drive both the torque change and endpoint variance minimizations. The generation of the desired activation concluded the off-line processing of desired signals; all following transformations were trial-by-trial simulations of the actual neural control of movement under the complications of signal-dependent noise.

The resultant desired neural activation was subject to signal-dependent noise (as in Harris and Wolpert 1998), the SD of which scaled with the desired activation. The canonical model featured signal-dependent noise with 30% SD and the noise was updated at 1 kHz, the refresh rate of the simulation. The canonical model also featured no feedback; therefore the movement followed only the descending, noisy neural activations. The activations were first filtered through the model muscles, consisting of two first-order filters, in series, with time constants 30 and 40 ms (van der Helm and Rozendaal 2000). These filters generated the muscle torque to move the arm.

### Endpoint variance

We used the time interval from the desired movement end (*t*_{f} = 500 ms) through an additional time window (Δ*t* = 50 ms) thereafter to calculate the variance of the simulated movement endpoint (Harris and Wolpert 1998) (3) where *s* is the two-dimensional hand position [*x, y*], the subscript d denotes desired, and the outer angled brackets denote an average over 100 trials of the noisy movement generation. The desired endpoint at the end of the movement was at the center of the target, 10 cm from the start.

### Sensitivity analysis

We generated our optimizations using a canonical model of the central controller, which mimics the transformations calculated by the CNS. This model presumed solely predictive control, generated a high level of signal-dependent noise, and produced the noise at 1 kHz, the same frequency as the sampling rate of the simulation. All of these presumptions affected the transformation of noisy activation into trajectory generation and could have affected the movements that minimize endpoint variance. We altered the canonical model to determine the robustness of the curvature of minimum endpoint variance trajectories.

We altered the model in four ways. Individually we added instantaneous, noiseless stiff and/or viscous feedback (termed mechanical feedback); added delayed, noisy stiff and/or viscous feedback (termed neural feedback); reduced the variance of signal-dependent noise from 30 to 5%; and slowed the frequency of noise generation from 1 kHz to 6 Hz.

The mechanical feedback generated error-dependent torques τ_{fback} = −**K**(*q* − *q*_{d}) − **V**(q̇−q̇_{d}), where the stiffness **K** = {−15, −6; −6, −16} N.m/rad and the viscosity **V** = {−2.3, −0.9; −0.9, −2.4} N.m.s/rad (Shadmehr and Mussa-Ivaldi 1994). Neural feedback consisted of delaying the mechanical feedback by 35 ms (van der Helm and Rozendaal 2000), then adding it to the desired neural activation such that the feedback is subject to the signal-dependent noise. The mechanical and neural feedback, respectively, considered possible contributions from intrinsic muscle stiffness and from neural feedback loops.

### Spline parameterization of trajectory space

The simulated controller depended on desired trajectory to produce the desired neural signal and ultimately the movement. We utilized quintic splines with two fixed knots at movement start and end and three movable knots equally spaced throughout the desired movement time. These knots enabled a parameterization of trajectory space; we aimed to minimize torque change and endpoint variance in each model system by varying the knot positions. As previously observed (Hamilton and Wolpert 2002), standard quintic splines do not feature enough degrees of freedom to bound the solution to have zero velocity, acceleration, jerk, and snap at the endpoints. Such quiescence is necessary to simulate our fourth-order system, consisting of second-order muscles driving second-order dynamics. To alleviate this problem, we used the Spline Toolbox of Matlab to add floating knots between the five specified knots. Each new knot required an additional quintic spline, but we did not specify the location of the floating knots. These floating knots added necessary degrees of freedom to enable additional boundary conditions, which in turn made possible the actual simulation of the fourth-order system.

### Toy model

To explore the fundamental mathematics and physics underlying our canonical model, we also generated a toy model. In this model, linear force actuators controlled a point mass that moved in the viscous curl field. We found the trajectory that minimized the time derivative of the linear force generation. We then implemented the second-order linear muscle model and signal-dependent noise, as described in the preceding text, and found the trajectory that minimized endpoint variance.

### Human data

We used data from subjects previously analyzed in Thoroughman and Shadmehr (2000). Twenty-four subjects made 500-ms, 10-cm reaching movements while holding a manipulandum (see Thoroughman and Shadmehr 2000 for details). A velocity-dependent force field ({0, 13; −13, 0} N.s/m) was active for 576 movements. Within target direction, we averaged over the last 128 overall movements where the performance had reached a plateau (Thoroughman and Shadmehr 1999). As described in the preceding text, quintic splines were generated to represent the desired trajectory. The knots were adjusted to minimize the root-mean-square difference between a quintic spline and the human movement.

### Data reduction

An emergent feature of our simulations was a C-shaped curvature throughout the movement. We therefore chose a fractional curvature metric that reflected the maximal perpendicular displacement (the component of displacement perpendicular to the target direction) divided by the distance between movement start and target (which was always 10 cm).

## RESULTS

### Simulations

We began our simulations by confirming the minimization of either torque change or endpoint variance generated straight movements for the unencumbered arm. We observed trajectories for both minimizations that were largely straight; across all eight directions (Fig. 3*A*), the mean curvature (maximal perpendicular displacement divided by distance to target) was 0.00142 for minimum torque change (MTC, gray lines) and 0.0135 for minimum endpoint variance (MEPV, black lines), and the maximal curvature among all directions was 0.0136 in the 135° direction for MTC and 0.0253 in the 0° direction for MEPV. The trajectories produced by the two hypotheses mimicked the hallmark straight hand path generated by people in short, unencumbered, horizontal, point-to-point arm movements and predicted by the minimum-jerk hypothesis.

Next we investigated whether external dynamics interacting with the arm would differentiate predictions made by minimizations of jerk, torque change, and endpoint variance. We simulated the passive dynamics of a manipulandum and the perturbing forces of a clockwise viscous force field. Because jerk measures the kinematics of movement, external dynamics do not alter the perfectly straight trajectories and bell-shaped velocity profiles of jerk minimization. We therefore investigated how a perfectly trained controller, with full knowledge of the robot and force field dynamic demands, would minimize torque change or endpoint variance in the execution of point-to-point reaches.

We discovered that MTC trajectories curved markedly counterclockwise when the viscous forces pushed clockwise (Fig. 3*B*, gray). The mean curvature across directions was 0.201; the maximum curvature among all directions was 0.386 in direction 0°. The MEPV trajectories featured even more pronounced curvature (Fig. 3*B*, black): 0.241 on average, 0.474 at most in direction 0°. We discovered that the dynamics generated by the robotic manipulandum were sufficient to distinguish the predictions these two hypotheses from minimum jerk.

### Sensitivity

We simulated a canonical model of the central controller, which mimics the transformations calculated by the CNS. This model presumed solely predictive control, generated a high level of signal-dependent noise, and produced the noise at 1 kHz, the same frequency as the sampling rate of the simulation. We individually introduced either mechanical or neural feedback (as detailed in methods), reduced the SD of the noise generation, or slowed the noise production to 6 Hz. These changes affected only the MEPV simulations, as MTC depended only on the desired movement. The trajectories predicted using one of the controllers is shown in Fig. 3*C*. This controller featured stiffness and viscosity that were implemented through neural activation, such that the feedback control signals were delayed and subject to signal-dependent noise. Although these changes altered fine details of the MEPV trajectories, all optimized trajectories in the viscous environment featured average curvature around or >0.20 and maximal curvature >0.35 (Table 1).

### Fundamental considerations of optimal solutions using toy model

Each of these sets of trajectories featured certain directions with marked curvature (such as 0°) and others with subtler curvature (such as 90°). These differences were caused by anisotropic dynamic parameters and by the control of a linked arm that possesses inherent strong nonlinearities (see appendix). We generated a toy model with linear force actuators to dispose of these anisotropies and nonlinearities and to discover fundamental properties of the optimal solutions. Note that the linearity and isotropy of the toy model possessed no information on target direction; the results apply equally well to any movement direction.

The trajectory that generated minimum force change was markedly curved opposite the direction of force (Fig. 4, gray), with a curvature of 0.19 (indicating that the maximal perpendicular displacement was 1.9 cm and that the distance from start to target was 10 cm). The trajectory that minimized endpoint variance was even more curved in the same direction (Fig. 4, black), with a curvature of 0.33. A heuristic explanation of these curvatures is that early in movement, the forces push equally hard regardless of movement direction. Midway through movement completion, these curved movements take advantage of the perturbing viscous force to push the mass toward the target; a straight trajectory instead requires the controller to counter the viscous force in its entirety. This savings allows for lesser force change and lesser activation, hence lesser signal-dependent noise, hence smaller endpoint variance. The full canonical model deviates from these curvatures because of the nonlinearities and anisotropies of the human and robotic arms.

### Human behavior

We last compared these optimizations with human performance. Subjects (*n* = 24) learned on 1 day the passive dynamics of the robot, then on a second day reached asymptotic performance after 500 movements of training in the clockwise viscous force field. Some subjects curved in a clockwise direction, others in a counterclockwise direction, and some generated S-shaped trajectories featuring both clockwise and counterclockwise curvatures. The MTC and MEPV trajectories featured prominent counterclockwise curvature; we aimed to test the null hypothesis that subjects moved the same way as these curved trajectories. We therefore calculated the maximal curvature in the counterclockwise direction for mean trajectories across subjects and the 95% confidence interval of the mean (Fig. 3*D*). Although these movements did feature counterclockwise curvature (Table 1), the curvature generated by the torque change optimization was eight times greater and the curvature generated by the endpoint variance optimization was 10 times greater. The torque change optimization and all variants of the endpoint variance optimization were significantly different from human performance (*P* < 0.001).

## DISCUSSION

Three prominent theories (Harris and Wolpert 1998; Hogan 1984; Uno et al. 1989) have posited that the CNS generates motor behaviors, such as arm movements, by optimizing a cost incurred during movement. These three strategies share a motivating principle: most movement tasks specify only a goal of movement without specifying the full time series of body positions and joint torques necessary to perform the task.

If the CNS indeed minimizes a single hypothesized cost in its motor behavior, then computational simulations of controllers that minimize that cost should mimic human behavior. Our first simulations determined trajectories for the unencumbered arm, supported against gravity, reaching in the horizontal plane. We confirmed the conclusions of the original proponents of the theories: optimization of torque change or endpoint variance replicates the straight trajectories classically (Morasso 1981) generated by people.

We sought external dynamics that more complexly altered the torques to be overcome during movement. People have previously shown the ability to learn external viscous environments (Lackner and DiZio 1994; Shadmehr and Mussa-Ivaldi 1994); this skill learning has been probed to investigate the transfer of motor memories over speeds (Goodbody and Wolpert 1998), movement direction (Gandolfo et at 1996), and posture (Shadmehr and Moussavi 2000). We were motivated by this ready adaptation to viscous forces to consider optimal behaviors in the simplest of previously used environments, viscosity that pushes perpendicular to the velocity of the hand. The additional dynamics of the viscosity led to optimal behaviors that clearly differentiated MTC and MEPV from the straight trajectories of MJ.

The MTC solutions produced marked curvature, both on average and in specific directions (Fig. 3*B*, gray). The curvature of these trajectories was counterclockwise, away from the direction of the forces. Torques necessary to produce this curvature were less than those necessary to produce the straight MJ trajectory. This is reasonable because the curl nature of the viscosity pushed equally hard regardless of the initial direction of movement. A movement beginning counterclockwise to the target can ride the forces toward the target rather than fully resisting the forces throughout movement. The MEPV solutions (Fig. 3*B*, black) were even more curved toward the counterclockwise. This result seems reasonable as activation depends on the second derivative of torque and therefore MEPV seeks an even smoother torque generation than MTC.

Underexplored is the precise physiological implementation of signal-dependent noise. What is the physiological sampling rate of this noise? Is it on the millisecond time scale of neuronal noise? Is it much slower, perhaps on the order of the single- to double-digit hertz that underlie the phasing of EMG? We investigated these possibilities by implementing signal-dependent noise in our canonical model that updated at 1 kHz and then a different model that low-pass filtered that noise at 6 Hz. The slower frequency would roughly correspond to the occurrence frequency of each phase of a classical triphasic activation pattern underlying simple movements (Basmajian and DeLuca 1985; Thoroughman and Shadmehr 1999). We see very similar curvature across noise update frequencies and therefore do not see sensitivity to the particular temporal dependence of noise generation.

### Human behavior

Trajectories generated by people reaching in novel viscous environments quickly converse onto trajectories resembling original, unperturbed movements (Lackner and DiZio 1994; Shadmehr and Mussa-Ivaldi 1994). Here we examined behavior as subjects adapt, over hundreds of movements, to viscous forces that push perpendicular to hand velocity direction (Thoroughman and Shadmehr 2000). Performance at asymptote actually slightly overcompensated for clockwise forces, as average curvature was 2% counterclockwise. Although this curvature was in the direction of the MTC and MEPV trajectories, the magnitude of the curvature in the simulations was 8–10 times greater. If people plan these movements based on one cost, this suggests that this costs depends much more strongly on kinematic factors than dynamic factors or the variance of the endpoint.

### Humans and optimal control

The behavior of subjects after tens of minutes of training may not indicate optimal principles underlying motor planning. The proponents of minimal variance planning have suggested that the establishment of optimal behavior may take many times more movements, over several training sessions or even months, as the brain needs to ascertain the variance of a chosen behavior rather than just its mean. Furthermore, our 24 subjects began training with control of the passive dynamics of the manipulandum. Subject therefore could have retained recent kinematic experience as a goal of movement in spite of countervailing optimal control considerations.

Three pieces of evidence argue against this consideration and bolster the notions that our subjects can properly express optimal planning and that additional training would not alter behavior. First, other subjects have received more extensive training in viscous environments over several days (Donchin and Shadmehr 2004). Subjects that consistently trained in the curl field, without catch trials, did not vary their overall curvature from 2% counterclockwise throughout 4 days of training. Although this evidence does not consider training that occurs over weeks or months, the consolidation of motor skills allowed by hours of time (Brashers-Krug et al. 1996) and REM sleep (Karni et al. 1994) does not lead to trajectories that approach the curvatures of MTC or MEPV.

Second, subjects have demonstrated the ability to quickly adapt to novel dynamic demands by altering their trajectories from straight hand paths. A group of subjects recently experienced a mass-spring system simulated by a robotic manipulandum (Dingwell et al. 2004). With certain combinations of masses and spring constants, subjects quickly learned to move their hands in markedly biphasic speed profiles, very unlike trajectories generated by moving the manipulandum alone. Dingwell et al. hypothesized that this deviation aimed to optimize the smoothness of object trajectory; a simulation of this optimization indeed mimicked human behavior. Humans were capable of quickly and flexibly adopting a wide variety of movement patterns, dependent on movement dynamics and in pursuit of optimality other than hand path smoothness.

The last piece of evidence arose from one of our own subjects. Although the average across subjects indicated only subtle curvature, 1 of our 24 subjects exhibited more marked curvature (Fig. 3*D*, blue). This particular subject generated curvature that outlaid the rest of the group; his movements on average 0.11 and maximally 0.20. The MTC and MEPV curvature was only twice as large than his curvature. The performance of this subject indicates that people can adopt strategies markedly different from their recent performance in the baseline task. The remaining 23 subjects, however, show that the central tendency of human performance is motivated by kinematic optimization much more strongly than optimization of dynamics or endpoint variance.

### Possible combination of costs

The three optimal control hypotheses considered here all presume that motor planning depends on a single cost. An alternate hypothesis is that the nervous system considers a multitude of costs (Winters 2000) to plan movements and that a chosen trajectory either minimizes some set combination of costs or satisfies inequalities that determine acceptable levels of each cost. We do not know how the nervous system simultaneously considers multiple costs or could optimize a composite cost, but here we discuss an exploratory example.

The three candidate costs have different units so may not be directly summed; in this exploratory example, we chose to divide each cost by the minimum value found by our canonical model, therefore nondimensionalizing each cost and enabling summation. We then explored two linear domains of combined cost: one that jointly considered jerk and torque change and one that jointed considered jerk and endpoint variance. In each case, we explored the cost space by sliding a parameter λ, which determined the fractional inclusion of jerk (= λ) and the other cost (= 1− λ); for instance, in the jerk-torque change exploration, when λ = 0.2, the composite cost equaled 0.2 * jerk + 0.8 * torque change. We explored these cost spaces in simulated movements toward 0°, the direction of maximal curvature in the MTC and MEPV predicted trajectories. For each λ we found the trajectory that minimized the composite cost, then calculated the individual costs incurred and the maximum curvature in each found trajectory.

In these two cost domains (Fig. 5, *A* and *B*), we found that jerk rapidly rose as the composite cost disfavored jerk but that the other cost (TC or EPV) rose linearly with shallow slopes as the composite cost favored jerk. The combination of these curves led to curvature of optimal trajectories that markedly dropped as the cost favored jerk (as λ increased). If we construct 95% confidence intervals in these cost domains that correspond to our human performance (Table 1), we find surprisingly low and broad intervals; λ can fall within 0.32–0.54 in the jerk-TC domain and within 0.39–0.65 in the jerk-EPV domain.

We stress that this is an exploratory example; different assumptions of how costs are combined would lead to very different results. In these cost domains, however, we draw two conclusions. First, the rising of TC and EPV costs are surprisingly shallow and therefore may not generate strong gradients to drive optimization that determines optimal behavior. Jerk, on the other hand, has a strong gradient and is a good candidate to drive optimization. Second, the strong difference in trajectories predicted by optimizing single costs belies the broad and robust consideration of TC or EPV that could be included in composite costs that, when optimized, mimicked human behavior.

We conclude that human subjects, in this task, while newly experiencing viscous forces, exhibit behavior that cannot be mimicked by the MTC or the MEPV hypotheses. We propose that further exploration of the multitude of motor tasks and dynamic demands encountered by people constitutes a fertile ground for careful characterization of how competing costs influence motor planning.

## APPENDIX

### Inverse dynamics

We used an equation of motion to generate joint torques from desired trajectories (A1) Subscripts s and e represent shoulder and elbow; subscripts *h* and *r* refer for human and robot. The composite inertial and centripetal/Coriolis matrices are (Bhushan and Shadmehr 1999) (A2) where the inertial and centripetal/Coriolis matrices of the arm are (A3) with The inertial *D*_{r} and centripetal/Coriolis *C*_{r} matrices of the robot resemble the ones in *Eq. A3*, having *p*_{1}, *p*_{2}, ṗ_{1}, ṗ_{2}, p̈_{1}, p̈_{2} and the masses and lengths of the robot's links instead (Fig. 1).

Given *x*, *y*, ẋ, ẏ, ẍ, and ÿ of an initial motor plan in Cartesian coordinates, we can compute *q*_{s,} *q*_{e}, q̇_{s}, q̇_{e}, q̈_{s}, and q̈_{e} (A4) where *l*_{1} is the length of the upper arm and *l*_{2} is the length of the forearm. Then (A5) (A6) The Jacobian of the human or robotic arm is (A7) where the *l*'s are replaced by *r*'s and the human joint angles by robot joint angles for the robot Jacobian.

The lengths and masses of the arm in the simulations were determined from height, *h* = 168 cm and weight, *w* = 59 kg. The anthropometric relations (Winter 1990) are The kinematic and dynamic parameters of the robot (Bhushan and Shadmehr 1999) are (A8) where *r*_{1} and *r*_{2} are the lengths of the proximal and distal links.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2007 by the American Physiological Society