Abstract
Few computational models have addressed the spatiotemporal features of unconstrained threedimensional (3D) arm motion. Empirical observations made on hand paths, speed profiles, and arm postures during pointtopoint movements led to the assumption that hand path and arm posture are independent of movement speed, suggesting that the geometric and temporal properties of movements are decoupled. In this study, we present a computational model of 3D movements for an arm with four degrees of freedom based on the assumption that optimization principles are separately applied at the geometric and temporal levels of control. Geometric properties (path and posture) are defined in terms of geodesic paths with respect to the kinetic energy metric in the Riemannian configuration space. Accordingly, a geodesic path can be generated with less muscular effort than on any other, nongeodesic path, because the sum of all configurationspeeddependent torques vanishes. The temporal properties of the movement (speed) are determined in task space by minimizing the squared jerk along the selected endeffector path. The integration of both planning levels into a single spatiotemporal representation simplifies the control of arm dynamics along geodesic paths and results in movements with near minimal torque change and minimal peak value of kinetic energy. Thus, the application of Riemannian geometry allows for a reconciliation of computational models previously proposed for the description of arm movements. We suggest that geodesics are an emergent property of the motor system through the exploration of dynamical space. Our data validated the predictions for joint trajectories, hand paths, final postures, speed profiles, and driving torques.
 forward control strategies
 pointtopoint arm movements
 geodesics
 minimal effort
 minimum jerk
 minimum torque change
