BM 1360

BRR BINNY SEETHARAMA

Journal of Biomechanics 000 (2000) 000}000

A method for measuring endpoint sti!ness during multi-joint arm
movements
E. Burdet    *, R. Osu, D. Franklin, T. Yoshioka, T.E. Milner, M. Kawato 
Department of Mechanical Engineering, National University of Singapore, 119260 Singapore
School of Kinesiology, Simon Fraser University, BC, Canada
Kawato Dynamic Brain Project, ERATO, JST, Kyoto, Japan
ATR Human Information Processing Research Laboratories, Kyoto, Japan
Received 24 March 2000; accepted 5 June 2000

F

Abstract

PR

O

O

Current methods for measuring sti!ness during human arm movements are either limited to one-joint motions, or lead tosystematic errors. The technique presented here enables a simple, accurate and unbiased measurement of endpoint sti!ness duringmulti-joint movements. Using a computer-controlled mechanical interface, the hand is displaced relative to a prediction of theundisturbed trajectory. Sti!ness is then computed as the ratio of restoring force to displacement amplitude. Because of the accuracy ofthe prediction (( 1 cm error after 200 m s) and the quality of the implementation, the movement is not disrupted by the perturbation.
This technique requires only  as many trials to identify sti!ness as the method of Gomi and Kawato (1997, Biological Cybernetics 76,
163}171) and may, therefore, be used to investigate the evolution of sti!ness during motor adaptation.  2000 Elsevier Science Ltd.
All rights reserved.

TE

D

Keywords: Human arm; Voluntary movements; Sti!ness measurement; Robotic interface

1. Introduction

U

N

C

O
R

R

EC

Endpoint sti!ness of the arm, resulting from the resistance of muscles to a small displacement of the hand, is an
important mechanical property of the musculoskeletal
system, characterizing resistance to disturbances encountered in the environment (Hogan, 1985). When the handis perturbed slightly, it tends to return to its original
position. Sti!ness can be computed as the ratio of restoring force to displacement amplitude.
Using this technique, the features of endpoint sti!nessduring posture have been investigated extensively(Mussa-Ivaldi et al., 1985; Tsuji et al., 1995; Gomi andOsu, 1998). However, measuring sti!ness during multijoint movements remains a technical challenge. A rigid,
powerful mechanical interface with computer-controlled
dynamics is also needed to control hand position during
movement (Gomi and Kawato, 1997).
* Correspondence address: Department of Mechanical and Production Engineering, National University of Singapore, 119260, Singapore.
Tel.: 65-779-1459;.
E-mail address: e.burdet@ieee.org (E. Burdet).

Most previous estimates of sti!ness during movement
have been realized by means of force perturbations.
Bennett and others (Bennett et al., 1992; Lacquaniti et al.,1993) have used stochastic force disturbances and measured the resulting change in hand position. This method
is relatively straightforward to implement, but results in
sti!ening of the joints due to muscle co-contraction (Milner, 1993) and can reduce or even abolish the stretch
re#ex (Stein and Kearney, 1995). For these reasons, Bennett (1993) later employed a displacement relative to themean undisturbed movement for determining joint sti!ness during elbow movements. Although the mechanicalinterface for displacing the hand must be precisely controlled, only a few trials are needed and data analysis isstraightforward (divide force by displacement amplitude).
Force impulses have been used to estimate sti!nessduring multi-joint movements (Gomi and Kawato, 1997).
However, the method su!ers from the limitation thata perturbation of the same amplitude, applied at di!erent
points in the trajectory or in di!erent directions, will
displace the hand by di!erent amounts. This is because
limb sti!ness depends on joint angles, angular velocity
and perturbation direction (Mussa-Ivaldi et al., 1985;

0021-9290/00/$ - see front matter  2000 Elsevier Science Ltd. All rights reserved.
PII: S 0 0 2 1 - 9 2 9 0 ( 0 0 ) 0 0 1 4 2 - 1

BM 1360 BRR BINNY SEETHARAMA
E. Burdet et al. / Journal of Biomechanics 000 (2000) 000}000

F

