J Neurophysiol
88: 991–1004, 2002; 10.1152/jn.00943.2001.

Short- and Long-Term Changes in Joint Co-Contraction AssociatedWith Motor Learning as Revealed From Surface EMG
RIEKO OSU,1,2 DAVID W. FRANKLIN,2,3 HIROKO KATO,4 HIROAKI GOMI,5,6 KAZUHISA DOMEN,7
TOSHINORI YOSHIOKA,1,2 AND MITSUO KAWATO1,2
1
Kawato Dynamic Brain Project, Japan Science and Technology Corporation and 2ATR Human Information Science
Laboratories, Department 3, 2-2-2 Hikaridai, Soraku-gun, Kyoto 619-0288, Japan; 3School of Kinesiology, Simon Fraser
University, Burnaby, British V5A 1S6, Canada; 4Intelligent Communication Laboratory, Nippon Telegraph and Telephone
Corporation Communication Science Laboratories, 2– 4 Hikaridai, Soraku-gun, Kyoto 619-0237, Japan; 5Human andInformation Science Laboratory, Nippon Telegraph and Telephone Corporation Communication Science Laboratories and6Core Research for the Evolutional Science and Technology Program, Japan Science and Technology Corporation 3–1
Wakamiya, Morinosato, Atsugi-city, Kanagawa-prefecture, 243-0198, Japan; and 7Rehabilitation Center, Hyogo College ofMedicine, Nishinomiya-city, Hyogo 663-8501, Japan
Received 15 November 2001; accepted in final form 12 April 2002

Osu, Rieko, David W. Franklin, Hiroko Kato, Hiroaki Gomi,
Kazuhisa Domen, Toshinori Yoshioka, and Mitsuo Kawato.
Short- and long-term changes in joint co-contraction associated withmotor learning as revealed from surface EMG. J Neurophysiol 88:
991–1004, 2002; 10.1152/jn.00943.2001. In the field of motor control,
two hypotheses have been controversial: whether the brain acquires
internal models that generate accurate motor commands, or whether
the brain avoids this by using the viscoelasticity of musculoskeletal
system. Recent observations on relatively low stiffness during trained
movements support the existence of internal models. However, no
study has revealed the decrease in viscoelasticity associated withlearning that would imply improvement of internal models as well as
synergy between the two hypothetical mechanisms. Previously observed decreases in electromyogram (EMG) might have other explanations, such as trajectory modifications that reduce joint torques. To
circumvent such complications, we required strict trajectory control
and examined only successful trials having identical trajectory andtorque profiles. Subjects were asked to perform a hand movement in
unison with a target moving along a specified and unusual trajectory,
with shoulder and elbow in the horizontal plane at the shoulder level.
To evaluate joint viscoelasticity during the learning of this movement,
we proposed an index of muscle co-contraction around the joint(IMCJ). The IMCJ was defined as the summation of the absolutevalues of antagonistic muscle torques around the joint and computedfrom the linear relation between surface EMG and joint torque. The
IMCJ during isometric contraction, as well as during movements, wasconfirmed to correlate well with joint stiffness estimated using theconventional method, i.e., applying mechanical perturbations. Accordingly, the IMCJ during the learning of the movement was computed for each joint of each trial using estimated EMG-torque relationship. At the same time, the performance error for each trial was
specified as the root mean square of the distance between the target
and hand at each time step over the entire trajectory. The time-series
data of IMCJ and performance error were decomposed into long-term
components that showed decreases in IMCJ in accordance with learning with little change in the trajectory and short-term interactions
between the IMCJ and performance error. A cross-correlation analysis
and impulse responses both suggested that higher IMCJs follow poor

performances, and lower IMCJs follow good performances within a
few successive trials. Our results support the hypothesis that viscoelasticity contributes more when internal models are inaccurate,
while internal models contribute more after the completion of learning. It is demonstrated that the CNS regulates viscoelasticity on ashort- and long-term basis depending on performance error and finallyacquires smooth and accurate movements while maintaining stability
during the entire learning process.

Address for reprint requests: R. Osu, ATR Human Information Science
Labs, 2–2-2, Hikaridai, Seika-cho, Soraku-gun, Kyoto 619-0288, Japan(E-mail: osu@atr.co.jp).

The costs of publication of this article were defrayed in part by the paymentof page charges. The article must therefore be hereby marked ‘‘advertisement’’in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

www.jn.org

INTRODUCTION

Muscle and peripheral reflex loops possess springlike properties that pull joints back to equilibrium positions by generating restoring forces against external perturbations. This viscoelasticity can be regarded as the peripheral feedback control
gain, which is adjustable by regulating muscle co-contractionlevels and reflex gains. It has been hypothesized that by exploiting this viscoelasticity (Mussa-Ivaldi et al. 1985), the CNS
can control the limbs by simply commanding a series of stableequilibrium positions aligned along the desired movement trajectory (equilibrium-point control hypothesis) (Bizzi et al.
1984; Feldman 1966; Flanagan et al. 1993; Flash 1987; Hogan
1984). This theory, however, requires that viscoelastic forces
increase as the movement speeds up, because the dynamic
forces acting on the multijoint links grow in rough proportion
to the square of the velocity. On the other hand, the alternative
hypothesis, referred to as internal model control, enables therealization of fast and accurate movements even with low
viscoelastic forces. Under this hypothesis, the CNS learns
internal models that simulate the dynamics of the musculoskeletal system and external environment and generates the required feedforward motor commands (Bizzi and Mussa-Ivaldi
1998; Kawato et al. 1987; Miall et al. 1993; Shidara et al.
1993).
It has been a matter of controversy whether the CNS relies

0022-3077/02 $5.00 Copyright © 2002 The American Physiological Society

991

992

OSU, FRANKLIN, KATO, GOMI, DOMEN, YOSHIOKA, AND KAWATO

on high viscoelastic forces without internal models or utilizes
acquired internal models with low viscoelastic forces (Gomiand Kawato 1996; Gribble and Ostry 2000; Gribble et al. 1998;
Katayama and Kawato 1993; Koike and Kawato 1993; Lacknerand Dizio 1994; Latash and Gottlieb 1991). Recent observations on relatively low stiffness levels during well-trained
movements support the existence of internal models (Bennett
et al. 1992; Burdet et al. 2000, 2001; Gomi and Kawato 1996).
On the other hand, reports that EMG is higher in a novel
environment than a normal environment (Basmajian and De
Luca 1985; Bernstein 1967; Milner and Cloutier 1993; Thoroughman and Shadmehr 1999) indirectly suggest that the viscoelasticity at the beginning of learning may not be as low asthat after extensive training. In other words, the CNS may relyon viscoelastic forces more heavily at the beginning of learningwhen internal models are poor, and it may gradually increase
the internal model contribution as learning proceeds, resulting
in decreases in the viscoelasticity. Although several studies
have tried to model such a dual strategy (Flash and Gurevich
1997; Gribble and Ostry 2000; Katayama et al. 1998; Wang etal. 2001), no previous experimental study has clearly proventhe existence of pure decreases in viscoelastic forces thatwould imply improvement of internal models. The observed
decrease in electromyogram (EMG) during learning in previous studies might have other explanations, such as trajectory
modifications leading to reduced joint torques or lower reflex
contributions due to the attenuation of external perturbations.
One way to circumvent such complications is to require strict
trajectory control and look only at successful trials having
identical trajectory and torque profiles.
In this study, we developed a novel method to evaluate
viscoelastic forces around the joint using EMG signals andinferred changes in viscoelasticity associated with learning.
We also investigated short-term interactions between performance errors and viscoelastic forces by time-series analysis.
Our findings suggest that the relative contributions of internal
model control and viscoelasticity to the final motor command
are adaptively regulated on a long-term and short-term basis.
METHODS

Experimental design
The experiments consisted of two parts. First, we proposed an index
of muscle co-contraction around the joint (IMCJ) computed fromsurface EMG and joint torques, and compared the IMCJ with stiffness
measured using the conventional method, i.e., applying mechanicalperturbations (method evaluation). Then, we performed a second
experiment in which we elucidated learning-associated changes inviscoelasticity using the proposed IMCJ. Six healthy subjects participated in the learning experiments (4 males and 2 females; ages 20 –36
years; 1 male was left-handed). Two of the six also participated in themethod evaluation experiments under the isometric condition. Threeother subjects participated in the method evaluation experiment under
the dynamic condition (2 males and 1 female; ages 29 –34). The
institutional ethics committee approved the experiments, and thesubjects gave informed consent prior to participation.

