 |
||||
|
||||
|
||||
|
||||
Previous Article | Next Article
The Journal of Neuroscience, October 1, 1999, 19(19):8573-8588
Electromyographic Correlates of Learning an Internal Model of Reaching Movements
Kurt A. Thoroughman and Reza Shadmehr Department of Biomedical Engineering, Johns Hopkins School of Medicine, Baltimore, Maryland 21205-2195
ABSTRACT
TOP
ABSTRACT
INTRODUCTION
Theoretical and psychophysical studies have suggested that humans learn to make reaching movements in novel dynamic environments by building specific internal models (IMs). Here we have found electromyographic correlates of internal model formation. We recorded EMG from four muscles as subjects learned to move a manipulandum that created systematic forces (a "force field"). We also simulated a biomechanical controller, which generated movements based on an adaptive IM of the inverse dynamics of the human arm and the manipulandum. The simulation defined two metrics of muscle activation. The first metric measured the component of the EMG of each muscle that counteracted the force field. We found that early in training, the field-appropriate EMG was driven by an error feedback signal. As subjects practiced, the peak of the field-appropriate EMG shifted temporally to earlier in the movement, becoming a feedforward command. The gradual temporal shift suggests that the CNS may use the delayed error-feedback response, which was likely to have been generated through spinal reflex circuits, as a template to learn a predictive feedforward response. The second metric quantified formation of the IM through changes in the directional bias of each muscle's spatial EMG function, i.e., EMG as a function of movement direction. As subjects practiced, co-activation decreased, and the directional bias of each muscle's EMG function gradually rotated by an amount that was specific to the field being learned. This demonstrates that formation of an IM can be represented through rotations in the spatial tuning of muscle EMG functions. Combined with other recent work linking spatial tunings of EMG and motor cortical cells, these results suggest that rotations in motor cortical tuning functions could underlie representation of internal models in the CNS.
Key words: motor learning; motor control; electromyography; internal model; computational modeling; human
INTRODUCTION
People learn to move novel objects along desired trajectories, in any direction, by simply practicing the task a few times. This adaptation is remarkable because of the computational complexity inherent to learning dynamics (Atkeson, 1989
). Adaptation of a neural internal model (IM), transforming desired trajectories of the hand into appropriate muscle activations, likely underlies this ability (Jordan, 1995
Previous work has reported EMG changes during learning of single-joint movements with different loads (Corcos et al., 1993
We measured EMG as subjects learned to reach while grasping a manipulandum, which created novel forces. We used a biomechanical model to transform EMG into a composite trace, corresponding to activation that specifically compensated for the imposed field. We found that with training, activation of the appropriate musculature gradually shifted from a delayed error-feedback response to a predictive feedforward response. To quantify how these changes corresponded to the formation of an IM, we analyzed the directional tuning functions of movement-initiating EMG. These spatial functions relate EMG to movement direction in hand-centered coordinates (Flanders and Soechting, 1990
Portions of this work have been presented previously in abstract form (Thoroughman and Shadmehr, 1998
MATERIALS AND METHODS
Experimental apparatus. Twenty-four right-handed subjects (12 women and 12 men), students in the Johns Hopkins School of Medicine, provided written consent to participate in this study. With their right hand, subjects grasped the end effector (handle) of a 2 df manipulandum mounted in the horizontal plane (Shadmehr and Brashers-Krug, 1997
Experimental protocol. Subjects made 10 cm point-to-point movements toward 8 mm square targets represented by boxes on the computer monitor. Subjects made movements in four directions (0, 45, 90, and 135°) away from the center of the work space and four directions (180, 225, 270, and 315°) back to the center; the order of the outward directions was determined pseudorandomly. If subjects completed the movement in 500 ± 50 msec, the box "exploded," and the computer generated a pleasing sound. If the cursor reached the box too slowly (in
550 msec), the target filled in blue; if the cursor reached the box too quickly (in
450 msec), the target filled in red. The only instruction provided was to explode as many targets as possible. We provided no instructions regarding straightness of movements or smoothness of trajectories.
Subjects made reaching movements in two or three different dynamic environments. One, termed the "null field," is the condition in which the torque motors do not create any forces; in the null field the subjects encounter only the inertial dynamics of the manipulandum. In the second and third environments, B1 and B2, torque motors produce an additional force as described by the equation:
where
(1) is the force produced at the end effector by the robot's motors,
is the instantaneous velocity vector of the hand, and B is a viscosity matrix. In the first environment, the null field, B0 equals [0 0; 0 0] N · sec · m
1. In the second environment, B1 equals [0 13;
13 0] N · sec · m
Subjects' training was divided into sets of 192 movements. Each set was followed by a 3 min rest period. On a first day of training, all 24 subjects completed three sets of movements in the null field. On a second day, all subjects completed one set in the null field and then three sets in the force field B1. Here we report data only from the second day.
After completion of the training in B1, three groups of eight subjects each completed two sets of movements. One group completed these two sets in the null field after the standard 3 min rest. A second group trained in B2 3 min after B1, and a third group trained in B2 6 hr after B1. During the 6 hr period, EMG electrodes were removed, and subjects left the laboratory. After their return, electrodes were reapplied on approximately the same arm positions as the earlier session, as marked by the earlier abrasion of the skin.
Displacement analysis. To facilitate averaging across movements, all movement data were temporally aligned on the basis of the movement onset, which was defined as when the tangential speed first crossed a threshold (0.03 m/sec). Perpendicular velocity was defined as the component of the velocity vector that pointed perpendicular to the line connecting the initial hand position and target position. Perpendicular displacement was computed by summing perpendicular velocity over the first 250 msec of the movement. In each direction we subtracted away the average perpendicular displacement generated in the null field (first set, day 2).
EMG normalization. Each subject's (bandpassed, RMS'd) EMG traces were normalized based on movement-initiating EMG recorded during the initial null field set. For each muscle m and in each direction of movement, we calculated a scalar activation am by averaging the EMG trace Am(t) over time interval between 50 msec before and 100 msec after the onset of movement (capital letters refer to full time series; small letters refer to scalar measures) and then averaging over all movements made in that direction during the initial null field set. We constructed an 8 × 1 vector
m for each muscle, each element representing a single direction. We then adjusted the scale of each muscle's EMG trace in each movement, producing the normalized trace An, using the equation:
where k is an index of movement number, and m is an index of muscle. Within each subject, therefore, the EMG values recorded during all movements (in the null field and in force fields) were normalized using the same linear transformation. This transformation facilitated comparisons of activation across subjects but preserved the training-induced changes of every subject's EMG.
(2) Polar analysis. We used a form of data analysis termed polar analysis to evaluate changes in EMG. We began the analysis by computing the movement-initiating activation for each subject, binned within each direction across eight movements. We then multiplied each scalar by the unit vector pointing in the direction of movement, so that these vectors, when plotted together, display am as radii pointing in the direction of movement corresponding to that activity. These polar plots summarized the function mapping target direction into initial EMG activity.
To determine the directional bias of each activation function, in each movement bin we added the eight am vectors together to form one resultant vector. This resultant vector had the same direction and is proportional in magnitude to the mean of the eight vectors (if a unit mass would have been placed on the end point of each vector, the mean of the vectors would have been the center of mass of the collection of points). By calculating the orientation of this resultant vector, we detected changes in the angular dependence of am as subjects learned dynamic environments.
Computational modeling: structure of the model. The purpose of the computational modeling was to predict the change in the pattern of joint torques that should result if an adaptive control system learned to completely compensate for the dynamics of the force field. To calculate the changes in joint torques attributable to adaptation, we simulated the dynamics of the four-link system of Figure 1 through the following coupled differential equations:
(3)
where I and G were inertial and coriolis/centripetal matrix functions, M was the torque field produced by the robot's motors, i.e., the environment,
(4) was the adaptive controller implemented by the motor system of the subject,
*(t) was the reference trajectory planned by the motor control system of the subject, Jr was the Jacobian matrix describing the differential transformation of coordinates from end point to joints,
were column vectors representing joint positions of the subjects and the robot (Figure 1), and the subscripts s and r denoted subject or robot matrices of parameters.
View larger version (25K):[in this window]
[in a new window]
Figure 1. Experimental apparatus. Subjects make reaching movements while grasping a 2 df manipulandum (robotic arm). The manipulandum could be programmed to produce a force field.
The robot's passive mechanical dynamics, Ir and Gr, as well as a typical subject's passive dynamics of the arm, Is and Gs, were previously estimated (Shadmehr and Brashers-Krug, 1997
K was a linear estimate of the subject's joint stiffness and was set to a constant K = [15, 6; 6, 16] N · m
(5) s, and model of the force field environment,
(
*). We assumed the desired trajectory,
We simulated Equations 3 and 4 for a configuration of the arm in which the hand made the first movement starting from x =
. This field can be written in joint coordinates of the robot or the subject, i.e., the M variable in Equation 3, using simple kinematic transformations (Shadmehr and Mussa-Ivaldi, 1994
The torques produced by the controller of Equation 3 on the shoulder and elbow joints were
s(t,
) and
(6)
where
[Fb, Ft, Fa, Fp]T is a vector of muscle forces, JmT is the transpose of the differential transformation from muscle lengths to joint angles, i.e., d
/d
There is no general consensus regarding the magnitude of muscle moment arms (i.e., matrix Jm in Eq. 6) for the human arm in the horizontal plane. Here we considered six different models (Table 1). Four of the matrices (models A-D) were derived presuming different levels of effective biarticulation of biceps and triceps: model A presumed no biarticulation, and models B-D presumed that the moment arms of biceps and triceps across the shoulder were one-half, the same, and twice the moment arms of the anterior and posterior deltoids, respectively. The fifth model (E) was derived from cadaver data (Wood et al., 1989
View this table: [in this window]
[in a new window]
Table 1. Computational representation of an internal model in terms of expected changes in the angular orientation of the EMG, as represented by a spatial function in hand centered coordinates, for each muscle Computational modeling: predicting rotation of directional bias. For each moment arm model, the simulation produced a time series of muscle forces appropriate to make movements in the null field. In each direction of movement i, we averaged over the first 150 msec of predicted force production (corresponding to the interval 50 msec before through 100 msec after movement initiation) to determine for each muscle a scalar variable fm(i). We performed polar analysis (see above) on these scalars to determine a resultant vector to summarize, in null field movements, the angular dependence of fm on movement direction. We then performed the same analysis on the forces appropriate to move in the force field B1. From these two polar analyses we determined, for each muscle, the difference between the orientations of fm resultant vectors in the null field and in force field B1.
The results for each model are shown in Table 1. We found that changing moment arm matrices tended to change the orientation of preferred directions in both null field and B1 but did not tend to change the angle between them. All models predicted approximately the same levels of rotation in the tuning functions of biceps and triceps. For the anterior and posterior deltoids, predicted rotations were very similar for all models except that of model D. The results in Table 1 show how much rotation of the spatial activation function, as represented in hand-centered coordinates, should occur as a subject learns an internal model for field B1.
Computational modeling: predicting time series of appropriate activations. Studies that have recorded EMG during single-joint movements have reported increases or decreases in the activation of individual muscles during learning. This analysis is appropriate for tasks with a single degree of freedom, but here we needed to consider movements with 2 df, actuated by multiple muscles, moving in several directions. To facilitate analysis of EMG across directions of movement, we derived an inverse muscle model, whose parameters were determined by subjects' EMG in the null field. The model was validated by predicting EMG patterns in the force field. These predictions were then used to remap EMG into composite traces based on their functional relevance.
To determine the muscle activations appropriate for our task, we first used Equations 3-6 and muscle model C to calculate the time series of muscle forces necessary to make a minimum jerk reaching movement, Fi,m,n(t), where i indicates movement direction, m muscle, and n force field (n = 0 signifies the null field, n = 1 the force field B1). Next, we transformed the modeled muscle force Fi,m,n(t) into appropriate muscle activations. To estimate the relationship between predicted muscle force and recorded EMG [Ai,m,n(t)], we found a piecewise linear fit between the force predictions in the null field and the recorded EMG from our subjects, using the following analyses. In the null field, for each muscle m and in each direction, we calculated the average predicted force and recorded EMG during the agonist period, from 50 msec before through 100 msec after the onset of the movement (fmag and amag), and during the antagonist period, from 200 msec through 300 msec into the movement (fmant and amant). We then fit the null field forces across all eight directions to the recorded null field EMG using the equations:
(7)
where
are 8 × 1 vectors reflecting recorded EMG and predicted force in all eight directions, and p are parameters of the fit.
We calculated the parameters to find the least-squared fit of the null field data and found that the fit was significant for all four muscles, during both the agonist and antagonist bursts (for all eight fits, r > 0.88; p < 1 × 10
View larger version (34K):[in this window]
[in a new window]
Figure 2. Data used for building a simple force activation model for the anterior deltoid and data used for testing the validity of the model. We fitted the predicted muscle forces (as computed by an inertial model of the arm and best estimates of moment arms) to subjects' mean EMG while moving in the null field during the agonist period (
View this table: [in this window]
[in a new window]
Table 2. Values of fitted parameters p in Equation 7 For each movement direction, each muscle, and both dynamic environments, we used the parameters p to map predicted force (f) into predicted EMG (e; note that a refers to actual activations) at baseline attributable to the reciprocal activation (base), during movement initiation (ag, from
(8)
For each movement direction i and force field n, x(i, m, n) for two of the muscles were positive; the size of these activations predicted the proportion that each of those muscles participated as an agonist in making a movement in the null field (for n = 0) or as an agonist for resisting the forces of field B1 (for n = 1). On the other hand, x(i, m, n) was negative for the remaining two muscles. These negative values predicted that the two remaining muscles produced a force opposite that of the imposed field and should be inhibited. We constructed 1 × 4 vector operators v,
(9)
to identify functional agonists (+) and antagonists (
View this table: [in this window]
[in a new window]
Table 3. Values of the vector operator v, as defined in Equation 9 Each vector operator predicts, in the four-dimensional space of muscle activation, in what direction activation should increase or decrease (at the peak of activity) to effectively make movements in a given field. To give an example, consider a movement toward a target at 90°. In the null field (n = 0), the vector operator v90,0,+ predicts that to initiate an accurate movement, the triceps and anterior deltoid should be activated above their baseline levels. Furthermore, if one considered the triceps and anterior deltoid activations to form an (x, y) plane, the model predicts that the average activity over the agonist period would point in the direction of the unit vector (0.859, 0.513). To move in the field B1 toward 90°, v90,1,+ predicts that biceps and anterior deltoid should be activated more strongly than in the null field. Midway through the movement (200-300 msec into the movement), the unit vector pointing in the direction of the most appropriate increase in activation would equal (0.869, 0.495).
Using model predictions to remap subjects' EMG. We used the vector operators (Table 3) to remap EMG recorded from the four muscles into composite traces, T:
where
(10) i(t)
The transformation in Equation 10 remapped the full time series of EMG traces using the predicted optimal activation over specific time intervals (the average over the agonist period for null field data, over midmovement for the increase appropriate for field B1). These remappings retained the natural temporal structure of the EMG traces for each muscle; the activation in each trace at each time point was weighed equally. When the remapping of the increase appropriate for B1 was instead determined by predictions during other time intervals, the direction of appropriate activation changed very little over the range from 50 msec before through 300 msec after the onset of movement (averaged across directions of movement, change in unit vector direction
The composite trace Ti,1,+ summarized EMG activity that we estimated as effectively countering the forces of field B1. This trace was averaged across directions and binned over 64 movements (three bins per target set of 192 movements). The averaged and binned traces were then fit with fourth-order polynomials. Of primary interest was the question of how this function evolved as subjects practiced in the force field. A major component of this change was an apparent shift in the peak of the function. Bootstrapping techniques (see below) were used to estimate 95% confidence intervals of timing of the peak.
Wasted contraction: a measure of co-contraction. Previous studies have measured "co-contraction" through the addition of the EMG traces of opposing muscles (Milner and Cloutier, 1993
For each subject, we first remapped the (previously scaled) EMG of each movement into composite traces (Eq. 10). At each sampling point, we then determined, in each agonist-antagonist pair of composite traces, which trace was smaller. Because this activation was cancelled by the larger, opposing activation, the time series of this smaller activity was termed wasted contraction. This amount was then subtracted from the larger trace, creating a time series of "effective contraction." The wasted contractions were averaged across force directions (n
0, 1 from above), averaged across movement direction, and binned over 32 movements (six bins per set). Each subject's wasted contraction was then scaled by the maximum effective contraction measured over bins in the null field and in B1. The wasted contraction traces were then averaged across subjects. This average wasted contraction was fit with a fourth-order polynomial; the maximum of this polynomial was calculated and compared during different stages of training. The 95% confidence intervals for the change in maximum wasted contraction were calculated by bootstrapping (see below).
Calculation of force produced by the human controller. We also used the structure and parameters of our model of the human-robot interaction to estimate the amount of force subjects create during movements. Using Equations 3 and 4 and using our estimates of the human and robot dynamical parameters, we transformed individual subjects' hand position and velocity into estimated torque produced by the subjects' muscles (the acceleration terms in these equations were estimated by filtering the velocity data through a second-order Savitsky-Golay filter). We then multiplied this torque by Js
Bootstrapping. Confidence intervals for the rotation of EMG resultant vectors, the timing of the first peak of the B1-appropriate EMG trace, and the maximum level of wasted contraction were determined by using bootstrapping (Fisher, 1993
RESULTS
EMG correlates of learning: a single direction of movement
Let us initially consider movements made toward a single direction. When subjects made movements in the null field at 135°, their EMG traces (Fig. 3, gray lines) revealed that the anterior deltoid was active early to initiate the movement, whereas triceps and posterior deltoid were active later to brake the movement.
View larger version (16K):[in this window]
[in a new window]
Figure 3. Hand paths and EMG recorded during movements. Left, Hand paths of a typical subject in her first and last (72nd) movements toward a target at 135° in the force field B1. The first movement is significantly perturbed from the straight line trajectory. With practice, movements become straight again. Dots are at 10 msec intervals. The arrow indicates 250 msec into the movement. This is the position at which we measured perpendicular displacement from a straight line. Remaining plots, Processed EMG (normalized units) in the null field (all 24 movements toward 135°, gray line), during initial stages of training in the force field (first 8 movements toward 135°, thin black line), and late in training in the force field (last 8 movements, thick black line). Movement begins at t = 0. Each EMG trace represents data averaged across all 24 subjects. The vertical dotted lines delimit the 150 msec interval over which EMG is averaged to calculate the scalar variable am, representing time-averaged agonist burst EMG.
Early in training in force field B1, when subjects attempted to make movements toward 135° they generated inaccurate hand paths that diverged from the target direction (Fig. 3, top trajectory). Midway through a typical movement, the activity in biceps (Fig. 3, thin black line) increased; the hand path was then corrected toward the target. The timing of the increased biceps activation suggested that this muscle was activated in response to the robot forces that extended the elbow during the reach. The increased activation of biceps was possibly an error-feedback response mediated via a stretch reflex mechanism. In contrast to these field-induced changes in the biceps, the time series of activations for anterior and posterior deltoids in the force field closely resembled the activations recorded in these muscles in the null field.
As subjects trained in field B1, the movements became straighter (Fig. 3). The EMG from these improved movements featured larger biceps activity very early, preceding the movement itself (Fig. 3, thick black line). The biceps then became less active later in the movement, during the same period in which it was active when subjects first trained in B1. Therefore, a major component of the adaptation process was a modification of the EMG patterns necessary to initiate the movement.
To summarize the changes in early EMG in movements toward 135°, we averaged EMG over the time interval from 50 msec before to 100 msec after the onset of movement. This scalar measure of activation was termed am. Figure 4 shows the evolution of am as subjects moved in the null field and then trained in B1. The figure also displays perpendicular displacement, a measure of movement error that is computed 250 msec into the movement. These data show that as subjects trained in B1, the curvature of the movements became significantly smaller (p < 5 × 10
View larger version (17K):[in this window]
[in a new window]
Figure 4. Perpendicular displacement (top) and subjects' am (time-averaged agonist burst EMG) during movements toward 135°. Subjects completed one set of movements in the well learned null field environment and then three sets in the force field B1. Perpendicular displacement is measured 250 msec into the movement, averaged across directions and subjects; these data reflect change from average perpendicular displacement in null field movements. x-axis tick marks indicate breaks between sets of movements. Each point includes data from eight movements (3 data points per set). Error bars reflect 95% confidence intervals of the mean.
In summary, the changes in EMG traces measured in movements toward 135° suggested that when subjects were first exposed to B1, they generated additional muscle activations late in the movement to correct for movement curvature. With further training, subjects learned to generate this added activation early in the movement, predicting rather than responding to the force field.
EMG correlates of learning across all directions of movement
To quantify the correlation between changes in EMG and the formation of an internal model, we wanted to extend our analysis from a single direction to many directions of movement. This extension is hampered by the fact that each muscle makes different contributions to movements in different directions. We circumvented this obstacle by using the predictions of our simulated controller. We considered the time pattern of EMG in the muscles for a given movement as a path in a four-dimensional space, each axis representing a muscle. Our model also predicted paths in this space: one path for moving in a null field and one for moving in a force field. The difference between these two predicted paths indicated activity in that space which should be increased, peaking midway through the movement, to move accurately in B1. In each movement direction, the increased activity projected into a single component of muscle activity (Eq. 9, Table 3; see Materials and Methods). The discovery of appropriate activations allowed us to project the EMG recorded from each subject, during each movement in the field B1, into functionally appropriate composite traces (Eq. 10).
We used this approach to map the four EMG traces into four composite traces. One of the resulting traces quantified the component of activations that was most appropriate for moving in the null field; another one quantified the component that most appropriately resisted the additional forces created by field B1. These were termed the null field and B1-appropriate EMG traces. This transformation remapped the EMG from directionally dependent individual muscles to directionally independent composite traces, permitting the comparison and averaging of these traces across directions.
We calculated the null field- and B1-appropriate EMG traces during the first null field set and the third set in field B1, respectively. The resulting traces were averaged over all movement directions and all subjects to arrive at a single trace. The results are shown in Figure 5 (solid lines). The left plot (Fig. 5A) shows the average null field-appropriate EMG trace produced by subjects as they reached in the null field. The peak of this trace was at
View larger version (15K):[in this window]
[in a new window]
Figure 5. Subjects' composite EMG traces (solid lines) and forces (dotted lines). A, The solid line is the null field-appropriate EMG averaged over all subjects and all movements in the initial null field set (192 movements). The dotted line is the average force in the direction of the target actually produced during the same movements. B, The solid line is the B1-appropriate EMG trace, averaged over all subjects and over all movements during the final set of training in force field B1. The dotted line is the average force produced perpendicular to the direction of target during the same movements. Forces in both plots were calculated from subject trajectories, using the model structure outlined in Equations 3 and 4. For the time series of EMG and force, maximum range of 95% confidence intervals of EMG and force were as follows: null field-appropriate EMG, ±5.3 units; parallel force, ±0.28 N; B1-appropriate EMG, ±6.2 units; perpendicular force, ±0.13 N.
The composite EMG trace shown in Figure 5A represents the component of activation appropriate for moving the hand toward a target, i.e., the agonist burst of activity. A second composite EMG trace (data not shown) produces force in the direction opposite the movement direction. This composite trace (defined by v0,
This close correspondence was maintained when we compared the B1-appropriate EMG with the pattern of forces that subjects had produced perpendicular to the direction of target during the training. Fig. 5B shows the B1-appropriate EMG and the corresponding forces, both averaged across subjects during the last set of training in the force field. With training, subjects had learned to alter their EMG early into the movement [Fig. 5B, peak of the solid line at 78 msec; 95% CIM, (67, 91) msec] so that the muscles produced a bell-shaped force pattern perpendicular to the direction of the target, effectively canceling the force field.
The relationship between the composite EMG traces and the subject-generated forces confirms the aptness of the transformation of EMG into composite traces. The vector operators in Table 3, created from the predicted activations of the biomechanical model, have elicited the components of EMG appropriate for moving in the null field and resisting the additional force produced by field B1. We will now examine the component appropriate for field B1 throughout the training in the force field, to determine how the neural activation of the appropriate musculature evolves with learning.
Within subjects, we averaged the B1-appropriate EMG across directions and across bins of 64 movements (three bins per set, nine bins spanning training in B1). The mean trace for each subject for each bin was averaged across subjects. The resulting traces are shown in Figure 6. Because the imposed force field was a function of hand velocity, and hand velocity was generally unimodal, we expected the EMG that effectively countered this force to be also unimodal. Indeed, we found that the B1-appropriate EMG was always unimodal. Furthermore, because the response to the imposed force field could initially be only through a delayed error feedback system (e.g., through spinal reflexes), we expected the initial B1-appropriate EMG to lag the imposed forces. This was also observed.
View larger version (19K):[in this window]
[in a new window]
Figure 6. The composite EMG trace appropriate for field B1 shifts with training. Top, Binned and averaged increase in the B1-appropriate EMG trace. Each trace is averaged over 64 movements and over all 24 subjects. The progression of plots from top to bottom shows EMG traces from early to late in training. The dashed lines in the top and bottom axes represent the magnitude of the force created by the viscous field B1, which remained consistent throughout training in B1. Although early in training the peak of the field-appropriate EMG lagged the imposed force, with training it preceded the imposed force. The timing of the peak was estimated by a fit to a fourth-order polynomial. The best estimate of the peak and its change with training is marked by the line that crosses the EMG traces from top to bottom. In the bottom figure, the timing of the peak of the EMG traces is plotted versus movement number. The dotted lines connecting data points represent the 3 min breaks between sets. Error bars reflect 95% confidence intervals of the mean (across all 24 subjects), as determined by bootstrapping. The dotted horizontal line represents the timing of the peak of the force created by the force field.
With practice in the force field, the peak of the B1-appropriate EMG trace shifted to occur earlier into the movement. To estimate the peak location of each EMG trace, the data were fitted to a fourth-order polynomial. We found that when subjects started to train in B1, the peak of the B1-appropriate EMG occurred 246 msec into the movement [95% CIM, (215, 274) msec]. This peak arrived 48 msec after the peak of the force created by the force field B1; the peak of the robot-imposed force coincided with the peak tangential velocity of the hand and on average occurred at 195 msec into the movement throughout training in B1. By the end of the third set of training in B1, subjects had modified their motor commands such that the peak of the B1-appropriate EMG occurred just 67 msec into the movement [95% CIM, (48, 86) msec; in t test of time shift, p < 1 × 10
The changes in timing of the B1-appropriate EMG activity revealed that subjects adapted their neuromotor output throughout all three sets of movements (Figure 6, bottom). This measure was much more sensitive than the kinematic measure of perpendicular displacement, which showed no learning during the second or third sets (compare Figure 8).
With training, subjects generated B1-appropriate EMG earlier in the movement. During the 3 min break between sets, however, the subjects partially unlearned this adaptation and began the next set with activation of B1-appropriate muscles shifted later than at the end of the previous set (Figure 6, bottom). In between the second and third sets of training in B1, this shift forward in time, i.e., the "forgetting," was significant (mean shift, 31 msec; p < 0.05). The shift of the B1-appropriate activation backward in time during the second and third sets and forward in time during breaks between sets indicate that the timing of this component of activation is a valid metric of learning that reveals more about the state of the adaptation than do kinematic measures alone.
Rotation of angular dependence of forces: simulation results
These results suggested that subjects initially responded to B1 by generating appropriate EMG quite late in the movement, perhaps because of an error-feedback mechanism, but with training they learned to activate the same muscles earlier in the movement. This transition suggested that reactive responses early in training effected proactive behavior late in training. We then focused our attention on the subjects' initial, purely feedforward, responses to the task: how do humans generate movement-initiating muscle activations based on the visual cue of a target, a desired end point of a movement? How does this neural computation, mapping stimulus to motor neuronal activation, change with learning?
To guide our investigation of feedforward muscle activation, we revisited the predicted forces generated by our ideal controller. We summarized each muscle's force production in both the null field and in B1 by averaging over the first 150 msec of the time series of force to create a scalar fm. We used polar analysis (see Materials and Methods) to visualize the angular dependence of fm (Figure 7, top). We then calculated resultant vectors (radial vectors in Figure 7, top), which indicate the directions of maximal force production and summarize the dependence of force on the direction of movement. Because we predicted muscle activations using a linear transformation of muscle force, the same polar analysis reveals the predicted directional bias of the muscle activations.
View larger version (28K):[in this window]
[in a new window]
Figure 7. Rotations in muscle-tuning curves as predicted from a computational model that assumes adaptation of an IM and tuning curves as recorded from subjects' EMG. Top, Movement-initiating muscle forces as predicted by the model for movements in the null field (gray) and in B1 (black). Polar plots display predicted movement-initiating force (fm) for each direction of movement. The radii of the scale circles (centered at the origin of each plot) are 20 N. The thick vector represents a tuning curve's resultant vector (preferred direction). The model predicts that the resultant vector should rotate between the null field and B1. The bars below each plot reveal the rescaling of the data (using Eq. 7 and the parameters in Table 2) into predicted muscle activations, in normalized units (nu) of EMG. Bottom, Tuning functions representing subjects' movement initiating EMG (am) during training in null field (gray lines), early force field (thin black lines), and late force field (thick black lines). Thick lines are resultant vectors (preferred directions). The error bars around individual data points reflect 95% confidence intervals of the mean activation.
The simulation results revealed that by learning to compensate for B1, the resultant vector of each muscle should rotate by a specific amount (Fig. 7, top). We found that the degree of rotation was similar over a wide range of moment arm models (Table 1) and remained constant over all levels of co-contraction. (For model C, the directions of the resultant vectors were biceps,
-motor neuronal output, and that when subjects learn to move in B1, the activation function rotates to produce appropriate muscle forces to accurately initiate movements.
Rotation of angular dependence of EMG with training in B1
Armed with the result that the simulated controller's muscle activation resultant vectors rotated between the null field and B1, we analyzed the EMG generated at the beginning of movements in both the well-learned null field and in B1. Just as in our analysis of movements made toward 135°, we averaged each EMG trace over the period from 50 msec before through 100 msec after the beginning of the movement to create the scalar am. We then averaged these scalars over bins of eight movements within directions and then averaged across all subjects. We compared am generated at three stages of training: at the end of the set in the null field, at the beginning of the first set of training in B1, and at the end of the third set of training in B1. The results are shown in Fig. 7, bottom row.
In the null field (Fig. 7, gray lines), each muscle was strongly activated in certain "preferred" directions and much less activated in the opposing directions. This directional bias to the activation function was summarized by using polar analysis to construct a resultant vector for each muscle. The resultant vectors pointed toward the direction (for biceps,
During the first 64 movements in field B1 (Fig. 7, thin black lines), subjects began to generate movement-initiating muscle activations that differed from those generated in the null field. The resultant vectors of the am rotated clockwise from the directions recorded in the null field. Through the end of training in field B1, the resultant vectors summarizing am continued to rotate in a clockwise direction (Fig. 7, bottom). To quantify this adaptation across subjects, we began by determining the orientation of individual subjects' resultant vectors during null field and B1 training. We rotated each subject's activation functions so that the mean null field resultant vectors pointed toward 0°. We then calculated the mean and 95% confidence intervals of the orientations of the am vectors during training in the null field and B1. The results are shown in Figure 8.
View larger version (16K):[in this window]
[in a new window]
Figure 8. Left, Perpendicular displacement of the hand 250 msec into the movement, averaged across movement directions, during movements in null field and force field B1. Remaining plots, Rotation of am resultant vectors during movements in force field B1 with respect to null field. Each point contains data from 64 movements. Each point is mean and 95% confidence interval (across subjects). Dotted lines represent 3 min breaks in training.
As subjects made reaching movements in the null field, the orientation of the EMG resultant vectors remained constant, and the curvature of movements remained small. In the first bin of 64 movements in