Fig. 1. Elastic force resulting from a position shift in the y-direction(movement of Fig. 3). The "gure shows the perturbed trajectory minus
the undisturbed one. The perturbation was simulated on free movements measured by the interface. The elastic force was computed during
the constant phase of the perturbation (note the di!erent time scale ofthe force plots). In the position and velocity plots, the dotted lines show
the error in 10 trials, the solid line their mean, and the dashed lines themean plus/minus standard deviation. In the force plots, the dotted lines
show the estimation resulting from Eq. (4), and the solid lines theirmean. The large dots correspond to the force that would result froma perfect trajectory prediction. In the Cartesian force plot, the # correspond to the exact value, and the `,a to the estimated value.

where

dF
dF du
!
K(x,x ,u),!
du dx
dx

EC

TE

Let q be the two-dimensional vector with the shoulder
and elbow as "rst and second coordinates, and x thecorresponding two-dimensional Cartesian hand position
vector. Let F (q(t),q (t),qK (t)) be the (endpoint) force required to move the arm along the trajectoryq(t),0)t)¹, and let F (q(t),q (t),u(t)) be the (endpoint)
force developed by the muscles activated by u to realize
this movement. The components of u correspond to all
the muscles involved in the movement. If there is nointeraction with the environmentD

2. Methods

F (q,q ,qK ),F (q,q ,u)

(1)

O
R

R

i.e., the sum of the internal forces acting on the musculoskeletal system is 0. If a force F is exerted on the hand (bythe mechanical interface), the internal forces are equal tothis external force
F (q,q ,qK )!F (q,q ,u),F.

(2)

Ne "eK ,0.

C

Let e represent a position perturbation in hand Cartesian
coordinates with a plateau phase (Fig. 1), during which
the following condition is satis"ed:(3)U

By linearizing Eq. (2) about the unperturbed trajectory,
we "nd that the change in force due to the perturbatione corresponds to
dF
e#Ke,
*F"
dx

(4)

PR

O

Gomi and Kawato, 1997; Bennett, 1993). Since the sti!ness depends on the displacement amplitude (Shadmehr
et al., 1993), bias may result. Furthermore, damping andsti!ness must be identi"ed together, so many trials areneeded with perturbations in multiple directions.
Consequently, we decided to adapt Bennett's method
to study the sti!ness of two-joint arm movements. Thisrequired consideration of nonlinear limb dynamics, andaccurate two-dimensional prediction of where the trajectory would have gone, had there been no perturbation, inorder to achieve a constant displacement. For one-joint
movements, the mean of past trajectories is an acceptableestimate of the unperturbed trajectory, whereas this is
generally not true for multi-joint movements because
they are less constrained.
This paper presents the key features of our techniquefor estimating sti!ness during movement: an accurate
prediction algorithm for multidimensional trajectories
and the implementation of a servo-controlled perturbation during movement. The method is tested in simulations and its e$ciency for measuring sti!ness without
disrupting motion is demonstrated in experiments conducted with two subjects.

O2

is the endpoint sti!ness. Note that K depends on re#exes.
*F corresponds to the mean force measured during perturbed movements minus the mean force required toperform unperturbed movements. We found in simulations that (dF /dx)e can contribute up to 10% of *F so itcannot be neglected.
Since the trajectory prediction is accurate for only
a short time interval (Fig. 1), a good position perturbation must be brief, and the constant position plateau
must be reached as quickly as possible. However, a very
fast or abrupt transition to the plateau requires themechanical interface to produce high forces within a brieftime, which can lead to vibration. A short transition
phase, which miminized vibration (Fig. 1), was achievedby using as servo command a sixth-order polynomialwith zero velocity and zero acceleration at the boundaries and zero end jerk (Burdet and Osu, 1999).
The method requires that there be an interval during
which the perturbation moves the hand at the velocity of
the unperturbed trajectory. To predict the unperturbed
trajectory, we noted that trajectories of repeatedBM 1360

BRR BINNY SEETHARAMA

E. Burdet et al. / Journal of Biomechanics 000 (2000) 000}000