Definition of the IMCJ
Gomi and Kawato (1996) and Burdet et al. (2000, 2001) measured
stiffness during multijoint arm movements using a high-performancecomputer-controlled mechanical interface [Parallel Link Direct-Drive

Air-Magnet Floating Manipulandum (PFM)] to displace the handslightly during each movement and measure the restoring force.
Unfortunately, these methods require many trials so they cannot be
used to observe progressive changes in the stiffness that accompanies
learning. Based on a report that the surface EMG is highly correlatedwith the static stiffness (Osu and Gomi 1999), and that the joint
stiffness is highly correlated with the joint torque (Gomi and Osu
1998; Hunter and Kearney 1982), we propose the following index forevaluating joint stiffness, using surface EMG instead of direct measurements.
If rectified surface EMG signals are assumed to be proportional toisometric muscle tension (Basmajian and De Luca 1985), the joint
torque can be expressed as the difference between the flexion torqueexerted by the flexor muscles (weighted muscle tension) and theextension torque exerted by the extensor muscles␶ s ⫽ c 1u1 ⫺ c2u2 ⫹ c5u5 ⫺ c6u6
␶ e ⫽ c 3u3 ⫺ c4u4 ⫹ c7u5 ⫺ c8u6(1)Here, ␶s and ␶e denote the shoulder joint torque and elbow joint
torque, respectively. ui denotes an individual muscle activity, which isassumed to be proportional to the rectified and averaged surface EMG
signals. u1 and u2 denote the activity of shoulder monoarticular flexorand extensor muscles, u3 and u4 denote the activity of elbow monoarticular flexor and extensor muscles, and u5 and u6 denote the activityof biarticular flexor and extensor muscles. The parameters ci include
both the moment arm and conversion factor from the muscle activity(rectified and averaged EMG) to muscle tension. The parameters ci are
all constants, as long as the moment arm is assumed to be constant.
Supposing that each muscle stiffness term is proportional to thecorresponding muscle torque (weighted muscle tension; ciui in Eq.1)(Gomi and Osu 1998; Hunter and Kearney 1982), we may be able touse the summation of muscle torques as a measure for the joint
stiffness. Therefore we define IMCJ as follows
S s ⫽ c1u1 ⫹ c2u2 ⫹ c5u5 ⫹ c6u6
S e ⫽ c3u3 ⫹ c4u4 ⫹ c7u5 ⫹ c8u6(2)Here, Ss and Se denote the IMCJ at the shoulder and elbow, respectively. Simultaneous increases in antagonistic muscle torques do notincrease the joint torque (Eq.1), but do increase the IMCJs because the
IMCJs are the summations of the absolute values of antagonisticmuscle torques (Eq.2).
At a given level of activity, muscle tension nonlinearly depends onlength and velocity. Accordingly, muscle stiffness also depends on lengthand velocity (Winters 1990), which means that surface EMG does not
linearly correlate with dynamic stiffness. Further, moment arms of somemuscles (e.g., pectoralis major, posterior deltoid, and brachioradialis)change during movements (Kuechle et al. 1997; Murray et al. 1995;
Winters 1990). However, we assumed linear length-tension, velocitytension curves, and constant moment arms, and approximated dynamic
muscle stiffness by the weighted summation of muscle activities usingparameters ci, estimated by isometric tasks. Although these simplifications were obviously wrong and might have caused error, the error was
the same for each identical movement during learning and therefore it
was possible to quantify changes in stiffness across time.
The validity of the IMCJ proposed above was assessed by thefollowing experiments under both isometric and dynamic conditions.

Evaluation of IMCJ as a good measure of stiffness
For the evaluation using isometric tasks, the subjects gripped thehandle of a force sensor and were instructed to produce a specified
force (0, 5, or 10 N) in a specified direction (16 directions in thehand’s x-y plane at even intervals) without co-contraction (Fig. 1A).
The current force vector applied by the hand to the handle and a small
cross indicating the target force were displayed on a computer mon-J Neurophysiol • VOL 88 • AUGUST 2002 • www.jn.org

CHANGES IN CO-CONTRACTION WITH MOTOR LEARNING

FIG. 1. A: experimental setup for estimating the electromyogram (EMG)torque relationship. B: experimental setup. Dotted curve connecting the start
position and target position denotes an actual trajectory, and the solid curve
denotes a target trajectory. Performance error was determined as the root meansquare of the distance between the hand and target position at each time step(the mean of the lengths of the gray lines).

itor. The right forearm of the subjects was fixed to a molded plastic
cuff tightly coupled to the handle and supported in the gravity direction by a beam. The wrist joint of the subjects was fixed by the cuff,
and only shoulder and elbow joint rotations in the horizontal plane
were permitted. During the experiment, each subject’s hand was kept
at the coordinate of [x,y] ⫽ [0.0,0.35] m. The subjects were required
to keep the head of the force vector on the target during eachexperimental set to preserve the constant external force. Additionally,
rectified and filtered surface EMG signals (moving average, 0.5 s) of
six muscles were displayed in a bar graph. A reference line wasmarked on the EMG bar graph. The reference line consisted of therectified and filtered surface EMG signals of six muscles that were
determined by requesting target force exertion before each experimental set. The subjects were also asked to keep the EMG bar graphthe same as the reference line so that the muscle activity would beconstant during each set. The stiffness was measured at the same timeby applying small perturbations using the PFM. The hand was slightly
pushed and pulled back in eight randomized directions within a brief
period (6 – 8 mm, 0.3 s, 8 directions, 3 times for each set). The subjects
were asked not to intervene voluntarily during the perturbations. The
details of the arm-impedance estimation method are provided elsewhere (Gomi and Kawato 1995, 1997; Gomi and Osu 1998).
The EMG was recorded from a shoulder monoarticular flexor(pectoralis major) and extensor (posterior deltoid), an elbow monoarticular flexor (brachioradialis) and extensor (lateral head of triceps
brachii), and a biarticular flexor (biceps brachii) and extensor (longhead of triceps brachii). The EMG signals were recorded using pairsof silver-silver chloride surface electrodes in a bipolar configuration.
Each signal was filtered [cutoff frequency, 25 Hz (low) and 1,500 Hz(high)] and sampled at 2,000 Hz. The EMG signals were rectified andaveraged for a period of 0.4 s before perturbations. This rectified andaveraged EMG was used as muscle activity ui in Eq.1. The force
exerted by the hand was measured by a force sensor attached to thehandle. The measured force was averaged for a period of 0.4 s. The
joint torque was calculated from the average force using a Jacobianmatrix.
The joint torque was decomposed into muscle torques using parameters estimated from a linear regression of the joint torque andmeasured EMG ui (Eq. 1). All of the parameters ci were estimated by
the least-square-error method. Then, each IMCJ was computed bysumming all the absolute values of flexor muscle torques and theabsolute values of extensor muscle torques related to each joint(Eq. 2). The computed IMCJs were compared with the stiffness
measured at the same time by applying perturbations (PFM-measured

993

stiffness). We estimated linear relationships between the PFM-measured stiffness ([Nm/rad]) and the IMCJs ([Nm]) using the acquired
data so we could convert the IMCJ unit to stiffness unit (rIMCJ).
Because the approximation of dynamic stiffness using IMCJ relies
on oversimplified assumptions such as constant moment arms, linear
length-tension, and velocity-tension curves, we evaluated whether
IMCJ is still applicable for the movement data despite these simplifications. We compared IMCJ, computed by applying isometrictorque-EMG relationship to dynamic EMG, to dynamic stiffness
measured simultaneously. To confirm that the IMCJ represents dynamic stiffness, it is important to obtain a wide variety of stiffness
estimates. However, estimating a single dynamic stiffness requires
many more trials than estimating a static stiffness. To acquire enoughvariety of stiffness values, we used data from three subjects acrossthree tasks measured on different days. Each subject learned threedifferent force-fields (null force-field, velocity-dependent force-field,
and position-dependent force-field) on different days. Before learning
each force-field, the isometric torque-EMG relationships were measured to estimate parameters ci in the same way as described above
except that no perturbation was given. Then, subjects performed
horizontal point-to-point movement away from the body in one of thethree force-fields. After enough training, stiffness during the movements was measured by applying small positional perturbation to thehand (Burdet et al. 2000). At the same time, EMG signals wererecorded from six arm muscles, and the corresponding shoulder andelbow IMCJs were computed using the estimated parameters ci. Eachdynamic stiffness was compared with the corresponding IMCJ. See
APPENDIX 2 for details.

Learning experiments
Six subjects participated in the learning experiments. The learning
experiments themselves consisted of three parts. Prior to the learningtask, the subjects executed isometric contraction tasks, enabling us toestimate the relationship between the surface EMG and joint torque(parameters ci in Eq.1) for the calculation of the IMCJs. Then, thesubjects learned reaching movements under strict trajectory control.
After the learning task, the subjects executed isometric contraction
tasks again to confirm that the state of the electrodes after the learning
was not different from that before the learning. From this, we confirmed that the electrode interface was not responsible for the observed changes in the surface EMG.
In the isometric contraction tasks prior to and following the learning task, each subject’s hand was coupled to a force-torque sensor.
The subject was instructed to produce a specified force (0, 5, 10, or 15
N for the prelearning trials and 0 or 10 N for the postlearning trials)
in a specified direction (16 directions). The hand position and the
instructions were the same as in the method evaluation experiment
except that no perturbation was given. The hand force and EMG
signals were recorded in the same way as in the method evaluation
experiments.
In the learning task, the subjects performed reaching movementswith the shoulder and elbow in the horizontal plane at the shoulder
level. Wrist movements were constrained by a brace. The learning
task consisted of moving the hand in unison with a target moving
along a specified trajectory (Fig. 1B). The specified trajectory was
curved inward, which was opposite to the natural curvature of spontaneous movements (Nakano et al. 1999). The average of 20 trials of
the subject’s own hand trajectory conforming to a 3.5-cm-wide inwardly curved path, performed during practice, was shifted 5 cm away
and used as the target trajectory (no constraint on time was given
during these preparatory 20 trials). The hand start position and end
position, which were shifted away from the original average trajectory, were located at [x,y] ⫽ [⫺0.24,0.37] m and at [x,y] ⫽ [0.21,0.39]m, respectively, i.e., a movement of about 0.5 m performed in approximately 0.5 s.
Obviously, the execution of such an unusual trajectory with a strict