Introduction
A pointing movement toward a target in threedimensional (3D) space defines a highly redundant task at the geometric, kinematic, and dynamic levels of control. An infinite number of possible hand paths can be selected for moving the hand to the target, and an infinite set of possible arm postures may be attained for every fixed hand location in 3D space. In particular, many arm configurations may be adopted at the end of the movement. Moreover, different speed profiles may be chosen along a given hand path, and diverse patterns of muscle activation may generate the same driving torques that ultimately move the hand toward the target.
The rules that the CNS applies to the control of movements are, in fact, poorly understood. One major difficulty results from the fact that the planning and control strategies cannot be directly accessed, and only some kinematic and dynamic movement properties can be measured under well defined experimental conditions. Invariant kinematic and dynamic properties observed during movement have provided some important insights into possible motor control strategies used by the CNS but led to a large number of different and apparently incompatible models (Hermens and Gielen, 2004).
Most existing models are based on the assumption that movements are planned before their execution. For example, it has been proposed that arm movements are planned in terms of hand Cartesian coordinates by maximizing motion smoothness (Flash and Hogan, 1985) or in terms of intrinsic joint coordinates by minimizing the squared change of joint torques (Uno et al., 1989; Nakano et al., 1999; Wada et al., 2001) and by minimizing peak value of kinetic energy (minimum peak work) (Soechting et al., 1995). In addition, stochastic models have been proposed assuming that the inherent noise in the motor system is minimized (Harris and Wolpert, 1998) or that the noise is optimally distributed among different degrees of freedom (Todorov and Jordan, 2002).
Experimental observations suggest that hand paths and arm postures are invariant with respect to the scaling of movement speed (Flash and Hollerbach, 1982; Boessenkool et al., 1998; Nishikawa et al., 1999) and change in external forces (Atkeson and Hollerbach, 1985; Flanders et al., 2003). It has been shown that geometrically defined movement features can be acquired by the human motor system (Sosnik et al., 2004), before the specification of kinematic attributes such as movement speed (Torres and Zipser, 2002, 2004). It is thus conceivable that separate planning constraints may be imposed at the geometric and temporal levels. Both levels may subsequently be integrated into a complete spatiotemporal representation.
The present study aimed at developing a model of the motionplanning strategy underlying the generation of 3D pointing movements. The model is based on Riemannian geometry and accounts for hand paths, final arm postures, handspeed profiles, and driving torques. It results in arm movements along geodesic paths that require less muscular effort than on any other, nongeodesic path, while maximizing smoothness. Our model thus reconciles different existing computational approaches into a single framework.
Materials and Methods
Theoretical background
A possible interpretation for the experimentally observed invariance of arm posture and hand path with changes in speed is the decoupling of geometric and temporal aspects in the trajectoryplanning process of 3D arm movements. The objective of the computational model presented here consists of the prediction of movement patterns as resulting from a model that assumes that the geometric and temporal levels are independently planned, and the comparison of the predicted movements with experimentally observed trajectories. In the first stage, the geometrical properties of the movement as expressed by the jointangular path in configuration space are determined. In the second stage, the temporal features of the movement in terms of hand speed are selected. In the third and final stage, the geometric and temporal features of the movement are integrated together into a unique spatiotemporal representation of the movement.
The approach for selecting the geometric properties of the movement is based on Riemannian geometry. Most readers will be familiar with Euclidean geometry, in which the squared distance between two nearby points that are separated by the vector dx = (dx_{1}, dx_{2}, dx_{3})^{T} is given by or written alternatively using a vector notation as
One also recalls that the shortest distance between two points in Euclidean space is a straight line. A train analogy is useful to further develop our basic model assumptions. We first observe that for a train, the speed is in principle independent of the shape of the specific railway track, similar to the decoupling of movement path from movement speed. We next consider the forces that act on the train in a hilly landscape. One might ask the following question: What is the optimal geometrical path between two stations A and B for laying a railway track in terms of forces that act on the train? This problem can be addressed by passing from Euclidean to Riemannian geometry.
In Riemannian geometry, the distance relationship is generalized with a metric tensor g(x) that encodes the geometry of the hilly landscape locally near the point x. The squared distance is given by Thus, Euclidean 3D space is a special case of a Riemannian manifold in which
The optimal path for the train is given by the straightest possible path through the hilly landscape, because for this path the centripetal forces acting on the train are minimized. The determination of the straightest path is a standard problem in Riemannian geometry and results in the computation of the geodesic path connecting station A with B.
In the following paragraphs, we use Riemannian geometry to model human pointing movements in space. Using the previous analogy, we identify the arm with the train and the hilly landscape with the Riemannian configuration space. Thus, we are interested in the answers to the following questions: What is the optimal geometric path of the arm between a given initial and final arm posture in configuration space in terms of forces acting on the arm? What might be the optimal speed for the arm to move along the chosen path?
Detailed mathematical analysis.
To further elaborate on these ideas, we have to provide some mathematical tools, which we present next. The configuration space of the arm, Q, defines a Riemannian manifold when endowed with a positive definite and symmetric metric tensor g with components g_{ij}, (i, j = 1, …, n). A point in configuration space is denoted by q with (local) coordinates q = (q^{1}, q^{2}, …, q^{n}) ∈ Q. Distances in arm configuration space are defined as in Equation 3, which can alternatively be written in terms of the coordinates as where we have used the letter σ for distance in the Riemannian manifold, and the letter s is preserved in the following for distance in Euclidean space (Eq. 1). We make use in the following of the summation convention stating that repeated lower and upper indices imply a summation from 1 to n; thus, Equation 5 can be rewritten as dσ^{2} = g_{ij}(q)dq^{i}dq^{j}. The metric tensor is symmetric, g_{ij} = g_{ji}, and assumed to be positive definite; i.e., it is g_{ij}q^{i}q^{j} > 0 for all q ≠ 0, and thus it is invertible. The inverse of the metric tensor g^{−1} has components (g^{−1})_{ij} = g^{ij} with In the Riemannian manifold, (Q, g), the length of a curve γ: [a, b] → Q is defined as the line integral where λ is an arbitrary parameter. Another quantity that is used in the following is the energy of a curve, which is defined as The name “energy” derives from the similarity of the integrand to the kinetic energy of a physical system. However, it is important to note that the energy of a curve is a geometrical quantity.
The curves that correspond to extremal paths in a Riemannian manifold are called geodesics. Geodesics are curves that are locally of minimal length. This means that any two points that are close enough are connected along the geodesics by the shortest possible path. It is not difficult to prove (do Carmo, 1992) that a curve of minimal length also minimizes the energy function (Eq. 8).
A necessary condition for an extremum of the energy function (Eq. 8), and thus for the length (Eq. 7), follows from the calculus of variation in the form of the Euler–Lagrange equations. The Euler–Lagrange equations applied to Equation 8 lead to the geodesic equation (i = 1, …, n) where the Christoffel symbols of the first kind are defined by We can rewrite Equation 9 using the “standard form” by multiplying it with the inverse of the metric to obtain where the Christoffel symbols of the second kind are defined by Γ_{jk}^{i} = g^{il}Γ_{jk}.
It is important to note that a geodesic curve is not only defined by the shape of its path, but also by its parameterization. It can be easily shown by using the geodesic Equation 11 that dσ/dλ = const, and thus, the parameter λ is linearly related to the arc length σ; i.e., it is λ = aσ + b with some constants a and b. A parameter λ with this property is called affine parameter, and thus geodesic curves are parameterized with respect to an affine parameter. In particular, the arc length itself defines an affine parameter (λ = σ) [the minimization of the energy function (Eq. 8) leads automatically to the parameterization of the geodesic path in terms of an affine parameter].
Of course, nothing prevents us from choosing any other, nonaffine parameterization (e.g., time t) along the path. However, such a choice implies that the geodesic Equation 11 is not satisfied. This can be shown if we define a new curve parameterization by the functional relation t = f(σ), where σ is the arc length in configuration space. Then the derivatives are related according to the chain rule by where a prime denotes differentiation with respect to σ. With the new parameterization, the geodesic equation in standard form becomes
The left side of Equation 13 defines the acceleration in curvilinear coordinates and the term on the right side can be interpreted as a generalized force. For t to be an affine parameter, f″ must vanish; i.e., σ and t must be linearly related or dσ/dt = const. In this case, the force term vanishes, and Equation 13 transforms again into the geodesic equation (Eq. 11). Because the expression dσ/dt defines a speed (if t defines time), we conclude that geodesic paths are forcefree paths of constant speed. Note that Equations 11 and 13 describe the same geodesic path but not the same geodesic curve.
To close this section, it is instructive to evaluate some of the above expressions in Euclidean geometry with coordinates x = (x_{1}, x_{2}, x_{3})^{T} and a metric defined by Equation 4. For example, the energy of the curve is given by The geodesic equation (Eq. 11) transforms to x″(λ) = 0, leading to the well known result that geodesic paths of Euclidean space are straight lines.
Implications for human trajectory formation.
How are these general mathematical considerations related to the model presented in the following sections? Our computational model assumes that the CNS adopts geodesic paths at the geometrical level as part of the trajectoryplanning process. It will be shown in the following that geodesic paths with respect to a suitably chosen metric in configuration space simplify the arm dynamics significantly. It is this link between geometry and dynamics that makes the use of Riemannian geometry so attractive. We then hypothesize that at the temporal level, the CNS selects a speed profile along the hand path in Euclidean task space, which in turn induces a nonconstant speed in configuration space (dσ/dt ≠ 0). The latter corresponds to a reparameterization of the geodesic paths in terms of a nonaffine parameter. This reparameterization does not change the shape of the geodesic path, but leads to a movement that is not forcefree, and therefore, torques at the joints are required to drive the arm toward its final configuration. We will show in the following that the torques are reduced along geodesic paths, because the sum of all configurationspeeddependent torques vanishes, whereas the remaining driving torques are to firstorder approximation linearly related to the hand acceleration. It is hypothesized in this paper that these paths are an emerging property of the system that seeks to reduce muscular efforts possibly based on proprioceptive feedback.
The next sections present the computational model. First, the forward and inverse kinematics of an arm with four degrees of freedom (DOFs) are presented as necessary prerequisites for the formulation of the computational model. This is followed by the details of the computational model and the description of the experimental protocol.
Forward kinematics
The human arm is approximated as a linkage of rigid bodies, and an arm configuration is parameterized by four joint angles q:= (θ, η, ζ, φ)^{T} ∈ Q, where we follow a parameterization of the arm configuration given in Soechting et al. (1995). Q denotes the configuration space. The first three angles describe the rotation around an ideally spherical shoulder joint, namely, the elevation angle θ, the azimuthal angle η, and the humeral angle ζ, whereas the flexion angle φ determines the rotation around the elbow joint (Fig. 1A). We neglected the three DOFs at the wrist, which is assumed to be fixated. The shoulder joint is located at the origin of the laboratory frame coordinate system (CS). Two fixed body frames, CS_{1} with axes x_{1}, y_{1}, z_{1} and CS_{2} with axes x_{2}, y_{2}, z_{2}, are attached to the upper arm and forearm segments, respectively, such that the zaxes are pointing along the longitudinal limb axes of the upper and forearm, and the x and yaxes are in the transverse directions. In the zero configuration [q = (0, 0, 0, 0)^{T}], the arm is fully extended in the direction of the (−z)axis, and the axis of the elbow joint is aligned with the xaxis. For the transformation of the arm from the zero configuration into an arbitrary configuration, the order of rotations has to be specified. In this work, an arm posture is defined as the sequence of the following rotations: (1) rotation around the x_{1}axis by angle π/2, (2) rotation around the y_{1}axis by angle η, (3) rotation around the x_{1}axis by angle θ − π/2, (4) rotation around the z_{1}axis by angle ζ, and (5) rotation around the y_{2}axis by angle φ, where the sign of the rotation angle is defined by the righthand rule. The forward kinematics map defines the elbow and hand locations in the laboratory frame CS as a function of the joint angles; i.e., x_{e} = F_{e}(q) for the elbow location and x_{h} = F_{h}(q) for the hand location, where x_{e} = (x_{e}, y_{e}, z_{e})^{T} and x_{h} = (x_{h}, y_{h}, z_{h})^{T}. For the chosen parameterization, we obtain and where l_{u} and l_{f} are the upper arm and forearm lengths, respectively. It should be noted that the immobilization of the wrist motion, as was done for the experimental data analyzed here, implies that the hand and the wrist follow the same kinematics as one rigid link. Finally, for a given hand path, x_{h} = x_{h}(λ), the change in hand path can be determined according to x′_{h}(λ) = J_{h}(q)q′(λ), where (J_{h})_{ij} = dF_{h},_{i}/dq^{j}, i = 1, 2, 3; j = 1, …, 4 is the hand Jacobian, and a prime denotes differentiation with respect to λ, where λ is a parameterization of the path. Similar relationships hold for the elbow path. The components of the elbow and hand Jacobian are given in Appendix A (the “hand Jacobian” defines the Jacobian of the whole arm system, whereas the “elbow Jacobian” defines the Jacobian of the upper arm segment).
Inverse kinematics
We determine next the inverse kinematic relation for the four DOF arm, which defines the jointangular vector as a function of the hand location. As a result of the excess of one DOF in configuration space (four DOFs in configuration space and three DOFs in external task space for an immobilized wrist), there exists no unique map from hand to configuration space. Indeed, for a fixed hand location, the arm can still rotate around an axis going through the shoulder and the hand, implying that the elbow location is constrained to move on a circle around this axis (Hollerbach, 1985). Changes in the elbow position along the circle do not influence the hand location, and therefore, such changes do not contribute to the task goal. If we denote the rotation angle of the elbow around this axis with α ∈ [0, 2π), we can parameterize the elbow locations for a fixed hand location x_{h} as x_{e} = x_{e}(α; x_{h}) (Fig. 1B). An explicit expression for the elbow location in terms of the angle parameter α and the hand location x_{h} is provided in Appendix B. Once the elbow and hand locations are fixed, the whole arm configuration is determined; i.e., there exists a relation of the form q = q(x_{e}, x_{h}). An explicit calculation assuming that x_{e} × x_{h} ≠ 0 and z_{e} < l_{u} leads to where we defined the function atan2(a, b):= atan() − sign(a)[1 − sign(b)](). The inverse kinematic relation follows then by inserting the expression for the elbow locations, x_{e} = x_{e}(α; x_{h}), in the relationship expressing the jointangular vector in terms of elbow and hand location, q(x_{e}, x_{h}).Thus, the inverse kinematic relation has the form q = q(α, x_{h}). Note that the four components of the jointangular vector are determined by the rotation angle α and by the three components of the hand location.
Obviously, not all rotation angles α ∈ [0, 2π) lead to a realizable arm posture. The subset of rotation angles that lead to postures inside the biomechanically admissible joint range are determined by using biomechanical jointrange models. We assume that only joint angles that satisfy the following inequalities lead to admissible arm configurations: where the external and internal humeral rotation of the upper arm, ζ_{ext} and ζ_{int}, describe the maximal clockwise and counterclockwise rotations, respectively, around the upperarm axis as a function of upperarm axis orientation (when viewed from the subject). The external and internal humeral rotations are derived from experimental data and are approximated by the curved surfaces shown in Figure 2 (Wang et al., 1998).
For each hand location in space, x_{h}, the interval of possible rotation angles, I(x_{h}) ⊃ α, defines the arm posture constraints. The corresponding jointangular vectors associated with realistic arm postures at hand location x_{h} are thus specified by the oneparameter family of vectors q(α) = q(x_{h}, α) with α ∈ I(x_{h}).
Description of the computational model
The implementation of the geometric and temporal levels in the trajectory planning process and the integration into a spatiotemporal representation are described in the following sections.
Geometric level: prediction of the jointangular path in configuration space.
The configuration space of the arm defines a Riemannian manifold when assigning a suitable (positive definite and symmetric) metric tensor in configuration space. A natural candidate for a metric in configuration space is defined by the (in general, nonEuclidean) kinetic energy metric M (or manipulator inertia matrix in the robotic literature); i.e., we set g_{ij} = M_{ij}, i, j = 1, …, 4 (Biess et al., 2001). The components of the kinetic energy metric are given by where K is the kinetic energy of the arm and q^{i}, i = 1, …, 4, denote the components of the jointangular vector q = (θ, η, ζ, φ)^{T}. The explicit form of the kinetic energy and the components of the kinetic energy metric for a four DOF arm are provided in Appendix C.
The kinetic energy metric depends on the current arm configuration and represents physically the instantaneous composite mass distribution of the whole arm linkage in the current arm configuration (Asada and Slotine, 1986). The use of this metric is motivated by studies of Soechting et al. (1995) and Flanders et al. (2003) that showed the significance of inertial properties of the whole arm for the selection of the arm configuration in 3D space. The observed variability of handpath curvature, which depends on the location in the task space, makes the kinetic energy metric attractive for the description of 3D pointtopoint movements. In addition, the kinetic energy metric is closely related to the dynamic equations of motion of the arm, as will be shown in Results. The line element associated with the kinetic energy metric is, according to Equation 5, We can then compute the geodesic path for the kinetic energy metric between an initial and final configuration, q_{0} and q_{f}, respectively. In this work, we define the affine parameter, λ, of the geodesic path as follows: if we denote with Σ the total arc length of the path in configuration space between the initial and final configuration, we define the (dimensionless) affine parameter by the normalized arc length λ = σ/Σ, which then takes values between 0 and 1. The geodesic equation for the kinetic energy metric is given as a result of Equation 9 by (g_{ij} = M_{ij}) where q′ = dq/dλ and the components of, the Coriolis matrix, C = (C_{ij}), i, j = 1, …, 4, are defined by
The specification of suitable boundary conditions is required for the solution of the geodesic equation (Eq. 31). We set For this choice of boundary conditions, the complete initial arm configuration, q_{0}, and the final hand location, x_{h} = x_{f}, are specified. However, because of kinematic redundancy, the given final hand location corresponds to a oneparameter family of possible final arm configurations, namely, all the accessible arm configurations that are generated by a rotation around the axis connecting the shoulder joint and the final hand location. The final arm posture is thus not uniquely determined. Instead, it is an outcome of the optimization procedure described below.
We remark that the computation of the solution of the geodesic equation (Eq. 31), subject to the set of boundary conditions in Equations 33 and 34, is not a trivial problem. Mathematically, it defines a nonlinear twopoint boundary value problem. A numerical solution for the computation of the geodesics with a high degree of accuracy was developed and presented by Biess et al. (2006), and we refer the reader to this work for additional details.
We describe next the computational steps to determine the geometric prediction of the model in terms of the optimal jointangular path in configuration space (Fig. 3). The configuration space Q is visualized as a twodimensional surface. Each point in configuration space corresponds to an arm posture. The input to the model is given by the initial arm configuration, q_{0}, and the final target location x_{f}, where the latter corresponds to a oneparameter family of final arm configurations. To extract this set of final arm postures, we compute first the circle C in task space, which contains the set of final elbow locations. The set of final elbow locations determines together with the fixed final hand location the set of possible final arm configurations. In configuration space, this set can be represented as a onedimensional manifold. Next, we determine the set of admissible arm postures. For this purpose, the final elbow positions on the circle C are discretized (by variation of the angle α in steps of 3°), and the final elbow locations inside (●) and outside (○) of the biomechanically admissible joint range are derived from biomechanical jointrange models. These sets correspond in configuration space to subsets of the onedimensional manifold. In the final step, we compute for each fixed value of the angle α ∈ I(x_{h}) the geodesic path between the initial arm configuration, q_{0}, and the final arm configuration, q_{f}(α) = q(α; x_{f}), α ∈ I(x_{f}). The solution of Equation 31 subject to Equations 33 and 34 thus defines a oneparameter family of geodesics in configuration space, γ_{α}: [0, 1] → Q, that all start at the same point but end at different points in configuration space. In extrinsic task space, this corresponds to a fixed initial arm posture, whereas many final postures may be compatible with the given final hand location. Among all geodesics that are compatible with the task constraints, we assume that the CNS selects the geodesic path with the minimal length in configuration space; i.e., we select the optimal path in joint configuration space as the geodesic path with the property where the energy of the curve is given according to Equation 8 by and α = α* is the rotation angle that minimizes the energy of the curve. Σ(α) denotes the arc length of path γ_{α} in configuration space. The optimal jointangular path is thus determined by the geodesic path γ_{α*}: q*(λ) = q(λ, α*), and the optimal hand path in task space follows from the forward kinematics map. The temporal properties of the movement are determined by ascribing a speed profile along the hand path, as presented in the next section.
Temporal level: prediction of the speed profile in Euclidean task space.
The timing of the movement is derived in task space, where distance is defined according to the Euclidean metric (Eq. 2). It is assumed that the speed profile along the hand path is determined by minimization of the squared jerk of the arc length s of the hand path integrated over the movement time: subjected to the boundary conditions where T is the measured total movement duration, S is the total Euclidean arc length of the hand path, and a dot denotes differentiation with respect to time. Note that Equations 37 and 38 correspond to a modified minimumjerk (MJ) model that is expressed in terms of the arc length along the hand path rather than in terms of the endeffector coordinates as in its original formulation. Remember that the solution of the original minimumjerk model resulted in straight hand paths with a bellshaped velocity profile. In the modified minimumjerk model, we only determine the time course (speed) along an a priori defined path that can have any smooth shape.
The modified minimumjerk model is motivated by the observation that speed profiles of pointtopoint movements in 3D are roughly bell shaped (Atkeson and Hollerbach, 1985). This choice will be also motivated retrospectively when discussing the dynamical properties of the resulting movement. The optimal solution for position and speed is given by (Flash and Hogan, 1985) The total movement time is an input parameter taken from the measurements, whereas the total arc length of the hand path is not known a priori but can be computed from the jointangular paths determined at the geometrical level. The arc length of the hand path in Euclidean task space follows from Equation 2 as which leads to where we expressed the hand vector in terms of the jointangular vector using the hand Jacobian. The total arc length of the hand path in task space follows then as S = s(1). It should be noted that the arc length of the hand path as a function of path parameter is a monotonically increasing function (ds/dλ > 0) and thus can be inverted, leading to a function λ = λ(s).
The relation between the speed in configuration space, λ̇, and the speed in task space, ṡ, is derived next. From Equation 41, we get where a dot and a prime denote differentiation with respect to time and arc length s of the hand path, respectively. For later purposes, we also derive the acceleration in configuration space, which follows by derivation of Equation 43 as Note that the kinetic energy of the arm is according to Equations 30 and 43 given by (σ = λ/Σ) Equation 45 is used in Appendix D for the derivation of the minimumpeak kinetic energy (minimum peak work) model, which follows as an outcome of our computational model.
Spatiotemporal level: integration of the geometric and temporal levels.
The spatiotemporal representation of the movement can be constructed from the selected motion patterns at the geometric and temporal levels by two successive transformations. First, the optimal jointangular path, q*(λ), is reparameterized with respect to the arc length of the hand path; i.e., q*(s) = q*(λ(s)). Then, the optimal time profile is imposed along the hand path, resulting in the optimal jointangular trajectory, q*(t) = q*(s*(t)), which defines the optimal movement kinematics. For brevity, we refer in the following to the model with these temporal and geometrical features as the geodesic (GEO) model. The implications of the GEO model for the arm dynamics are presented in Results.
Experiments
Subjects.
Four righthanded male volunteers (age range, 18–32 years) completed a series of natural unconstrained arm movements of two different types within a single session. None of the participants reported having any clinical symptoms or any history of motor, sensory, or neurological disorders. All participants gave their signed consent to participate in the experiment, as requested by the institutional ethics committee.
Procedure.
Subjects sat in front of a projection screen (a parafrontal plane relative to their torso) with the shoulder at a distance of 1 m. They were strapped to a chair by appropriate belts that minimized shoulder displacements and fixated the torso throughout the experiment. The height of the chair was adjusted such that the right shoulder was aligned with a fixed reference point in space. The subjects were naive to the purpose of the experiment. They were instructed to point toward visual targets that were randomly presented in different locations of the 3D task space. All experiments were performed under dimmed lighting conditions such that individuals were able to see their arm. Targets consisted of 5 cm balls that were projected on the screen. For this purpose, a 3D virtual reality was generated using stereovision. For additional details, we refer to Liebermann et al. (2006). The 3D positions of infrared emitting diodes (IREDs) were recorded at a rate of 100 Hz with an accuracy <1 mm using two motion tracking cameras (Optotrak; Northern Digital, Waterloo, Ontario, Canada). The IREDs were attached to the arm segments via exoskeletal metal frames. The center of the frame placed over the acromion (a triangle of 10.7 cm base and 11 cm height) was the assumed center of shoulder rotation. The center of the exoskeletal frame placed on the distal end of the upper arm (18.5 cm height, 16.2 cm width) was used to measure the elbow joint location. A third frame (11.4 cm height, 11 cm width) was centered and attached to the wrist. The latter was braced to eliminate any rotations of the radiocarpal joint (i.e., supination–pronation movements were only possible by rotating the forearm about the radioulnar joint).
Experimental protocol.
Data were collected during the performance of radial and frontoparallel movements. Radial movements were defined as the hand motion in the inward and outward directions relative to the body's longitudinal axis. In the radial condition, targets were presented at fixed distances of 0.6 and 0.8 m from the location of the shoulder center, at different heights. Pointing always started from a central position (0.3 m frontally and at 0 m horizontally and vertically relative to the shoulder). A set of radial movements consisted of 78 trials from the initial to the final targets and vice versa (forward and backward movements). Each set of movements was repeated five times, resulting in a total number of 78 × 2 × 5 = 780 movement trials.
Frontal plane movements were defined as movements between targets lying on frontoparallel planes relative to the body. Targets were presented in four frontoparallel planes at fixed distances of 0.30, 0.45, 0.60, and 0.75 m from the shoulder center. Each plane contained four initial targets and 12 final targets presented at different heights (−0.4, 0, and 0.4 m relative to the shoulder). It is important to note that the end effector was not constrained to move in the virtual plane between an initial and a final target. A set was comprised of 24 pointing movements from a fixed initial target position to several final targets in the same frontal plane and their reversals. The initial target within one of the four planes was changed four times, and each set was repeated five times, resulting in a total number of 5 × 16 × 24 = 1920 movements. In both movement conditions, the subjects could choose freely the initial arm configuration at the initial target.
Data analysis.
During the data collection process, markers were not always visible to the cameras, and data were lost. However, only if more than onethird of the total number of samples obtained for the IRED markers were lost during a movement was the trial excluded. In all other cases, a data recovery process was initiated as follows. Geometric and temporal information was used for recovering data. If one marker of an exoskeleton at a certain time frame was missing, it could be simply restored by geometrical considerations based on the rigidity of the exoskeletons. If marker signals were lost over a whole interval of frames, where the maximal allowed contiguous missing segment was set to 15 frames = 150 ms, the missing markers were reconstructed by extrapolating from all visible markers before and after the missing time frames. Raw displacement data obtained from the markers were filtered (Butterworth, second order). The cutoff frequency was chosen depending on whether position, velocity, or acceleration data were generated (position: 6.0 Hz; velocity: 5.5 Hz; acceleration: 5.0 Hz) (Giakas and Baltzopoulos, 1997). The movement onset and offset times are determined by optimally fitting (leastsquare error) a superposition of minimumjerk speed profiles to the absolute hand velocity profile (Lee et al., 1997). Not more than two minimumjerk submovements were needed to extract the movement onset and offset times. This procedure was only used to reconstruct the segments of the experimental speed profiles near the beginning and the end of the movement (approximately the segments of the speed profiles for which v < 0.15v_{max}) and not to fit a minimumjerk speed profile to the data.
The estimation of the elbow locations from the marker positions of the exoskeletons, which are attached to the limbs, is an essential part of the data analysis. Skin movements and random measurement noise lead to errors in the measured joint locations. In particular, errors occur in the localization of the elbow joint because of the difficulties in attaching the exoskeleton such that no relative motions between arm and grid occur while moving. However, a reliable elbow position is required for the determination of the joint angles and the length of the limbs. The estimation of the elbow position was composed of three steps. First, the marker coordinates were transformed to a coordinate system that moved with the upper arm. Second, an optimal estimate of the elbow location in the upper arm coordinate system was determined (Gamage and Lasenby, 2002). In this coordinate system, only the forearm moves, and the elbow location is a fixed point that can be determined by geometric methods. The estimated elbow locations in the shoulderfixed coordinate system followed then by backtransformation.
Results
The results of the Riemannian approach were compared with the predictions of the Euclidean approach, where a Euclidean metric is assumed in task and configuration space. Thus, the metric M(q) is replaced by the identity matrix. The geodesic paths in configuration space are then simply straight lines given by This path is the optimal solution to the squaredjoint derivative cost (SJD), The temporal prediction for the SJD model was assumed to result from the minimumjerk model for the arc length, as specified in Equations 37 and 38.
In addition, the results of the two previous models were compared with the predictions of the minimum torquechange (MTC) model. This model minimizes the squared change of joint torques integrated over the total movement time (Uno et al., 1989). Boundary conditions for the MTC model were set according to Equations 33 and 34, supplemented by zero jointangular velocities and accelerations at the beginning and end of the movement. As for the geodesic model, the final arm posture was not predefined, but rather resulted as an outcome of the optimization. For the solution of the MTC model, an optimization method based on the parameterization of jointangular trajectories was used (Biess et al., 2006).
Simulated and observed data are compared in this section. Results are presented in four subsections that addressed predictions at the geometric, temporal, spatiotemporal, and dynamic levels.
Geometric level: hand paths and final postures
Measured hand paths were persistently curved depending on the location in the workspace. However, the curvature was small for all types of movements. The curvature of the hand paths was assessed by evaluating a global curvature index defined as CI = S/L, where S denotes the total arc length of the hand path and L the Euclidean distance between the initial and the final targets. For straight hand paths, the curvature index is CI = 1, whereas for a semicircular path, the index is CI = π/2 ≈ 1.57. The last columns of Figures 9 and 10 show the distributions of the curvature index for radial and frontoparallel movements, respectively. As can be seen from these figures, most of the movements result in quasistraight hand paths with a curvature index CI ≤ 1.05.
Typical examples of the measured and predicted hand paths for the three models (GEO, SJD, and MTC) in the radial condition are shown in Figure 4 and in the frontal plane condition in Figure 5. Each row shows the projections of the three dimensional hand path of one movement into the xy, xz, and yzplanes. These randomly chosen examples from the data of our four subjects show that the MTC model consistently deviates from the observed paths in all three planar projections of the 3D hand movements. This was observed regardless of the movement condition or direction. In the same vein, it may be observed that the GEO and SJD models are close to each other. For the evaluation of the path predictions, a (global) handpath deviation index (HPDI) based on the maximal value of the minimal Euclidean distances between measured and predicted paths was calculated. Such an index was used to quantify differences such as those observed in Figures 4 and 5. The measured path was divided into M = 40 equidistant segments, which defined M − 1 inner points. The minimal Euclidean distances from the inner points to the predicted path, R_{i}, i = 1, …, M − 1, were determined, and the ratio of the maximal value R = max_{i}_{=1,… M}_{−1}{R_{i}} and the distance between the initial and final targets L was defined as the handpath deviation index, HPDI = R/L (Fig. 6A).
In all presented cases, the path predictions of the GEO model were closer to the measured hand paths, as indicated by the mean HPDI scores presented in Table 1, than the path predictions derived from the two other models. The paths resulting from the GEO and SJD models were similar, whereas the latter led to slightly larger HPDI values. In contrast, the hand paths resulting from the MTC model were strongly curved in most trials, and thus, large HPDI values were obtained for the MTC model.
In the present experiments, the arm was free to adopt any configuration among a large number of possibilities available during pointing to visual targets, in particular, at the final target. However, reproducible postural patterns were commonly observed within a limited set. The present model was designed to predict those arm postures that bring the hand to its final position by the shortest geodesics that connects the initial and the final arm posture in configuration space. The predicted versus the measured final joint angles of the shoulder are shown in Figure 7 for the radial movements and in Figure 8 for movements in a frontal parallel plane, together with the R^{2} values for the GEO, SJD, and MTC models. The figures show that the largest deviations from the perfect fit, where all data points lie on the diagonal, occur for the MTC model, whereas the GEO and the SJD model lead to similar and less spreadout distributions around the regression line. The analysis of the R^{2} values resulted in the highest score for the GEO model, although the torsional angles ζ were not always predicted successfully. For example, in the case of subject 1 in the radial condition and subject 3 in the frontalplane condition, Figures 7 and 8 show respectively for each of these subjects and conditions that the torsion angle was not predicted very well. The origin of this large discrepancy for these subjects remained unclear.
To further evaluate the predicted final arm postures, a relative error measure of the form of the final posture deviation index (FPDI) was determined as follows. At a given target location, the arm can rotate around an axis going through the shoulder and the final hand location. The rotation angle around this axis is given by the angle α (Fig. 1B). The measured and predicted final arm postures can thus be defined by the angles α_{exp} and α_{pre}, respectively. The deviation index for each final arm posture follows then as the ratio of the absolute difference between the predicted and measured angles, Δα = α_{exp} − α_{pre}, to the total anatomically accessible range, Δα_{tot} = I(α; x_{f}) at the final hand location x_{h} = x_{f}. Note that the total anatomically accessible angular range depends on the hand location and can be estimated from biomechanical jointrange models (Fig. 6B). As shown in Table 1, the FPDI measures were smallest for the GEO model, followed by the predictions of the SJD cost. In comparison, large values were obtained for the MTC model.
Temporal level: speed profiles
We first analyzed the characteristics of the observed tangential handspeed profiles and investigated whether our modeling assumptions on the temporal level were justified. For this purpose, the normalized speed profiles, ṽ(τ) = v(τ)/〈v〉, were analyzed, where the average velocity is given by 〈v〉 = S/T and τ = t/T defines normalized time. The normalization guarantees that all speed profiles are independent of the distance traveled and satisfy the following relationship: ∫_{0}^{1}ṽ(τ)dτ = 1. To assess deviations from the symmetric speed profile, we defined an asymmetry index as the ratio (S_{l} − S_{r})/S, where S_{l} is the distance traveled up to peak time of velocity (acceleration phase) and S_{r} is the distance traveled from peak time of velocity until the end of the movement (deceleration phase). S denotes the total distance traveled. The temporal evolution of hand paths according to the minimumjerk description assumes that the peak amplitude of the normalized speed profile should invariably reach a value of 1.875 at time τ = 0.5. In addition, the bellshaped minimumjerk speed profile has an asymmetry index of zero.
The distributions around the mean values (by subjects and movement types) for peak amplitude, peak time, and asymmetry index are given in Figure 9 for radial movements and in Figure 10 for frontoparallel movements. These distributions suggest that the minimumjerk model provides a reasonable fit to the mean values of amplitude, peak time, and asymmetry factor, although the amplitude for 3D arm movements may be better predicted by the minimumsnap model (snap = fourth derivative of handposition vector), which would lead to a value of 2.186 instead of 1.875 resulting from the minimumjerk model (Richardson and Flash, 2002). Measured acceleration times were on average slightly shorter than predicted by the minimumjerk model, as indicated by the peak times of the handspeed profiles. The mean asymmetry factor of the speed profiles was small and negative, showing that the distance traveled during the deceleration phase was slightly longer than the distance covered during the acceleration phase. Overall, the MJ model formulated for the arc length of the hand path is a good approximation, although it cannot explain all temporal features of the measured hand movements.
Randomly selected examples of predicted speed profiles resulting from the MJ model along the hand path and the MTC model, superimposed on measured speed profiles, are shown in Figure 11 for the radial movement type and in Figure 12 for movements on the frontal plane. The predicted MJ profiles are in good agreement with the experimental data, although the predicted peak amplitude slightly undershoots the experimental value, and the observed small deviations from the symmetric bellshaped form cannot be explained by the MJ model. In contrast, the peak amplitudes of the MTC speed profiles are in general too small and show double peaks, in disagreement with the observed data.
The speed profiles were further examined by using a global error measure as in Nakano et al. (1999). The error measure was based on the area that the normalized speed profile encloses with the time axis. When comparing the predicted and measured normalized speed profiles, the compounded area surrounded by the two profiles was divided into a common and a noncommon portion. The ratio of the noncommon area to the whole area defined a speed deviation index (SDI) for each profile (Fig. 6C). The SDI values for the different models are listed in Table 1. Note the large difference for the SDI measures between the GEO and the MTC model and the fact that the SDI values for the GEO and the SJD model are identical.
The results of the mixeddesign ANOVAs (3 models × 2 conditions of movement) with repeated measures on the last factor were performed using the different error measures as the dependent variables. The results showed that regardless of the error variable used to assess the modelobserved differences, the frontal and radial movement conditions did not significantly differ from each other. The interactions shown in Figure 13 did not achieve statistical significance either. However, a major effect of models was found using path, speed, and posture error measures (p < 0.001 in all cases). Post hoc pairwise comparison showed that the model effects were caused by the larger discrepancies when the MTC model was assessed regardless of the movement type. The post hoc tests also showed that the MTC model differed significantly (p < 0.001) from the GEO and SJD models with respect to the path, speed, and posture error measures. However, the latter two were similar in terms of path error and identical in speed error measures. The GEO and SJD models differed from each other in terms of the posture error (post hoc pairwise comparison shows a borderline value p = 0.0537), with an advantage for the former model, suggesting that the GEO model was the best to predict the posture of the arm. Finally, it is worthwhile mentioning that within the narrow error margins measured for path, speed, and posture, subjects differed from each other. However, all subjects showed the same relative differences with respect to the predictions of the GEO, SDJ, and MTC models. Therefore, the present results suggest that the MTC model is not successful in predicting path, speed, or posture during pointtopoint hand movements in 3D. The GEO model proposed in the current study yielded the best results.
Spatiotemporal level: jointangular positions
The model presented here assumes that the two preplanned aspects of the movement are integrated at some stage into one spatiotemporal representation. The jointangular trajectories are an outcome of this integration process and are computed as described at the end of the modeldescription section.
Figures 14 and 15 show several examples selected at random of predicted and measured jointangular trajectories for the three shoulder angles and the flexion angle at the elbow joint. Note that the arm configuration at the end point was not specified a priori in the computational model, and thus, the predicted and measured final joint angles may differ. In particular, several examples show large deviations for the humeral angle ζ in the frontal plane movement conditions. Calculations of the forward kinematics from the predicted joint angles confirm that the large differences in the predicted versus measured humeral rotation are compensated by much smaller differences in other joint angles, such that the discrepancy is indeed restricted to the null space and the boundary condition of reproducing the final hand position has been met. One explanation for this discrepancy might be that this degree of freedom is least controlled by the system. Note that in the case of a fully extended arm, the humeral angle ζ is identical to the rotation angle α (up to a constant) around the axis going through the shoulder location and the final hand location. Recall that the angle α parameterizes the null space of the endeffector location, and thus, changes in the angle α do not contribute to the task goal.
Dynamic level: driving torques
The computational model derived at the temporal and geometric levels leads to interesting properties at the dynamic level. These properties will be derived next.
First, we analyze the driving torques needed to generate the movement along the geodesic paths. The dynamic equations of motion for the arm are given by where N denotes the vector of gravitational torques. Frictional torques are neglected. The joint torques, τ, generated by muscle forces are divided into the driving torques, τ_{d}, and the torques, τ_{g}, that counterbalance the external torques generated by the gravitational field. Thus, we assume that gravitation does not contribute to the driving torques along the path in configuration space, but leads to a static posture maintenance at each location in task space. Therefore, we decompose the torques as follows: With these assumptions, the dynamic equations (Eq. 48) transform into As pointed out in previous research works (Flash and Hollerbach, 1982; Atkeson and Hollerbach, 1985), the separation of torques into posture maintenance torques and gravityindependent driving torques might be used by the motor system to simplify arm dynamics. The advantage of such a strategy derives from the scaling properties of the driving torques under a change of speed. Consider a trajectory that results from a given trajectory by time scaling t̃ = kt with some constant k.
Then, hand speed scales according to For k < 1, the movement speed is increased, whereas for k > 1, the movement is slowed down. The driving torques scale according to Equation 51 as For example, to move the hand twice as fast along a given hand path, the driving torques have to be multiplied by a factor of four.
We further analyze the driving torques by assuming that the path in configuration space is a priori specified and the jointangular trajectory can be represented in the form q(t) = q(λ(t)), where λ = σ/Σ is the normalized arc length in configuration space. Using the chain rule, the equations of motion in Equation 51 can be written as where a prime and a dot denote differentiation with respect to λ and time, respectively.
For movements along geodesic paths, the term in brackets on the left side vanishes according to Equation 31. The vanishing expression consists of inertial torques that depend on the speed in configuration space, τ_{2} = M(q)q″λ̇^{2}, and the torques that are commonly denoted as the centrifugal and Coriolis torques (Flash and Hollerbach, 1982), τ_{3} = C(q, q′)q′λ̇^{2} = C(q, q̇)q̇; i.e., along geodesics, it is τ_{2} + τ_{3} = 0. The remaining term τ_{1} = M(q)q′λ̇ describes inertial torques that depend linearly on the acceleration in configuration space. Thus, the arm dynamics for movements along geodesic paths with nonconstant speed is governed by Note that the movement along geodesic paths results in a much simpler dynamics than given in Equation 48. We remark further that the Coriolis and centrifugal interaction torques, τ_{3}, can still be present for arm movements along geodesic paths as long as they compensate the speeddependent inertial torques τ_{2}.
To analyze further the dynamical properties of movements along a geodesic path, we compare Equation 54 to Newton's equation of motion for a mass point m along a prespecified path, x(s), where s is the Euclidean arc length along the path. Newton's equation of motion, F = mẍ, transforms with ẋ = sẋ′ and ẍ = s̈x′ + ṡ^{2}x″ to where a dot and a prime denotes differentiation with respect to time and arc length s, respectively. The vectors t(s) = x′(s) and n(s) = x″(s)/κ(s) are the tangential and normal vector to the path, respectively, and κ(s) denotes the curvature of the path. Note the similar structure of Equations 54 and 56. If the mass point moves on a straight path (geodesics in Euclidean space) where x″ = 0, accelerations act according to Equation 56 only in directions tangential to the path, whereas accelerations normal to the path disappear.
Similar results hold for movements of the arm according to Equation 54. For further analysis, it is useful to multiply Equation 54 with the inverse of the metric, leading to For arm movements along geodesic paths (“straight” paths in Riemannian configuration space), the expression in the square brackets of Equation 57 disappears, and thus all accelerations act in direction tangential to the path (i.e., along q′), whereas all nontangential accelerations vanish. Movements along geodesics thus require less muscular effort, because only forces that induce acceleration in direction tangential to the path in configuration space have to be provided by the muscles. Any deviation from the straight (geodesic) path requires additional muscle forces. These forces are needed to keep the arm on the nongeodesic path (“curved” path in Riemannian configuration space) comparable with the centripetal force needed to keep a mass point on a (curved) circular path in Euclidean space.
Several examples of predicted and measured driving torques in the radial and frontoparallel movement conditions are shown in rows in Figures 16 and 17, respectively. The measured driving torques were obtained by evaluating the left side of Equation 51 using the experimental data. The predicted driving torques were derived from Equation 55 using the prediction of the GEO model. The first column in Figures 16 and 17 shows the total driving torque. The middle column shows the inertial torque τ_{1}, and the last column depicts the sum of all configurationspeeddependent torques, τ_{2} + τ_{3}. Note that according to the GEO model, τ_{2} + τ_{3} = 0, and thus the inertial torque τ_{1} is equal to the total driving torque τ_{d}. The measured torque profiles resulted in larger variability around the predicted curves. It is worth noting that the sign and size of the torque amplitudes were well accounted for by the model but the measured torque profiles showed larger fluctuations. We can think of several reasons for this discrepancy. First, the path does not always follow a geodesic path, and thus, Equation 55 does not hold exactly. Second, observed speed profiles show deviations from a symmetric profile, and thus, the assumption of MJ speed profile along the hand path is too stringent and cannot account for the small fluctuations in the torque profiles. Finally, the extraction of torques from position measurements is not trivial, and noise may be induced during data processing. Also note that the sum of measured torques, τ_{2} + τ_{3}, nearly vanishes in the presented examples, which suggests that most of the driving torques can be attributed to the torques τ_{1}. These results are in good agreement with the GEO model.
The implication of Equation 55 in task space can be derived by introducing the operational force F at the end effector and by reparameterization with respect to arc length, s, of the Euclidean hand path. The driving torques τ_{d} are related to the operational force F by the following fundamental relation (Kathib, 1987): Inserting Equation 58 into Equation 57 and multiplying both sides with the hand Jacobian leads for geodesic paths to where we have used the relation x′ = J_{h}q′, the definition of the mobility matrix W = J_{h}M^{−1}J_{h}^{T}, and the fact that the expression in square brackets in Equation 57 disappears along geodesic paths.
Next, we reparameterize Equation 59 with respect to Euclidean arc length s by inserting the function λ = λ(s) and Equation 44 into Equation 59. With dx/dλ = (dx/ds)(ds/dλ) = (dx/ds)(1/[λ′(s)]), we get where t = dx/ds is the unit tangent vector to the hand path and a prime now denotes differentiation with respect to s.
The second term in the square brackets of Equation 60 depends on the curvature of the hand path. For quasistraight hand paths, however, the term can be neglected, because for such paths, λ″(s) ≈ 0; i.e., where a_{t} is the tangential acceleration.
Note also that for quasistraight hand paths, the driving torques along geodesic paths in configuration space are according to Equation 55 well approximated by resulting in a linear relation between the hand acceleration s̈ and the driving torques. The predictions of the GEO model are compatible with the assumption of quasistraight hand paths as shown in Figures 4 and 5 as well as the experimental measured hand paths that lead according to Figures 9 and 10 to quasistraight handpath distributions of the curvature index.
We conclude from Equation 61 that the tangential acceleration in task space is identical to the product of the mobility matrix, which measures the ease of accelerating the hand in a certain direction, and the actual force. Conversely, for movements along quasistraight hand paths that were derived from geodesic paths in configuration space, the product of mobility matrix and operational force (effective acceleration) is always pointing in direction tangential to the hand path. The total acceleration of the end effector follows then as a = a_{t} + a_{n}, where a_{n} is the normal acceleration given by a_{n} = κṡ^{2}n. The vector n denotes the normal vector and κ the curvature of the hand path. It should be noted that the tangential and normal vectors as well as the curvature of the hand path are known quantities that follow from the projection of the geodesic path into task space.
We will next analyze what the selected temporal feature (speed) of the computational model implies at the dynamical level. We first consider the minimum torquechange model (Uno et al., 1989; Nakano et al., 1999; Wada et al., 2001). We analyze what this criterion would imply assuming that the movement is along a geodesic path and the speed along the hand path, ṡ(t), is the only free variable to be determined. Along a geodesic path, the cost associated with the minimum torquechange model is where we have used Equation 62 and defined A(t) ≡ M(q(s(t)))q′(s(t)). We noticed that for geodesic paths, Ȧ(t) ≈ 0, and thus the change of driving torque is to a first approximation proportional to the jerk of the arc length: We conclude that the minimization of the squared driving torque change integrated over the movement is roughly equal to a minimization of the squared jerk of the arc length of the hand path integrated over the movement time. The minimumjerk model and the minimum torquechange model are thus compatible when restricted to geodesic paths in configuration space. Moreover, as shown in Appendix D, the chosen motor patterns at the geometric and temporal levels are equivalent to a minimization of the peak value of the kinetic energy. Thus, the minimumpeak kinetic energy model as suggested by Soechting et al. (1995) is an additional outcome of our model.
A second, alternative solution for the speed determination along the geodesic path may consist of a cost in task space in the form of the squared change in effective acceleration integrated over the movement time, thus: The cost can be rewritten according to Equation 61 as where we used the relationships t′(s) = κ(s)n(s) and t(s) · n(s) = 0. Because for quasistraight hand paths, κ (s) ≪ 1, the second term in the last expression of Equation 66 can be neglected, and the cost transforms into the modified minimumjerk model that defined the temporal properties of the movement in the GEO model.
Discussion
In this study, we derived a computational model for 3D pointtopoint movements that reconciles different existing models into one computational framework. The observed invariance of path and posture with changes in speed led to the assumption that geometric and temporal aspects of movement are decoupled. This resulted in a model that predicted geometric and temporal properties independently. Geometric properties were determined in configuration space that, when equipped with a suitable metric, defined a Riemannian manifold. The chosen metric was the kinetic energy metric, which relates closely to the arm dynamics. The optimal jointangular paths were assumed to follow geodesic paths in this metric space. The timing of the movement was selected in task space by assigning a minimumjerk speed profile along the hand path.
At the dynamic level, we have shown that geodesic paths in configuration space correspond to arm movements with less muscular effort, because only forces that induce acceleration in direction tangential to the geodesic path have to be provided by the muscles. For any other (nongeodesic) path, additional muscle forces have to be provided to generate the configurationspeeddependent interaction torques. Movements along geodesic paths with quasistraight hand paths require driving torques that are linearly related to the tangential acceleration. It follows that the arm dynamic equations along geodesic paths are significantly simplified.
Moreover, we have shown that geodesic paths together with a minimumjerk speed profile imposed along the extrinsic hand path corresponded to jointangular trajectories with minimal peak kinetic energy over all final accessible arm postures at the target location. Our model is thus compatible with the minimum peakwork model suggested by Soechting et al. (1995), and predicts not only the final arm posture as presented by that model, but also the entire hand path and the jointangular trajectories.
Our study also provides a dynamic interpretation of the minimumjerk model by showing that a maximal smoothness criterion imposed on the handspeed profile leads to a near minimization of the squared change of driving torques along the selected geodesic path. The observed doublepeaked torque profiles may be the result of the simplified dynamics along geodesics. Two prominent and apparently contradictory models, the minimumjerk model (Flash and Hogan, 1985) and the minimum torquechange model (Uno et al., 1989; Nakano et al., 1999; Wada et al., 2001), could thus be reconciled within the framework of geodesics.
We do not propose that the CNS has an internal representation of Riemannian geometry or geodesics. We hypothesized that geodesics are an emerging property of the system that seeks to reduce muscular efforts possibly based on proprioceptive feedback. Along geodesic paths, the sum of configurationspeeddependent inertial torques and centrifugal and Coriolis torques disappears, and only inertial torques that depend on the acceleration in configuration space remain. It is important to note that centrifugal or Coriolis interaction torques along geodesic paths are still possible provided that the net torques, i.e., the sum of the torques τ_{2} and τ_{3}, disappear. These findings are interesting in light of evidence suggesting that during multijoint movements, muscle activation and neural control strategies have evolved to compensate for interaction torques (Flash and Hollerbach, 1982; Sainburg et al., 1993; Gribble and Ostry, 1999). An equilibrium trajectory control strategy (Feldman, 1966, 1986; Flash, 1987) may enable movement by continuously shifting the arm's equilibrium point along the predefined geodesic path.
Validity of the proposed model
The predictions of the model for jointangular paths, final posture, hand paths, and speed profiles were in good agreement with our experimental data.
Among the three tested models, the geodesic model accounted best for the observed hand paths and final arm postures. This seems to support previous results demonstrating the role of inertial properties in the selection of the final arm posture (Soechting et al., 1995; Flanders et al., 2003) and the independence of final posture from movement speed and gravity (Nishikawa et al., 1999). In contrast, the MTC model resulted in strongly curved hand paths and unrealistic arm postures. The overall fit of the minimumjerk speed profile was very good, although observed speed profiles frequently showed slightly larger amplitudes and deviations from the symmetric shape. The speed profiles computed from the MTC model resulted in toosmall amplitudes that showed in several cases unrealistic double peaks. The sign and size of torque amplitudes were well accounted for by the model, but all the fine details of the drivingtorque profiles could not be explained with the present model.
Simulation studies showed that the application of the geodesic model to planar movements resulted in strongly curved hand paths. This poses the question of whether there is a fundamental difference in the planning and control of movements in 2D versus 3D space. One possible explanation for this discrepancy may lie in the fact that planning of movements constrained to a horizontal plane is more strongly influenced by visual feedback and by the wish to move along straight lines as opposed to unconstrained 3D pointtopoint movements and movements in the vertical plane (Atkeson and Hollerbach, 1985; Desmurget et al., 1997).
Evidence for an independent control of temporal and spatial aspects of movement in CNS
Neurophysiological evidence suggests that spatial aspects of hand paths and their evolution over time may be planned at different independent levels. For example, movement direction toward visual targets (Georgopoulos et al., 1988) and changes in movement direction (Georgopoulos et al., 1989) are encoded by directionally tuned cells within the primary motor cortex (Kettner et al., 1988), suggesting that the motor cortex may be one of the specific sites for processing spatial information before the actual movement execution (Georgopoulos, 1995; Scott, 2000).
More complex spatial information related to movement (e.g., during drawing geometrical patterns) seems to be mediated via neural activity in prefrontal and premotor regions (Moran and Schwartz, 1999; Averbeck et al., 2002). In a recent study (Torres and Andersen, 2006), the space–time separation has been shown during the learning of an obstacleavoidance task in monkeys. It has been proposed that the geometric planning stage is performed by the posterior parietal cortex.
On the other hand, timing is a fundamental feature of motor skill that seems to be widely represented in the CNS (Rao et al., 1997). Previous neurophysiological and brain mapping studies have described a highly distributed network of cortical and subcortical areas that code temporal aspects of the movements (Harrington et al., 1998; Bengtsson et al., 2005; Buhusi and Meck, 2005; Xu et al., 2006). These subcortical areas may include both the cerebellum and the basal ganglia. For example, for simple movements, temporal information is encoded in the anterior and posterior regions of the cerebellar cortex via climbing fibers [e.g., eyelid conditioned responses and compensatory eye movements (Raymond et al., 1996); discontinuous movements of the hand (Spencer et al., 2003); finger tapping movements (Ivry and Spencer, 2004)].
Last, integration of different features of movement, such as space and time, may be confined to premotor and frontal brain regions, which interact with the cerebellum, cortical (e.g., motor cortex), and subcortical (basal ganglia) motorrelated areas (Sakai et al., 2000, 2002; Garraux et al., 2005).
Effect of gravity
Our model assumed that gravity has no effect on the driving torques. In previous studies, it has been reported that hand paths of movements in the sagittal plane tend to be more curved in the upward than in the downward direction (Atkeson and Hollerbach, 1985; Papaxanthis et al., 2003). It has been concluded that this effect can be assigned to the gravitational field that acts, respectively, against and along the movement direction, respectively.
Within the framework of our model, however, the observed asymmetry in curvature for upward versus downward movements is caused by the different initial arm configurations in upward and downward movements and not by the presence of the gravitational field.
To conclude, we propose that the geometrical movement properties are acquired in nonEuclidean configuration space, independently from the temporal ones, based possibly on proprioceptive feedback, which leads the arm to reduce muscular effort. Once such properties are internalized, they may form part of a trajectory formation strategy that is implemented using feedforward control. Such a strategy integrates the two independent levels into a third spatiotemporal level. We showed that if the controller adopts such a solution, the arm joints move with near minimal torque change and minimal peak work while bringing the hand to the required location. Our data validated this assumption. In the light of these results, neurophysiologists may be encouraged to determine in more detail the specific brain regions involved in the space–time separation process during movement. Finally, our study showed the use of Riemannian metric spaces in the description of human arm movements, and it resulted ultimately in a computational model that reconciled seemingly opposing models contrasted in the literature.
Appendix A
The components of the elbow Jacobian for the four DOF linkage arm model defined by the joint angles q:= (θ, η, ζ, φ)^{T} ∈ Q can be derived from the forward kinematic relations x_{e} = F_{e}(q). The components of the elbow Jacobian are then defined as J_{e},_{ij} = ∂F_{e},_{i}/∂q^{j}, (i = 1, 2, 3; j = 1, …, 4), leading to the nonzero components Similarly, for the hand Jacobian J_{h},_{ij} = ∂F_{h},_{i}/∂q^{j}, we obtain the nonzero components where l_{u} and l_{f} denote the length of the upper and forearm, respectively.
Appendix B
For a flexed arm and a fixed hand position x_{h}, the arm can still rotate around an axis going through the shoulder and the fixed hand position (Fig. 1B). The allowed elbow positions are contained in a circle C around this axis (Hollerbach, 1985). The circle can be obtained as the intersection between the sphere S_{1}: x^{2} = l_{u}^{2} around the origin with radius l_{u} and the sphere S_{2}: (x − x_{h})^{2} = l_{f}^{2} around the hand location x_{h} with radius l_{f}, giving the elbow location as a function of the rotation angle α ∈ [0, 2π) of the intersection circle C = S_{1} ∩ S_{2}. After algebraic manipulation, the elbow location can be expressed as with the radial distance to the center of the intersection circle u=(1/2d)(l_{u}^{2}−l_{f}^{2}+d^{2}), intersection circle radius and d = x_{h}. Furthermore, it is ϕ = atan2(y_{h}, x_{h}), ϑ = asin(z_{h}/d), e_{x} = (1, 0, 0)^{T}, e_{r}(α) = (0, cos α, sin α)^{T}, where atan2(a, b):= atan (a/b) − sign(a)[1 − sign(b)](π/2). The angles α, ϕ, and ϑ are defined in Figure 1B. R_{y}, R_{z} define rotation matrices around the y and zaxes, respectively, and are given by
Appendix C
The kinetic energy of a four DOF arm in terms of the coordinates q = (θ, η, ζ, φ)^{T} ∈ Q is determined by with and the constants I_{1} = I_{u},_{x} + m_{u}a_{u}^{2} + m_{f}l_{u}^{2}, I_{2} = I_{u},_{z}, I_{3} = I_{f},_{x} + m_{f}a_{f}^{2}, I_{4} = I_{f},_{z}, A = m_{f}l_{u}a_{f} (Biess et al., 2001). The parameters m_{i}, I_{i},_{x}, I_{i},_{y}, I_{i},_{z}, l_{i}, r_{i}, (i = u, f) denote mass, principal moments of inertia around the x, y, and zaxes of the bodyfixed coordinate systems, length, and distance to the center of mass of the upper and forearm, respectively. It was assumed that the (transversal) x and y components of the moment of inertia for the upper arm and forearm, respectively, are the same; i.e., I_{u},_{x} = I_{u},_{y} and I_{f},_{x} = I_{f},_{y}. The kinetic energy is a quadratic form in the joint velocities and can be written as M ∈ ℝ^{4×4} is the kinetic energy metric (manipulator inertia matrix) with nonzero components and M_{ij} = M_{ji}.
Appendix D
The selection of the motion pattern on the temporal and geometrical as suggested in this work leads to the minimization of the peak value of kinetic energy and is therefore in accordance with the model suggested by Soechting et al. (1995). This result can be derived from our computational model as follows. For each path in configuration space, γ_{α}, α ∈ I(x_{f}), that connects a given initial arm configuration with the set of accessible final arm configurations compatible with a given hand location x_{f}, the kinetic energy is given by the following (derived from Eq. 45): where Σ is the total arc length of the path in configuration space, x′(λ) denotes the derivative of the hand path, and s is the Euclidean arc length of the hand path.
We determine next the single path out of the oneparameter family of geodesics that has minimal peak kinetic energy. The necessary condition for the peak value of the kinetic energy is where t = t* denotes the peak time. If we assume that a minimumjerk profile is imposed on the temporal level, Equation 98 leads to resulting in τ* = 1/2, where τ is normalized time τ = t/T. The peak of the kinetic energy thus occurs in the middle of the movement, which is in reasonable agreement with the experimental data (Figs. 9, 10). The peak value of the kinetic energy, K̂, as a function of the parameter α, follows then as where S(α) is the total length of the hand path corresponding to path γ_{α}. For hand paths that are quasistraight, as obtained in the computational model, it is S(α)/x′(λ, α) = ∫_{0}^{1}x′(λ, α)dλ/x′(λ, α) ≈ const, and the peak value of kinetic energy is proportional to the squared arc length of the path in configuration space. The geodesic with the shortest arc length that connects the initial and the final arm configurations results thus in the minimization of the peak value of kinetic energy.
Footnotes

This work was supported in part by the Minerva Foundation, Germany, and in part by the DIP (German–Israeli Project Cooperation). A.B. acknowledges support from a Dov Biengun research fellowship and thanks Dr. J. Reingruber for stimulating discussions. T.F. is the incumbent of the Dr. Hymie Moross Professorial Chair.
 Correspondence should be addressed to Armin Biess, Department of Mathematics, Weizmann Institute of Science, 76100 Rehovot, Israel. armin.biess{at}weizmann.ac.il