velocity pro"les were generated (second panel of Fig. 2).
The candidate best matching the current trajectory wasthen selected, and used for predicting the unperturbed
trajectory (third panel of Fig. 2). The best match was
determined by minimizing the following recursive function once the velocity had crossed a threshold of 0.03 m s:
d "(v !w )#ad
,
(5)
GI
I
GI
G I\
where v is the velocity of the current time sample k, w
IGI
is the ith velocity candidate at time k and a"0.94
is a forgetting factor which limits matching to the last
100 m s. The movement was decomposed parallel andperpendicular to the main movement direction and candidate velocity pro"les were selected independently in thetwo directions.
The template velocity pro"le, computed as the mean of
10 past movements, was updated after each trial. The
standard deviation was computed with the most recent
trial replacing the earliest of the 10 stored trials. If thestandard deviation was larger than 1.05 times the previous
standard deviation, then the most recent trial was excluded, otherwise it was included (Burdet and Osu, 1999).
A horizontal planar parallel manipulator, powered by
two DC direct-drive motors controlled at 2 kHz, wasused to apply the perturbation (Gomi and Kawato,1997). The motor shaft angles were measured using optical encoders (409600 pulse/rev) while a force-torquesensor (Nitta Corp. No. 328) mounted on the handle
measured the force exerted on the hand to a resolution of
0.06 N.
The speed and acceleration were obtained from theposition signal using, respectively, second- and thirdorder Butterworth "lters with a cut-o! frequency of 50
Hz, such that di!erentiation and "ltering were realized
simultaneously. A feedforward term compensated for theinertia (i.e., the mass matrix multiplied by the acceleration) and the friction, but not for the velocity dependent
forces. A relatively small force ((5 N) was required toperform arm movements (Fig. 3).
The current trajectory was predicted at 1 kHz. The
position perturbation was realized using Cartesian
PD-trajectory control. The commanded trajectory corresponded to the predicted trajectory plus the perturbation. Note that there was no trajectory control before orafter the perturbation. The perturbation had an amplitude of 8 mm and was 300 m s in duration (Fig. 1).
Transition phases of up to 100 m s were necessary tominimize vibration. The e!ect of vibration was avoided
by restricting the sti!ness identi"cation to a 60 m s interval centered at the midpoint of the plateau phase. Notethat the measured sti!ness corresponds to a temporalmean over these 60 m s.
The prediction algorithm was tested by simulating the
perturbation on free movements measured with theinterface. The arm was modeled as a double pendulumwith lengths (0.3,0.3) m and masses (1.93,1.52) kg for the

C

O
R

R

EC

TE

D

PR

O

O

F

movements performed under the same conditions weresimilar ("rst panel of Fig. 2). As the velocity pro"le of any
movement was roughly similar to the mean velocity
pro"le of the preceding trials, the shape of the meanpro"le was used as a template for the current trajectory.
The velocity pro"les were low-pass "ltered at 50 Hz,
scaled to produce the same amplitude (in the main movement direction) and truncated at both ends using a 0.03
m s velocity threshold (in the main movement direction)prior to averaging. To avoid distortion of the template
due to di!erences in the times of the peak velocity, themean was computed independently for the accelerationand deceleration phases. The template velocity wasscaled over a range of 11 amplitudes and 11 time shifts,
chosen to minimize the prediction error over 20 unperturbed movements. In this way, 11;11"121 candidate

3

U

N

Fig. 2. Steps of the prediction algorithm. The "rst panel shows thevelocity of 10 trials in the main movement direction (measured data).
The mean of the 10 velocity measurements is computed, scaled inamplitude and shifted in time to generate 121 velocity pro"le candidates. Some of these candidates are plotted in the second panel. At any
time, the closest candidate can be used to predict the future of thecurrent trajectory. This is illustrated in the third panel. In this graph,the solid lines correspond to the current velocity, and the dashed lines
to a prediction carried out every 0.1 s.

BM 1360 BRR BINNY SEETHARAMA
E. Burdet et al. / Journal of Biomechanics 000 (2000) 000}000

D

Fig. 4. The same movement than in Fig. 3, with a perturbation of
amplitude 8 mm in the y direction. The solid line corresponds to themean free movement, and the dotted lines to 10 perturbed movements.

sometimes positive, sometimes negative, and distributed
symmetrically around the unperturbed velocity (Fig. 1).
The mean prediction error is close to 0, and the mean
force is close to the force corresponding to perfect trajectory prediction. This indicates that the prediction algorithm and the sti!ness estimation are unbiased andaccurate.
Tests with a perturbation of amplitude 0 showed thatthe mechanical interface did not disrupt voluntary movement while driving the hand (Fig. 3). The subjects wereunable to detect the perturbation nor did it appear to
a!ect the trajectory importantly. The distance between
the means of 10 free movements and of 10 constrained
movements was smaller than the standard deviation of
the free movements for the position, as well as for theforce. Similar consistency was found with the 8 mm
amplitude perturbation (Fig. 4). The force was #at, i.e.
there was no oscillation, in the 60 m s during which
sti!ness was identi"ed.
The sti!ness values obtained for the two subjects at the
midpoint of the movement, listed in Table 1, are of
the same order as those of Gomi and Kawato (1997). The
deviation is small enough to test hypotheses about thein#uence of the dynamic environment on sti!ness

TE

Fig. 3. Position and force measured during movements of amplitude
0.4 m and duration 1 s to the right (Cartesian coordinates relative tothe shoulder). The solid and dashed lines correspond, respectively, tothe mean of 10 free movements plus/minus standard deviation. The
dotted lines are the mean and the mean plus/minus standard deviationof 10 movements with the hand constrained along the prediction of the
movement between 0.3 and 0.6 s (relative to a 0.3 m s velocity threshold).

PR

O

O

F

4

U

N

C

O
R

R

EC

upper arm and the lower arm, respectively. The corresponding dynamics F and corrective term (dF /dx)ewere computed numerically. The sti!ness matrix K,[30 10; 10 20] N m was taken from (Gomi and Kawato,1997).
The implementation was tested on horizontal movements performed by two subjects. The experiments wereapproved by the institutional ethics committee, and thesubjects gave informed consent prior to participation.
A perturbation which constrained the hand for the predicted trajectory (0 amplitude) was "rst tested (Fig. 3),
and then a perturbation of amplitude 8 mm (Fig. 4).
Finally, sti!ness was identi"ed using 40 movements with
perturbation directions randomly chosen from the set+k(2p/8), k"128, (Burdet et al., 1999).

3. Results

The simulation showed that the prediction error issmall ((1 cm after 200 m s). The prediction error isBM 1360

BRR BINNY SEETHARAMA

E. Burdet et al. / Journal of Biomechanics 000 (2000) 000}000
Table 1
Sti!ness measured in two human subjects for the movement of Fig.
3 compared with the approximate value of a typical subject in Gomi
and Kawato (1997)(N m/rad)Subject ETI

Subject OSU

Gomi97

K
Q
K
QC
K
CQ
K
C

40.3$3.0
21.0$2.3
12.9$3.0
21.3$2.3

23.6$3.1
7.4$2.2
9.5$3.1
14.7$2.2

30
10
10
20

gradual changes of the trajectory, it may enable us to
examine the evolution of sti!ness during motor adaptation.

Acknowledgements

The sti!ness, measured in Cartesian space, has been transformed into
joint coordinates K ,[K K ; K K ] as in McIntyre et al. (1996).
O
Q QC CQ C
The $ de"ne 95% con"dence intervals.

The experiments were performed at ATR. We thank H.
Gomi and T. Yoshioka for having set up the testbed, andthe reviewers for their valuable remarks. This research
was supported by ERATO/JST, Japan, JISTEC, Japan,
and by grants from the Swiss National Science Foundation and the Natural Sciences and Engineering ResearchCouncil of Canada.

References
Bennett, D.J., 1993. Torques generated at the human elbow joint inresponse to constant position errors imposed during voluntary
movements. Experimental Brain Research 95, 488}498.
Bennett, D.J., Hollerbach, J.M., Xu, Y., Hunter, I.W., 1992. Timevarying sti!ness of human elbow joint during cyclic voluntary
movement. Experimental Brain Research 88, 433}442.
Burdet, E., Osu, R., 1999. Development of a new method for identifying
muscle sti!ness during human arm movements. Report 1-21,
Kawato Dynamic Brain Project, ERATO, JST, Japan.
Burdet, E., Osu, R., Franklin, D., Milner, T.E., Kawato, M., 1999.
Measuring sti!ness during arm movements in various dynamic
environments. 1999 ASME Annual Symposium on Haptic Interfaces and Virtual Environments for Teleoperator Systems, Nashville,
USA, 421}428.
Gomi, H., Kawato, M., 1997. Human arm sti!ness and equilibriumpoint trajectory during multijoint movement. Biological Cybernetics 76, 163}171.
Gomi, H., Osu, R., 1998. Task-dependent viscoelasticity of humanmultijoint arm and its spatial characteristics for interaction withenvironments. Journal of Neuroscience 18, 8965}8978.
Hogan, N., 1985. The mechanics of multi-joint posture and movement
control. Biological Cybernetics 52, 315}331.
Lacquaniti, F., Carrozzo, M., Borghese, N.A., 1993. Time-varying
mechanical behavior of multijointed arm in man. Journal of Neurophysiology 69, 1443}1464.
McIntyre, J., Mussa-Ivaldi, F.A., Bizzi, E., 1996. The control of stablepostures in the multijoint arm. Experimental Brain Research 110,248}264.
Milner, T.E., 1993. Dependence of elbow viscoelastic behavior on speed
and loading in voluntary movements. Experimental Brain Research93, 117}180.
Milner, T.E., Cloutier, C., 1993. Compensation for mechanically unstable loading in voluntary wrist movement. Experimental Brain Research 93, 522}532.
Mussa-Ivaldi, F.A., Hogan, N., Bizzi, E., 1985. Neural, mechanical, andgeometric factors subserving arm posture in humans. Journal of
Neuroscience 5, 2732}2743.
Shadmehr, R., Mussa-Ivaldi, F.A., Bizzi, E., 1993. Postural force "eldsof the human arm and their role in generating multijoint movements. Journal of Neuroscience 13 (1), 45}62.
Stein, R.B., Kearney, R.E., 1995. Nonlinear behavior of muscles re#exesat the human ankle joint. Journal of Neurophysiology 73, 65}72.
Tsuji, T., Morasso, P.G., Goto, K., Ito, K., 1995. Human hand impedance characteristics during maintained posture. Biological Cybernetics 72, 475}485.