J Neurophysiol • VOL 88 • AUGUST 2002 • www.jn.org

994

OSU, FRANKLIN, KATO, GOMI, DOMEN, YOSHIOKA, AND KAWATO

time course requires learning. The applied shift from the original
average trajectory also enhanced the requirement of learning. To
ensure that the subject would learn the accurate geometry and timecourse of the target trajectory, only the hand trajectories close to thetarget trajectory (⬍4 cm at each time step) were regarded as successful trials. These operations enabled us to acquire movements having
identical trajectory and torque profiles both at the beginning and at theend of learning, which was necessary to prove the existence of puredecreases in the viscoelasticity and implies improvement of internal
models. The current hand and target positions were displayed on a
CRT. After each trial, feedback of the resulting movement wasprovided to the subject by replaying the target and hand movementson the CRT, providing temporal and positional error information.
Hand positions within 4 cm or over 4 cm from the target weredisplayed in different colors so that the subject could learn his/her
weak points. The performance error for each trial was specified as theroot mean square of the distance between the target position and theactual hand position at each time step over the entire trajectory. The
number of trials required was 96 –120. Position data were acquired
using the OPTOTRAK system. Surface EMG signals were recorded
from six muscles involved in shoulder and elbow movements in thesame way as in the method evaluation experiment.

Computation of dynamic torque
The position data obtained during the learning trials were digitallyfiltered by a fourth-order Butterworth filter with an upper cutoff
frequency of 10 Hz. Derivatives of the position data were calculated
by successively applying a three-point local polynomial approximation. Ballistic components of the movements were extracted using thecurvature as a threshold to determine the beginning and the end ofeach movement [500 (1/m)] (Pollick and Ishimura 1996). Dynamic
torques were calculated through the dynamics equation of a two-joint
arm model using the position data and link parameters estimated fromthe link length for each subject (the data of an adult man’s arm
measured with a 3-dimensional scanner as a standard). The mass of
the links was adjusted for each subject by changing the standard value
proportional to the link length of the subject. The inertia moment of
the links was adjusted by changing the standard value proportional to
the third power of the link length of the subject. Viscosity coefficientswere estimated from the absolute average torque for each movement
using the equation in Gomi and Osu (1998). We averaged the absolutedynamic torques across whole movement durations determined bycurvature criteria (average dynamic torque).

IMCJ during learning
Since IMCJs were correlated with the PFM-measured stiffness
during isometric contraction tasks and dynamic tasks (see RESULTS),
we calculated the IMCJ during the learning of movements. First,
parameters ci in Eq.1 were estimated for each subject from the EMG
signals and joint torques in the isometric contraction tasks executed
before the learning of the movements. The EMG signals during the
movements were rectified and averaged across the entire movement
duration determined by curvature criteria for each muscle. Further, toroughly examine which part of each movement duration is responsible
for the change, each movement duration determined by the curvature
criteria was divided into the first half and the latter half, and the EMG
signals were rectified and averaged over either the first half or thelatter half of each movement duration. Then, the estimated parameters
ci were applied to the rectified and averaged EMG signals during
movements ui to compute the average torques of individual muscles(ciui). The IMCJs of the shoulder and elbow were computed as thesummations of the average absolute torques of individual muscles
according to Eq. 2.

Because the parameters ci were computed from isometric data and do
not take into account the changing moment arm or velocity-tension
relation, they do not accurately reproduce muscle torque during movements. The EMG signals required to generate certain muscle torques
were larger during the movements than during the isometric contraction,
probably due to muscle tension-shortening-velocity characteristics. Accordingly, the IMCJs during the movements computed from isometric
EMG-torque relationships were corrected according to EMG-torque relationships during the movements. Namely, IMCJs were scaled based on
the ratio of dynamic torque to EMG-estimated torque. The dynamic
torque was, as described above, computed from actual movement trajectories using an arm inverse dynamics model and then rectified andaveraged (average dynamic torques). The EMG-estimated torque wascomputed according to Eq.1 by applying the isometric EMG-torquerelationship (ci) to the EMG during the movements (ui) and rectified andaveraged in the same way as the absolute dynamic torque. The correction
for the movements from the isometric condition was made for each jointof each subject. We further re-scaled IMCJ by converting the unit of
IMCJ ([Nm]) to the unit of stiffness ([Nm/rad]), using a linear relationship between the PFM-measured stiffness and the IMCJ estimated in themethod evaluation experiment (rIMCJ).
No matter how strictly we constrained the movements, the joint
torque might have slightly differed from trial to trial, which might
contribute to the stiffness values (Gomi and Osu 1998). To extract therIMCJ independent of the joint torque that implies improvement ofinternal models, we subtracted the torque-dependent components fromthe rIMCJ. Assuming that the average rIMCJ was linearly dependent
on the average dynamic torque, we expressed the average rIMCJ as
the summation of the weighted average dynamic torque, a constant,and residuals that could not be explained by the joint torque. The
parameters (the weight and the constant) were linearly estimated by
the least square error method. Then, the torque dependent components
were subtracted from the rIMCJ. We called this residual component
the torque-independent rIMCJ.

Bayesian multivariate feedback model for statistical analysisof time-series data
Progressive changes in viscoelasticity and performance error may bedescribed by a dynamical system with stochastic noises. To examine theproperties of the system, we can apply time-series analysis to the rIMCJ
and performance error during learning. To draw inferences from thetime-series data, we need to select a suitable hypothetical model torepresent the data. Having chosen a model, it becomes possible toestimate parameters and use the fitted model to enhance our understanding of the mechanism generating the series. Accordingly, we set up astatistical model whose structure was designed assuming interactions
between viscoelastic force and performance error (Fig. 2).
The observed rIMCJ and performance error at a certain trial number
were assumed to consist of the following three components: 1) asmooth and long-term change in the mean level of the rIMCJ andperformance error, expressing a gradual decrease with the progress oflearning (nonstationary trend components); 2) short-term fluctuating
components depending on previous trials, describing interactions between rIMCJ and performance error and able to be expressed as anauto-regressive (AR) model (cyclical components); and 3) observation noise. The rIMCJ and performance error can be expressed asfollows
R共n兲 ⫽ r共n兲 ⫹ tr共n兲 ⫹ ␻r共n兲 ␻r ⬃ N共0,10⫺4兲(3)E共n兲 ⫽ e共n兲 ⫹ te共n兲 ⫹ ␻e共n兲 ␻e ⬃ N共0,10⫺4兲(4)Here, n (1 ⱕ n ⱕ N) denotes the current trial number. Current rIMCJ
R(n) is composed of cyclical component r(n), trend component tr(n),and observation noise ␻r(n). Current performance error E(n) is composed of cyclical component e(n), trend component te(n), and obser-

J Neurophysiol • VOL 88 • AUGUST 2002 • www.jn.org

CHANGES IN CO-CONTRACTION WITH MOTOR LEARNING

995

which were themselves generated based on the estimated model(Monte Carlo simulation, see APPENDIX 1).
RESULTSHigh correlation between measured stiffness and IMCJ

FIG. 2. Statistical model whose structure was designed assuming interactions between viscoelasticity and performance error. See the text for details.

vation noise ␻e(n). N(m, ␴ 2) is the normal distribution with mean mand variance ␴2. Trend components are modeled in the form of thefollowing second-order stochastic difference equations
t r 共n兲 ⫽ 2 䡠 tr共n ⫺ 1兲 ⫺ tr共n ⫺ 2兲 ⫹ ␧tr共n兲␧tr ⬃ N共0, ␴tr2 兲(5)t e 共n兲 ⫽ 2 䡠 te共n ⫺ 1兲 ⫺ te共n ⫺ 2兲 ⫹ ␧te共n兲␧te ⬃ N共0, ␴te2 兲(6)Here, ␧tr(n) and ␧te(n) denote system noise. The order of the two is
selected to extract low-frequency components as the mean-nonstationary trends.
The following information regarding the short-term interaction
between the rIMCJ and performance error was incorporated into themodel describing cyclical components: 1) current rIMCJ may change
according to the performance errors in the previous trials (feedback
from performance error to rIMCJ), and 2) A high rIMCJ is assumed
to decrease the current performance error but is unlikely to have aneffect on subsequent performance errors (instantaneous response of
performance error to rIMCJ). Therefore the cyclical components in
Eqs.3 and 4 are described by the following special form of a multivariate auto-regressive model allowing instantaneous responses

冘M

r共n兲 ⫽

冘

We confirmed that the IMCJs were actually linearly correlated with the stiffness measured directly by applying mechanical perturbations (PFM-measured stiffness) during isometricforce regulation tasks. We first confirmed that the joint torquecould be linearly reconstructed from EMG signals in accordance with our previous studies (Gomi and Osu 1998; Osu andGomi 1999). Figure 3A compares the measured joint torqueand joint torque reconstructed from EMG signals of subject
H.S. The coefficients of determination for the two subjects
were 0.976 and 0.971. Therefore, first of all, the joint torque
was well predicted from EMG signals. Then, IMCJs werecomputed as the summations of absolute torques of individual
muscles (Eq.2). Figure 3B shows the relationship between the
PFM-measured joint stiffness and corresponding IMCJ for
both subjects. The open circles denote shoulder joint stiffness
and the crosses denote elbow joint stiffness. The thin marks
denote data from subject Y.K. and the thick marks denote data

M

Arr共m兲r共n ⫺ m兲 ⫹m⫽1

Are共m兲e共n ⫺ m兲 ⫹ ␧r共n兲

␧r ⬃ N共0, ␴r2兲

(7)

m⫽1

冘
M

e共n兲 ⫽ Q 䡠 r共n兲 ⫹Aee共m兲e共n ⫺ m兲 ⫹ ␧e共n兲

␧e ⬃ N共0, ␴e2兲

(8)

m⫽1

Here, r(n) is assumed to depend on the rIMCJ, performance errors of
the previous M trials [r(n ⫺ m), e(n ⫺ m), 1 ⱕ m ⱕ M], and system
noise ␧r(n). e(n) is assumed to decrease with current rIMCJ [Q ⫻ r(n),Q ⬍ 0] and to depend on the performance errors of the previous M
trials [e(n ⫺ m), 1 ⱕ m ⱕ M] and system noise ␧e(n).
This model is formulated as an extended Bayesian multivariate
feedback (BMF) model (Kato and Kawahara 1998; APPENDIX 1). For
the system analysis, the special form of the cyclical components wastransformed into an ordinary form of a multivariate AR (MAR) model
under the assumption that the noise sequences are mutually independent. This model can be represented in state space form and a Kalmanfilter algorithm can be applied to calculate the likelihood of the model.
Parameters Arr, Are, Aee, Q, ␴r, ␴e, ␴tr, and ␴te were estimated by the
maximum likelihood method for each subject (Ishiguro and Akaike1989). The order M of the model was selected by Akaike’s Information Criterion (AIC) (Akaike 1974).
The impulse responses of the system enable us to describe how theperformance error and rIMCJ interact with each other. The responseof the rIMCJ obtained by providing a unit impulse input to a performance error reveals how the CNS utilizes the information of previous
performance levels to modify a subsequent rIMCJ. The impulse
responses of the model were calculated based on estimated model
parameters (Akaike and Nakagawa 1972; Ishiguro et al. 1999).
The reliability of the estimated parameters of the model was confirmed by reapplying the system analysis to the simulated data sets,

FIG. 3. A: regression results for subject H.S. Horizontal axis denotes thejoint torque predicted from Eq. 1, and the vertical axis denotes the measured
joint torque. Left: shoulder torque; right: elbow torque. B: relationship between
the joint stiffness of the shoulder and elbow measured using perturbations(horizontal axis) and the corresponding index of muscle co-contraction aroundthe joint (IMCJ) (vertical axis) for subjects Y.K. and H.S. The open circles
denote the shoulder joint stiffness and the crosses denote the elbow joint
stiffness.

J Neurophysiol • VOL 88 • AUGUST 2002 • www.jn.org

996

OSU, FRANKLIN, KATO, GOMI, DOMEN, YOSHIOKA, AND KAWATO

from subject H.S. A linear relationship was observed between
the IMCJs and PFM-measured joint stiffness. The correlation
coefficients for the two subjects were 0.891 and 0.882. These
results suggested that the IMCJ can closely predict the magnitude of the joint stiffness. The following linear relationshipbetween the IMCJs and joint stiffness was estimated by using
the least square error method. Because the relationships between the IMCJs and PFM-measured joint stiffness were similar for both the shoulder and elbow of both subjects, the slope
and the intercept were estimated using all of the datarIMCJ 共关Nm/rad兴兲 ⫽ 6.27 䡠 IMCJ 共关Nm兴兲 ⫹ 4.61

IMCJs ([Nm]) during learning movements were converted
into rIMCJs ([Nm/rad]) using the above linear relationship.
Note that the IMCJs were corrected for movements from theisometric condition for each joint of each subject before theconversion to rIMCJs, as explained in METHODS.
We also confirmed that the IMCJ is applicable for the
dynamic condition. Figure 4 compares joint stiffness and thecorresponding IMCJ. Each asterisk represents a value from oneof the three subjects under one of the three force-fields. Even
across subjects and tasks, we still observed good linear relationships between IMCJ and dynamic stiffness (r ⫽ 0.85 forshoulder and 0.78 for elbow). We may suppose that, if themeasurements were limited to the same day and the same
subjects, the reliability of the IMCJ would be even better than
the results obtained here. At least within similar movement
trajectories, the current method works well to quantify the
relative change of joint co-contraction during the dynamic
condition.

FIG. 5. Changes in the profiles of the joint torque and rIMCJ for subject
Y.M. From top to bottom are the shoulder and elbow torques (␶s, ␶e) andshoulder and elbow rIMCJs (Ss, Se). From left to right are the initial 4 trials (allunsuccessful), the initial 4 successful trials, and the final 4 successful trials.

Invariant EMG-torque relationships before and after
learning
As each experiment took a few hours, the impedance of theelectrode interface might have changed with the passing oftime. To verify that changes in the electrode interface were notresponsible for observed changes in the surface EMG, we
compared the relationships between the EMG levels andtorques in the isometric contraction tasks prior to learning withthose after learning for each subject. If the EMG levels inexerting the same joint torques were considerably lower after
learning than before learning, the decreases in the EMG observed during the learning could not be ascribed to the effectsof the learning. Fortunately, the observed EMG-torque relationships after learning were not detectably different fromthose before learning for any subject. The joint torque could bereconstructed from the EMG in isometric contraction tasks
after learning by using the parameters estimated from the
EMG-torque relationships before learning. The coefficient ofdetermination for the seven subjects was 0.923 ⫾ 0.059 (SD).
The slope of the regression line for the seven subjects was
1.008 ⫾ 0.099. The high coefficients of determination and theslope values close to one suggested that the relationships
between the torques and EMG levels were preserved even after
extensive trials. Therefore the observed decreases in the EMG
levels after learning could not be ascribed to the long-termchanges of the electrode state or muscle fatigue.
Long-term decrease of rIMCJ

FIG. 4. A and B: relationship between the IMCJ of the shoulder (A) andelbow (B) (horizontal axis) and the corresponding joint stiffness measured
using perturbations (vertical axis) during movements.

Figure 5 shows the changes in the joint torque and rIMCJ
time profiles during learning for subject Y.M. The first andsecond rows show the shoulder and elbow torques, respectively. The torques were calculated using the dynamics equation of a two-joint arm model. The third and fourth rows show
shoulder and elbow moving-averaged rIMCJ respectively. The
moving-averaged rIMCJ was calculated by applying estimatedparameters ci to EMG signals that were rectified and averaged
using a 0.1-s moving-average window. The left column shows

J Neurophysiol • VOL 88 • AUGUST 2002 • www.jn.org

CHANGES IN CO-CONTRACTION WITH MOTOR LEARNING

997

the profiles of the initial four trials, which were all unsuccessful; the middle column shows the profiles of the initial four
successful trials (early stage of learning); and the right column
shows the profiles of the final four successful trials (late stageof learning). At the very beginning of the learning, the subjects
failed to meet the task requirements, and as a result, theirtorque profiles were variable from trial to trial. However, evenat the early stage of learning, the subjects soon managed toachieve several successful trials. The torque profiles of these
initial successful trials were nearly identical to the torqueprofiles of the final successful trials. The applied strict constraint on the trajectory worked well to acquire rIMCJ datawith similar torque profiles. As shown by the profiles, therIMCJ in the successful trials decreased although the changes
in the torque profiles were small. The decreases in the rIMCJ
were more evident in the latter half of each movement duration.

FIG. 7. A–C: (A) Shoulder and (B) elbow rIMCJs and (C) performance error
across successful trials normalized for each subject (z score: mean 0, SD 1). Eachcolor corresponds to a subject (subject Y.O., magenta; subject K.D., green; subject
N.H., cyan; subject H.S., blue; subject Y.K., yellow; subject Y.M., red).

FIG. 6. A–C: (A) Shoulder and (B) elbow rIMCJs and (C) performanceerror across all trials for each subject. Each color corresponds to a subject(subject Y.O., magenta; subject K.D., green; subject N.H., cyan; subject H.S.,
blue; subject Y.K., yellow; subject Y.M., red).

Figure 6, A–C, shows changes in shoulder and elbow rIMCJ
and performance errors accompanying learning, averaged
across the entire movement duration. Each color corresponds
to a subject (Y.O., magenta; K.D., green; N.H., cyan; H.S.,
blue; Y.K., yellow; Y.M., red). The solid curves denote trend
components extracted by applying the second-order trendcomponent model expressed in Eqs. 5 and 6. Both the rIMCJ
and performance errors fluctuated in the short term. In four offive subjects, the mean levels of rIMCJ for both shoulder andelbow gradually decreased as the learning proceeded. For theother two subjects (N.H. and Y.K.), the mean levels of rIMCJ
were rather low at the beginning of learning, but at the sametime, the performance errors were comparatively large. For
these subjects, the hand fell short of the target at the beginningof the learning, resulting in smaller rIMCJs with larger performance errors.
Figure 7, A and B, shows relative changes in shoulder andelbow rIMCJ in successful trials (dotted curves), with thesuperimposed trends (solid curves) extracted by applying thesecond-order trend component model. Because the magnitude
of rIMCJ was different between subjects, it was normalized foreach subject to 0 mean and 1 SD. Each color corresponds to a
subject. Figure 7C shows performance errors in successful
trials (dotted curves) with the trends (solid curves). As ex-

J Neurophysiol • VOL 88 • AUGUST 2002 • www.jn.org

998

OSU, FRANKLIN, KATO, GOMI, DOMEN, YOSHIOKA, AND KAWATO