N

C

O
R

R

EC

O

PR

TE

The method presented in this paper requires fewer
trials to accurately estimate sti!ness during movement
than methods which have been used previously, including methods based on stochastic perturbations (Bennett
et al., 1992; Lacquaniti et al., 1993). Force perturbations,such as those used by Gomi and Kawato (1997), aresimpler to implement than our position perturbation.
However, both sti!ness and damping must be estimated,making the identi"cation more di$cult, since damping
estimates tend to be inaccurate. Therefore, our techniquefocused on sti!ness identi"cation only. As a consequence,
we are able to identify sti!ness accurately using onlya third as many trials as required by their method.
Although the neuromuscular system is nonlinear, thelinear estimation of Eq. (4) is exact if the perturbation
amplitude is small enough. In practice, the sti!ness value
depends on the perturbation amplitude (Shadmehr et al.,
1993; Milner and Cloutier, 1993). However, in contrast to
the force perturbation method (Gomi and Kawato, 1997),
this value is not biased by the perturbation direction.
Because it is not possible to realize a position perturbation briefer than 300 m s, our sti!ness estimate,measured between 120 and 180 m s after the beginningof the perturbation, includes a contribution due to re#exes. Future work will investigate the in#uence of re#exeson the sti!ness.
Once the technique introduced in this paper is implemented, it enables a simple and intuitive examination of
sti!ness during movement: sti!ness, which corresponds
to the force divided by the (constant) displacement, canbe inferred directly from the force plots. This could not bedone with methods based on force perturbations (Bennett et al., 1992; Lacquaniti et al., 1993; Gomi and
Kawato, 1997). Since our technique is also robust toD

4. Discussion

O

F(Burdet et al., 1999). When sti!ness was measured on
3 di!erent days, the sti!ness values varied by less than
10%.

U5