pected, the performance errors were small and nearly constantand identical for all successful trials. Because of individual
differences in performance, the number of successful trials
were different between subjects. The success rate of eachsubject was 70% for Y.O., 81% for K.D., 54% for N.H., 79%
for H.S., 48% for Y.K., and 35% for Y.M. In these successful
trials, every subject showed a negative correlation between thetrial number and rIMCJ except for one subject (Y.K., denotedin yellow). Even for subject N.H. (denoted in blue), whoshowed an increase in rIMCJ at the beginning of the learning,rIMCJ for successful trials showed a significant decrease inboth shoulder and elbow. The results suggest that the rIMCJ
decreased with little change in the performance error.
Table 1 shows correlation coefficients between the successful trial number and the average torque-independent rIMCJ(see METHODS for definition). For the shoulder joint, five of six
subjects showed a significant negative correlation between thetrial number and rIMCJ. For the elbow joint, four of six
subjects showed a significant negative correlation between thetrial number and rIMCJ. Therefore, for the majority of thesubjects and joints, the rIMCJ was found to decrease withlearning independent of the joint torque. Consequently, learning enabled the generation of similar trajectories with lesscontribution from the viscoelasticity.
Short-term interaction between rIMCJ and performanceerror
As shown in Fig. 6, both the rIMCJ and performance error
went up and down frequently throughout the learning process.
Looking at the figure, these fluctuations appear to be simple
white noise, but they might have some temporal interactions.
They might reflect some dynamic system underlying the learning process. To examine the temporal relationship between therIMCJ and performance error, the cross-correlation was computed according to the following deterministic cross-correlation sequence. The values were normalized so that, for an
autocorrelation, the sample at zero lag would be 1.0

冘

N⫺兩m兩

C er 共m兲 ⫽e共n兲r共n ⫹ m兲

(9)FIG. 8. Average cross-correlations between performance error and rIMCJ andtheir 95% confidence intervals, in the latter half of each movement duration,
averaged across 6 subjects. Trend components were presubtracted from normalized values. A positive correlation value for a positive time lag means a positive
correlation between the performance error and subsequent rIMCJ.

rIMCJ in the latter half of each movement duration, averaged
across all subjects, and their 95% confidence intervals. Significant positive-correlation values were observed at the lags of
one and two trials for the shoulder, and at the lag of two trials
for the elbow. A positive correlation value for a positive time
lag means a positive correlation between the performance error
and subsequent rIMCJ, that is, a higher rIMCJ follows a poor
performance, and a lower rIMCJ follows a good performance.
For the elbow, we also found a significant negative correlation
at zero lag, suggesting an instantaneous response of the performance error to the rIMCJ. The negative correlation values
for zero lag indicated that the performance was good when therIMCJ was high, while the performance was bad when therIMCJ was low. The observed positive correlation values atlags of a few trials suggested the existence of short-term
interactions between the viscoelasticity and performance error.
These observations were found in the latter half of the movement duration but were not significant in the first half of the
movement duration.
To further examine the properties of the interactions, weanalyzed the time series data for each subject using the BMF
statistical model described in Eqs. 3– 8 and in Fig. 2. Theinteractions between rIMCJ and performance error were examined for each subject by computing the impulse response ofn⫽1Here, N denotes the total number of the trials, and m denotes
the lag in the trial number between the rIMCJ and performanceerror. e(n) and r(n ⫹ m) denote the cyclical components of theperformance error and rIMCJ normalized for each subject. The
long-term trend components described by Eqs.5 and 6 were
presubtracted from the performance error and rIMCJ. Figure 8
shows cross-correlations between performance error and
TABLE 1.

Correlation coefficients between the successful trial
number and rIMCJ for each subject
Subject

Shoulder

Elbow

Y.O.
K.D.
N.H.
H.S.
Y.K.
Y.M.

⫺0.709**⫺0.660**⫺0.415**⫺0.442**0.175 n.s.
⫺0.418**⫺0.644**⫺0.486**⫺0.378**0.071 n.s.
0.075 n.s.
⫺0.520**** P ⬍ 0.01.
n.s., not significant: rIMCJ, index of muscle co-contraction around the joint.

FIG. 9. Estimated response of rIMCJ when a performance error impulse was
fed into the system for each subject. Solid curves denote estimates obtained fromreal data. Broken curves denote estimates obtained from simulated data generated
by the Monte Carlo method (Marsaglia and Bray 1964). Both the real data andsimulated data showed similar positive impulse responses for 5 of 6 subjects(subjects Y.O., K.D., N.H., H.S., and Y.K.), suggesting that the estimated performance error–rIMCJ interactions were reliable.

J Neurophysiol • VOL 88 • AUGUST 2002 • www.jn.org

CHANGES IN CO-CONTRACTION WITH MOTOR LEARNINGTABLE 2.

Maximum likelihood and AIC values for each AR orderin cyclical components of measured and simulated data for subject
Y.O.

999

and for subject N.H., it was two. The reliability of the statisticalmodel is described in APPENDIX 1.
DISCUSSION

Measured Data

IMCJ

AR OrderMaximum Likelihood

AIC

1
2
3
4
5

⫺173.185⫺171.047⫺167.591⫺163.263⫺162.038

362.369
364.094363.182360.527364.075

Simulated Data
AR OrderMaximum Likelihood

AIC

1
2
3
4
5

⫺169.350⫺165.547⫺150.505⫺143.564⫺143.249

354.700
353.094
329.009321.127326.499

AIC, Akaike’s Information Criterion; AR, auto-regressive.

the system. To confirm the reproducibility of the model, we
executed simulations by the Monte Carlo method (see APPENDIX
1). We processed only the rIMCJ values for the latter half ofeach movement duration, where significant cross-correlations
were found.
Figure 9 shows the estimated response of the shoulder
rIMCJ when a performance error impulse was input into thesystem. For five of six subjects, similar impulse responses wereobtained from both real data and simulated data generated by
the Monte Carlo method, suggesting that the observed transfer
characteristics between the components were significant. Thatis, for these subjects, the shoulder rIMCJ showed a significant
positive response to performance error input, which indicatesthat rIMCJ was positively correlated with the performanceerror levels in the immediately preceding trials. The results
suggested that large performance errors will lead to a greater
shoulder viscoelasticity, whereas small performance errors willlead to a lower viscoelasticity in subsequent trials. Interactions
like these were not so obvious for the elbow rIMCJ and when
the shoulder rIMCJ was restricted to the first half of eachmovement duration. Therefore previous performance mainly
affected the viscoelasticity during the braking phase, and theshoulder viscoelasticity was more sensitive to this preceding
performance than the elbow viscoelasticity.
The reproducibility of the estimated models and the reliabilityof the estimated values were examined from several aspects.
Table 2 shows maximum likelihood and AIC values for each ARorder in the cyclical components for subject Y.O. The top tableshows the values when real data were used, while the bottom tableshows the values when simulated data generated by the Monte
Carlo method were used (see APPENDIX 1). In the case of thissubject, an order of four (where the AIC value was minimal) was
selected in the real data. The AIC value was again minimal at theorder of four in the simulated data. The same AR order selected
for both real data and simulated data demonstrates the reproducibility of an estimated model. All of the subjects here showed thesame AR order for both the real data and simulated data. The ARorder selected for subjects K.D., H.S., Y.K., and Y.M. was one,

We proposed a novel method for evaluating muscle cocontraction levels around the joint using surface EMG signals,
and we confirmed that the values correspond well to the joint
stiffness values in isometric contraction tasks as well as in
dynamic tasks. The advantage of this method over the directcomparison of raw EMG signals, or EMG signals normalizedby maximum voluntary contraction (MVC), is that the values
have a physically meaningful unit. Because the magnitudes of
raw EMG signals will change drastically for a number ofreasons, such as the state of the electrodes, the distance between the electrodes, the configuration of muscles, skin condition, etc., they do not directly relate to a physical quantity,
such as a force or stiffness. Even with MVC normalization, thecontribution of each muscle is hard to assess. Therefore adding
the raw EMG or normalized EMG of different muscles together has no physical meaning, because the weights among
the multiple muscles are quite arbitrary. In contrast, in thepresent method, EMG signals are converted with reference tothe generated joint torque so that they represent an absolute
quantity. Because these converted values successfully represent a physically meaningful unit, they enable arithmetic operations such as addition and subtraction among operations ofmultiple muscles. This gives a measure for the net joint stiffness composed of multiple muscles as well as the relative
contribution of each muscle to the joint stiffness. The present
method can be used in a practical manner because the only
requirement for computing the IMCJ is to measure the joint
torque and surface EMG signals. Measuring the relationshipbetween the joint torque and EMG is rather easy comparedwith that between the stiffness and EMG in terms of bothmechanical and procedural demands (Osu and Gomi 1999).
Although the EMG-stiffness relationship during movements
is quantitatively different from that during isometric contraction because of nonlinear components (such as length-tension
or velocity-tension curves) or changes in moment arms according to the posture, the IMCJ successfully reproduced the relative change of joint co-contraction during movements (Fig. 4).
The linear assumptions for force/velocity curves might haveresulted in overestimating the agonist muscle torques compared with antagonist muscle torques. However, small errors inweighting antagonistic muscle pairs would not have severely
affected the obtained results because none of the individual
muscles showed a tendency distinctively opposite to the observed trend during learning. The assumption of constant moment arms is also over-simplified, especially for shoulder muscles. Pectoralis major moment arm might have changed 20%
and posterior deltoid moment arm might have changed 30% for
the movements examined here (Kuechle et al. 1997). Biceps
and triceps moment arms seemed to be rather constant, while
the brachioradialis moment arm might have changed 20 –30%
for the movements examined here (Murray et al. 1995). Suchchanges might have caused some error but would not haveseverely affected the results as long as we were observing
almost identical trajectories within each subject, on the sameday, and just looking at relative changes.

J Neurophysiol • VOL 88 • AUGUST 2002 • www.jn.org

1000

OSU, FRANKLIN, KATO, GOMI, DOMEN, YOSHIOKA, AND KAWATO

Even if the linearity assumption was wildly violated inreality, the proposed rather simple model and methods are still
useful for the purpose of comparison and give us a better index
for net joint stiffness than conventional normalization andsummed EMG from multiple muscles. Previous estimation of
dynamic torque and stiffness from EMG using a neural network model (Koike and Kawato 1995), and our recent estimation of dynamic stiffness from EMG signals with simultaneous
measurements of stiffness by PFM, support the validity ofusing EMG signals for the stiffness estimation (Franklin et al.
2000). In summary, if torque levels are well corrected formovements, the present method is quite effective in assessing
relative changes in the joint viscoelasticity during movements
within a single subject.
Integration of internal model control and viscoelasticity
We observed a long-term decrease in rIMCJ in the learningof planar arm movements. The learning enabled the subjects togenerate the same trajectory with less rIMCJ. The present work
succeeded in observing pure changes in the viscoelastic force
by requiring strict trajectory control and examining only successful trials. Such observations would be difficult in forcefield learning (Milner and Cloutier 1993; Shadmehr andMussa-Ivaldi 1994; Thoroughman and Shadmehr 1999), because force-fields severely perturb hand trajectories. Thoroughman and Shadmehr (1999) observed decreases in wasted
contraction, that is, the amount of activation canceled by opposing muscles, in the learning of a novel force-field. However, wasted contraction does not necessarily reflect the viscoelasticity of the system. Even with a decrease in the wasted
contraction, the increased effective contraction that compensates for the applied force-field can provide the necessary
viscoelasticity to stabilize the movement. Accordingly, theirfindings do not directly demonstrate the change of weight on
the viscoelastic force. By giving strict constraints on trajectories and comparing stiffness values among movements withsimilar trajectories, we were able to demonstrate the decreaseof viscoelastic force for the first time.
Given that there is no change in the movement trajectory nor
in the joint torque, the only possible explanation for the decrease in the viscoelastic forces (peripheral feedback controlgain) is an increase in the contribution of the feedforward
component (internal model control). When learning a new task,
unpredictable dynamic forces acting on the multijoint links(e.g., interaction force, Coriolis force, or centrifugal force)might perturb the movements. As implicated by the equilibrium-point control hypothesis, viscoelastic force produced byboth intrinsic muscle properties and spinal reflex will counteract such unpredictable perturbation (Flash 1987). As the learning proceeds, the CNS acquires the internal models that predictand generate necessary motor commands to compensate for the
perturbation (Imamizu et al. 2000). The CNS weighs the viscoelastic force more strongly at the beginning of learning when
the internal models are poor, and it gradually increases theinternal model contribution as the learning proceeds. The
short-term interaction between the performance error and
rIMCJ suggests that the CNS actively regulates the viscoelasticity.
We propose an integration of the two theories, that is, theequilibrium-point control hypothesis and internal model– con-

FIG. 10. A: integration of internal model control and viscoelasticity-dependent control. Thick green arrow shows that a high viscoelastic force instantaneously reduces the performance error. Thick red arrow shows that the performance error results in transient increases in viscoelastic forces insubsequent movements. Thick blue arrow shows that an improvement of theinternal models results in gradual decreases in both viscoelastic forces andperformance error. Summation of the output from these 2 control systems is the
final motor command sent to the limbs. B: changes rIMCJ (green), performance error (red), and accuracy of the internal models (blue) associated withlearning predicted by the integrated model. As the internal models improve, theoutput from the internal models increases.

trol hypothesis, on the assumption that learning improves theinternal models (Fig. 10A). The final motor command is thesummation of the parallel outputs from the “feedforward controller using the internal models” and “feedback controller
supported by the viscoelasticity”; furthermore, the contributionof each output to the command is regulated by changing theperipheral feedback gain, that is, the magnitude of the viscoelastic force. The CNS monitors the improvement of theinternal models and the performance. The CNS relies on theviscoelasticity when the internal models are imperfect or theenvironment is unstable, while it utilizes the internal models
after they improve and the environment is stable (Fig. 10B). In
addition to such long-term interactions, our findings further
suggest short-term viscoelasticity-performance interactions ona trial-by-trial basis. If a movement is currently inaccurate, thecontribution of the viscoelasticity is increased within several
trials to improve the performance. If a movement is currently
accurate, on the other hand, the contribution of the internal
models is increased within several trials. Accordingly, the CNS
can learn the internal models without loss of movement accuracy by using viscoelasticity-dependent control at the onset oflearning and employing off-line feedback of the performance
to regulate the viscoelasticity.
Such a strategy can be related to the common engineering
technique of solving a problem by slowly shifting task parameters from domains where the solution is easy to domains
where the solution is difficult (continuation method). In robot

J Neurophysiol • VOL 88 • AUGUST 2002 • www.jn.org

CHANGES IN CO-CONTRACTION WITH MOTOR LEARNING

learning, for example, higher learning rates were obtained by
gradually increasing the speed or decreasing the stiffness as thelearning progressed (Katayama et al. 1998; Sanger 1994).
Behavioral studies have suggested that the internal models inthe brain are relatively local in their generalization. They areaccurate around previously learned trajectories but inaccurate
for unexplored trajectories; the internal models are not parametric with a global generalization capacity (Gandolfo et al.
1996; Ghahramani et al. 1996; Imamizu et al. 1995; Kitazawa
et al. 1997; Thoroughman and Shadmehr 2000). An internal
model that achieves a novel desired trajectory cannot be acquired from those trajectories that differ greatly from thedesired one. By sufficiently increasing the viscoelastic forces,
however, trajectories around the desired one can be repetitively
practiced, even at the onset of learning. An increased viscoelastic force at the onset of learning is also effective in learning
schemes that utilize motor command errors read from thefeedback controller as learning signals for internal models(Kawato et al. 1987, 1993).
Active and predictive control of viscoelastic force
Whether the CNS actually regulates the joint viscoelasticity
or whether this parameter is simply an incidental by-product of
the overlapping activity of agonist and antagonist muscles hasremained an unsettled problem (Gomi 1996; Smith 1996; Thoroughman and Shadmehr 1999; Van Galen et al. 1996). The
present results showed that performance error can explain
future rIMCJ, especially in the shoulder, during deceleration;
this supports the idea that the CNS actively and predictively
controls the viscoelasticity. Viscoelasticity regulation by theoff-line feedback of the performance may assist in learning
movements. Such short-term interactions are not yet implemented in other biological models recently proposed, in which
the feedforward compensation of the external force is parallelwith the co-contraction mechanism to stabilize the movement(Flash and Gurevich 1997; Gribble and Ostry 2000; Katayama
et al. 1998; Wang et al. 2001).
Van Galen and his colleagues (Van Galen and DeJong 1995;
Van Galen and Schomaker 1992; Van Galen et al. 1990)
hypothesized that stiffness control is an effective means ofaccuracy control. They stated that a high feedback gain will
reduce the endpoint error. The CNS may increase the viscoelastic forces to execute more accurate movements. This isin good agreement with the present results, i.e., after unsuccessful trials, the rIMCJ increased to result in successful trials.
Increases were more frequent in the latter half of each reaching
movement duration, and this suggests the intention of thesubjects to increase the endpoint accuracy.
The hypothesis that an increase in stiffness improves theendpoint accuracy, however, may appear to contradict theobserved long-term decrease in rIMCJ. Although the rIMCJ
was lower at the end of learning than at the beginning oflearning, the accuracy of the movements increased significantly after extensive learning. This hypothesis of stiffness
control may also be at odds with Harris and Wolpert (1998),
because stiffness increases require larger motor commands,
which in turn introduce higher levels of noise that causes larger
trajectory errors. One explanation is that the CNS adopts dual
strategies to improve the accuracy and chooses which strategy
to use depending on the accuracy of the internal models. When1001

the internal models are immature, the CNS tries to increase theviscoelastic force and utilizes muscles and neural feedback tocompensate for unexpected interaction forces. Once an internal
model is acquired, the CNS tries to decrease the motor command amplitude by decreasing the viscoelastic force, reducing
the noise and increasing the accuracy. An increase in viscoelastic force is effective for a temporary increase in accuracy
against perturbations; any muscle activation should be decreased to maintain a good performance over the long run.
We can think of several advantages in controlling an arm
with low viscoelasticity: low muscle activation increases efficiency and prevents muscle fatigue, and a compliant limb
reduces the potential for injury when it contacts an object in theenvironment. However, viscoelasticity is indispensable for stable control. Without viscoelasticity, spatiotemporal disturbances in the environment, such as objects to be manipulated,
force-fields, or even changes in trajectories, can perturb limbs,
resulting in instability. It therefore makes sense computationally to increase the limb compliance when moving in a stableenvironment, while adaptively increasing the viscoelastic force
in response to instability in the environment (Lacquaniti et al.
1993). This study suggests that the CNS is equipped with aningenious mechanism of learning and controlling movements,
which regulates viscoelasticity on a short- and long-term basis,
depending on performance error, and finally acquires smooth
and accurate movements while maintaining stability during theentire learning process.
APPENDIX1

Advantages of the BMF model
The BMF model has the following advantages.
1) Multivariate mean-nonstationary time-series data can be decomposed into its nonstationary trend components and other mutuallyrelated stationary cyclical components described by an AR modelthrough one single procedure.
2) The AIC can be used to evaluate the accuracy of the decomposition.
3) The structure of the mutual relationships among the stationary
cyclical components (e.g., impulse responses) can be derived fromestimated AR coefficient matrices by applying system analysis techniques (Akaike and Nakagawa 1972; Kato and Kawahara 1998).
4) The reliability of the estimated parameters of the model can be
confirmed by reapplying the system analysis to the simulated data
sets, which are themselves generated based on the estimated model(Monte Carlo simulation).

Evaluation of the model
To examine the reproducibility of the estimated model, we generated data by using Gaussian noise as the system-noise input in themodel estimated from the measured data (the Monte Carlo method),and re-estimated the model parameters by fitting the AR model to thesimulated data. Examining empirical data generated through Monte
Carlo simulation is an effective way of investigating the properties
and reproducibility of such a model, especially when samples arelimited and it is impossible to get another data set to test the model(Lutkepohl 1993). Actually, simulated data can be produced by applying artificially generated Gaussian white noise as input to the noise
components of the model estimated from the measured data. If theassumptions about the whiteness and Gaussian distribution of thenoise are satisfied, the simulated data should reflect the stochastic
characteristics of the measured data. If the assumptions are not ade-

J Neurophysiol • VOL 88 • AUGUST 2002 • www.jn.org

1002

OSU, FRANKLIN, KATO, GOMI, DOMEN, YOSHIOKA, AND KAWATO

FIG. A1. A and B: predictability of (A) rIMCJ and(B) performance error for 1 typical subject (subject
Y.O.) using the model. Plots show normalized data(solid curves) and 1-step predictions (dashed curves)
calculated from the proposed statistical model. C and D:
distribution of prediction errors for (C) rIMCJ and (D)performance error with fitted Gaussian curves for subject Y.O. E and F: autocorrelation function of prediction
error for (E) rIMCJ and (F) performance error with the
95% confidence interval (dash-dot lines) for 6 subjects.

quate, the white Gaussian noise input will not be able to simulate themeasured data. Therefore by reapplying the same system-analysis
procedure to the simulated data, we checked whether the simulated
data had physiological properties identical to the measured data.
Similar estimation results obtained from both the measured and simulated data sets verified the adequacy of the estimators, and at thesame time, the assumption of the white Gaussian noise. In our simulation, Gaussian random numbers were computed using the Marsaglia method (Marsaglia and Bray 1964) and input to each noise
component of the model. The parameters were again estimated byfitting the BMF model to the simulated data. Impulse responses werecomputed from the parameters obtained by the simulation data andcompared with the impulse responses obtained from the real data.
Similar impulse responses and the same AR order obtained from both
real data and the simulated data set confirm the validity of the modeland the estimators (Fig. 9 and Table 2).
Figure A1, A and B, shows the prediction capability of the modelfor the (A) rIMCJ and (B) performance error for subject Y.O. The
solid curves denote normalized data and the dashed curves denote
one-step predictions computed from the proposed statistical model.
The one-step predictions demonstrate the predictive part of the model
described as an auto-regressive form by Eqs. 7 and 8. Differences
between one-step predictions and the data (prediction error) demonstrate the stochastic nature of the system, mainly described by thesystem noise ␧r(n) and ␧e(n) in Eqs. 7 and 8. A similarity of a one-step
prediction with the data demonstrates the predictability of the model(goodness of fit), although they did not perfectly fit one anotherbecause of the stochastic nature of the system.
As mentioned above, in deriving the properties of the estimators,
white Gausssian noise was assumed as the system noise input to themodel. If the assumed model were appropriate, the interactions would
be fully extracted by the model and only the stochastic components
would remain as the prediction error. In this case, if the BMF model
structure and estimated model parameters are appropriate, the prediction error should be white and Gaussian as assumed. On the other
hand, if the assumed model were poor or inappropriate to describe theunderlying system, the interactions that could not be modeled would

appear in the prediction error. We confirmed the validity of theassumption by an autocorrelation function and the distribution ofprediction errors. Figure A1, C and D, shows distributions of prediction error for the same subject. The observed normally distributed
prediction error demonstrated that the Gaussian noise assumption was
satisfied by the data. All of the other subjects showed similar Gaussian
distributions of prediction error. Figure A1, E and F, shows autocorrelations of the prediction error for all subjects. The dash-dot lines
denote 95% confidence intervals. The lack of correlations at lags ⬎0demonstrates that the assumption of the whiteness of the noise wasalso satisfied. The white and Gaussian prediction errors observed forall subjects confirmed the validity of the assumed model and thereliability of the estimated parameters and impulse responses.
APPENDIX 2

The details of dynamic stiffness estimation and IMCJ
computation
For estimation of dynamic stiffness and corresponding dynamic IMCJ,
the subject performed a series of reaching movements while coupled to
the PFM. The right forearm of the subject was firmly attached to the PFM
using a molded plastic cuff. The cuff immobilized the wrist joint, permitting movement only at the shoulder and elbow. The right forearm was
supported on a beam fixed to the handle of the PFM. The subject
performed a reaching movement from a start target (located 31 cm infront of the shoulder) to an end target (located 56 cm in front of theshoulder). The subject was required to produce movements of 600 ⫾ 100
ms duration. The force-fields applied were null force-field, velocityFx
13 ⫺18ẋ
dependent force-field⫽ ␣, and positionFy
18 13
ẏ
Fx
␤x
⫽
. Here ␣ and ␤ were chosen
dependent force-field
Fy02
depending on the subject’s force ( ⁄3 ⱕ ␣ ⱕ 1; 300 ⱕ ␤ ⱕ 500). The
endpoint stiffness of the arm was determined by applying displacement ineach of eight directions during the mid-portion of a movement. The冉冋 册 冋
册冋 册冊
冉冋 册 冋 册冊

J Neurophysiol • VOL 88 • AUGUST 2002 • www.jn.org

CHANGES IN CO-CONTRACTION WITH MOTOR LEARNING

displacement onset occurred 150 ms after the movement onset. Eighty
trials without displacement were randomly intermingled with the displacement trials for a total of 160 trials. The details of the stiffness
estimation and the experimental setup are presented elsewhere (Burdet etal. 2000, 2001). The endpoint stiffness was converted to joint stiffness. At
the same time, EMG signals were recorded from six arm muscles. The
EMG signals of 80 trials without displacement were rectified and averaged for each muscle for a period of 300 ms during the mid-portion of amovement, including the portion where the perturbations were applied in
the displacement trials. The means over 80 no-displacement trials wereused to compute shoulder and elbow IMCJs.
We thank Dr. Y. Tohkura of Nippon Telegraph and Telephone Corporation
for continuing encouragement, and Dr. T. E. Milner of Simon Fraser University for improving the manuscript.
This work was supported by the Human Frontier Science Program, theSpecial Coordination Fund to Brain Science, and the Telecommunications
Advancement Organization of Japan.

REFERENCESAKAIKE H. A new look at the statistical identification. IEEE Trans Automat
Control 19: 716 –723, 1974.
AKAIKE H AND NAKAGAWA T. Statistical Analysis and Control of Dynamics
Systems. Tokyo, Japan: Scientific Publishers, 1972.
BASMAJIAN JV AND DE LUCA CJ. Muscles Alive. Baltimore, MD: Williams andWilkins, 1985.
BENNETT DJ, HOLLERBACH JM, XU Y, AND HUNTER IW. Time-varying stiffnessof human elbow joint during cyclic voluntary movement. Exp Brain Res 88:
433– 442, 1992.
BERNSTEIN N. The Coordination and Regulation of Movements. Oxford: Pergamon Press, 1967.
BIZZI E, ACCORNERO N, CHAPPLE W, AND HOGAN N. Posture control andtrajectory formation during arm movement. J Neurosci 4: 2738 –2744, 1984.
BIZZI E AND MUSSA-IVALDI FA. Neural basis of motor control and its cognitive
implications. Trends Cog Sci 2: 97–102, 1998.
BURDET E, OSU R, FRANKLIN DW, MILNER TE, AND KAWATO M. The central
nervous system stabilizes unstable dynamics by learning optimal impedance. Nature 414: 446 – 449, 2001.
BURDET E, OSU R, FRANKLIN DW, YOSHIOKA T, MILNER TE, AND KAWATO M.
A method for measuring endpoint stiffness during multi-joint arm movements. J Biomech 33: 1705–1709, 2000.
FELDMAN AG. Functional tuning of nervous system with control of movement
or maintenance of a steady posture. 2. Controllable parameters of themuscles. Biophysics 11: 565–578, 1966.
FLANAGAN JR, OSTRY DJ, AND FELDMAN AG. Control of trajectory modifications in target-directed reaching. J Mot Behav 25: 140 –152, 1993.
FLASH T. The control of hand equilibrium trajectories in multi-joint arm
movements. Biol Cybern 57: 257–274, 1987.
FLASH T AND GUREVICH I. Arm trajectory generation and stiffness control
during motor adaptation to external loads. In: Self-Organization, Computational Maps and Motor Control, edited by Morasso PG and Sanguinetti V.
Amsterdam: Elsevier, 1997, p. 423– 482.
FRANKLIN DW, OSU R, BURDET E, KAWATO M, AND MILNER TE. Learning tostabilize unstable dynamics: III. EMG correlates. Soc Neurosci Abst 26: 171,
2000.
GANDOLFO F, MUSSA-IVALDI FA, AND BIZZI E. Motor learning by field approximation. Proc Natl Acad Sci USA 93: 3843–3846, 1996.
GHAHRAMANI Z, WOLPERT DM, AND JORDAN MI. Generalization to local
remappings of the visuomotor coordinate transformation. J Neurosci 16:
7085–7096, 1996.
GOMI H. Is stiffness a byproduct or a target? Behav Brain Sci 19: 450 – 451,
1996.
GOMI H AND KAWATO M. Task dependent stiffness of human multi-joint arm
during point-to-point movement. NTT Basic Res Labs Inform Sci Res Lab,
ISRL-95-104, 1995.
GOMI H AND KAWATO M. Equilibrium-point control hypothesis examined bymeasured arm stiffness during multijoint movement. Science 272: 117–120,
1996.
GOMI H AND KAWATO M. Human arm stiffness and equilibrium-point trajectory during multi-joint movement. Biol Cybern 76: 163–171, 1997.

1003

GOMI H AND OSU R. Task-dependent viscoelasticity of human multijoint arm
and its spatial characteristics for interaction with environments. J Neurosci
18: 8965– 8978, 1998.
GRIBBLE PL AND OSTRY DJ. Compensation for loads during arm movements
using equilibrium-point control. Exp Brain Res 135: 474 – 482, 2000.
GRIBBLE PL, OSTRY DJ, SANGUINETI V, AND LABOISSIERE R. Are complex
control signals required for human arm movement? J Neurophysiol 79:1409 –1424, 1998.
HARRIS CM AND WOLPERT DM. Signal-dependent noise determines motor
planning. Nature 394: 780 –784, 1998.
HOGAN N. An organizing principle for a class of voluntary movements.
J Neurosci 4: 2745–2754, 1984.
HUNTER IW AND KEARNEY RE. Dynamics of human ankle stiffness: variation
with mean ankle torque. J Biomech 15: 747–752, 1982.
IMAMIZU H, MIYAUCHI S, TAMADA T, SASAKI Y, TAKINO R, PUTZ B, YOSHIOKA
T, AND KAWATO M. Human cerebellar activity reflecting an acquired internal
model of a new tool. Nature 403: 192–195, 2000.
IMAMIZU H, UNO Y, AND KAWATO M. Internal representations of the motor
apparatus: implications from generalization in visuomotor learning. J Exp
Psychol Hum Percept Perform 21: 1174 –1198, 1995.
ISHIGURO M AND AKAIKE H. DALL: Davidon’s Algorithm for Log Likelihood
Maximization. Tokyo: The Institute of Statistical Mathematics, 1989.
ISHIGURO M, KATO H, AND AKAIKE H. Ardock: An Auto-Regressive Model
Analyzer. Tokyo: The Institute of Statistical Mathematics, 1999.
KATAYAMA M, INOUE S, AND KAWATO M. A strategy of motor learning using
adjustable parameters for arm movement. Proc 20th Annu Int Conf IEEEEng Med Biol Soc, 1998, p. 2370 –2373.
KATAYAMA M AND KAWATO M. Virtual trajectory and stiffness ellipse duringmulti-joint arm movement predicted by neural inverse models. Biol Cybern
69: 353–362, 1993.
KATO H AND KAWAHARA H. An application of the Bayesian time series modeland statistical system analysis for F0 control. Speech Commun 24: 325–359,
1998.
KAWATO M, FURUKAWA K, AND SUZUKI R. A hierarchical neural-networkmodel for control and learning of voluntary movement. Biol Cybern 57:
169 –185, 1987.
KAWATO M, GOMI H, KATAYAMA M, AND KOIKE Y. Supervised learning for
coordinative motor control. In: Computational Learning and Cognition,
SIAM Frontier Series, edited by Baum EB. Philadelphia, PA: Society for
Industrial and Applied Mathematics, 1993, p. 126 –161.
KITAZAWA S, KIMURA T, AND UKA T. Prism adaptation of reaching movements:
specificity for the velocity of reaching. J Neurosci 17: 1481–1492, 1997.
KOIKE Y AND KAWATO M. Virtual trajectories predicted from surface EMG
signals. Soc Neurosci Abstr 19: 543, 1993.
KOIKE Y AND KAWATO M. Estimation of dynamic joint torques and trajectoryformation from surface electromyography signals using a neural networkmodel. Biol Cybern 73: 291–300, 1995.
KUECHLE DK, NEWMAN SR, ITOI E, MORREY BF, AND AN KN. Shoulder muscle
moment arms during horizontal flexion and elevation. J Shoulder Elbow
Surg 6: 429 – 439, 1997.
LACKNER JR AND DIZIO P. Rapid adaptation to Coriolis force perturbations ofarm trajectory. J Neurophysiol 72: 299 –313, 1994.
LACQUANITI F, CARROZZO M, AND BORGHESE NA. Time-varying mechanicalbehavior of multijointed arm in man. J Neurophysiol 69: 1443–1464, 1993.
LATASH ML AND GOTTLIEB GL. Reconstruction of shifting elbow joint compliant characteristics during fast and slow movements. Neuroscience 43:
697–712, 1991.
LUTKEPOHL H. Introduction to Multiple Time Series Analysis. Berlin: SpringerVerlag, 1993.
MARSAGLIA G AND BRAY TA. A convenient method for generating normal
variables. SIAM Review 6: 260 –264, 1964.
MIALL RC, WEIR DJ, WOLPERT DM, AND STEIN JF. Is the cerebellum a Smith
predictor? J Motor Behavior 25: 203–216, 1993.
MILNER TE AND CLOUTIER C. Compensation for mechanically unstable loadingin voluntary wrist movement. Exp Brain Res 94: 522–532, 1993.
MURRAY WM, DELP SL, AND BUCHANAN TS. Variation of muscle momentarms with elbow and forearm position. J Biomech 28: 513–525, 1995.
MUSSA-IVALDI FA, HOGAN N, AND BIZZI E. Neural, mechanical, and geometricfactors subserving arm posture in humans. J Neurosci 5: 2732–2743, 1985.
NAKANO E, IMAMIZU H, OSU R, UNO Y, GOMI H, YOSHIOKA T, AND KAWATO
M. Quantitative examinations of internal representations for arm trajectory
planning: minimum commanded torque change model. J Neurophysiol 81:
2140 –2155, 1999.

J Neurophysiol • VOL 88 • AUGUST 2002 • www.jn.org

1004

OSU, FRANKLIN, KATO, GOMI, DOMEN, YOSHIOKA, AND KAWATOOSU R AND GOMI H. Multi-joint muscle regulation mechanisms examined bymeasured human-arm stiffness and EMG signals. J Neurophysiol 81: 1458 –
1468, 1999.
POLLICK FE AND ISHIMURA G. The three-dimensional curvature of straightahead movements. J Motor Behavior 28: 271–279, 1996.
SANGER TD. Neural network learning control of robot manipulators using
gradually increasing task difficulty. IEEE Trans Robotics Automat 10:
323–333, 1994.
SHADMEHR R AND MUSSA-IVALDI FA. Adaptive representation of dynamics
during learning of a motor task. J Neurosci 14: 3208 –3224, 1994.
SHIDARA M, KAWANO K, GOMI H, AND KAWATO M. Inverse-dynamics modeleye movement control by Purkinje cells in the cerebellum. Nature 365:
50 –52, 1993.
SMITH A. Does the cerebellum learn strategies for the optimal time-varying
control of joint stiffness? Behav Brain Sci 19: 399 – 410, 1996.
THOROUGHMAN KA AND SHADMEHR R. Electromyographic correlates oflearning an internal model of reaching movements. J Neurosci 19: 8573– 8588,
1999.

THOROUGHMAN KA AND SHADMEHR R. Learning of action through adaptive
combination of motor primitives. Nature 407: 742–747, 2000.
VAN GALEN GP AND DEJONG WP. Fitts’ law as the outcome of a dynamic noise
filtering model of motor control. Hum Mov Sci 14: 539 –571, 1995.
VAN GALEN GP, HENDRIKS AW, AND DEJONG WP. What behavioral benefit
does stiffness control have? An elaboration of Smith’s proposal. Behav
Brain Sci 19: 478 – 479, 1996.
VAN GALEN GP AND SCHOMAKER LRB. Fitts’ law as a low-pass filter effect ofmuscle stiffnesss. Hum Mov Sci 11: 11–21, 1992.
VAN GALEN GP, VAN DOORN RR, AND SCHOMAKER LR. Effects of motor
programming on the power spectral density function of finger and wrist
movements. J Exp Psychol Hum Percept Perform 16: 755–765, 1990.
WANG T, DORDEVIC GS, AND SHADMEHR R. Learning the dynamics of reaching
movements results in the modification of arm impedance and long-latency
perturbation responses. Biol Cybern 85: 437– 448, 2001.
WINTERS JM. Hill-based muscle models: a systems engineering perspective. In:
Multiple Muscle Systems, edited by Winters JM and Woo SL-Y. New York:
Springer-Verlag, 1990, p. 69 –93.

J Neurophysiol • VOL 88 • AUGUST 2002 • www.jn.org