Download Neurophysiology of Motor Cortex: Directional Coding and Muscle Activation and more Study notes Kinematics in PDF only on Docsity!
In Progress in Motor Control III, chap 6, pp 125-166, Latash and Levin (eds), Human Kinetics (2003)
On the role of primary motor cortex in arm movement control
Emanuel Todorov
Department of Cognitive Science
University of California, San Diego
September 11, 2002
Abstract
The role of primary motor cortex (M1) in arm movement control has been the subject of lasting and often heated debates. Early views of low-level muscle control have been discredited by numerous correlations with higher-level parameters related to hand kinematics. The abstract directional coding hypothesis advanced to explain these observations has also been discredited, this time by numerous correlations with lower-level muscle- or joint-related parameters. Taken at face value, the contradictions accumulated over the years may seem to suggest that no coherent account of M1 function can ever be given, and one has to accept a catalog of disparate empirical findings in place of an explanation. This is in principle possible, but, at the risk of disagreeing with some of the neurophysiologists whose findings we find it necessary to re-interpret, we prefer the following more optimistic approach: assume that many of the observed phenomena are in fact epiphenomena, and search for simple underlying mechanisms that give rise to the seemingly contradicting results. To uncover such mechanisms, we need models. The present chapter summarizes a series of recent studies as well as unpublished results aimed at building a coherent model of M1. A coherent model requires an answer to the following question: if the M1 population output encodes just one thing, what could that thing be? The model we arrive at is related to earlier views of muscle-based control. It accounts for later results by incorporating the properties of the musculo-skeletal system: M1 neurons are assumed to activate groups of muscles while taking into account the state-dependence of muscle force production and effects of multijoint mechanics. The M signal needed to compensate for these peripheral properties turns out to be correlated with behavioral parameters on multiple levels of description, quantitatively explaining many results which appear unrelated to the pattern of muscle activation. The chapter ends with an extended discussion of several important issues (origins of directional tuning, muscle synergies, redundancy of population codes, spinal cord function, and equilibrium-point control) from the point of view of our model.
1. Introduction
Primary motor cortex (M1) plays a fundamental role in the control of voluntary arm movements, as evidenced by the profound deficits following M1 lesions. Precisely what that role is has been the subject of lasting debates. These debates are fueled by the numerous correlations found between the activity of M1 neurons and various behavioral parameters. The earliest experiments in awake behaving monkeys (Evarts, 1968) showed that M1 firing is better correlated with isometric force than limb position. Given the dense projection from M1 to the spinal cord, sometimes directly to motoneurons, it is natural to interpret these early observations as evidence that M1 controls muscle force or activation (Evarts, 1981). It was later observed, however, that in multijoint movement tasks the majority of M1 neurons encode not the acceleration (which is proportional to force), but the direction (Georgopoulos et al., 1982) or velocity (Schwartz, 1994) of hand movement. In addition the M1 directional tuning curves were quite broad, which, in the context of single- joint studies that dominated the field at the time, seemed inconsistent with low-level muscle control. These observations, and the ability of the population vector (PV) to recover movement direction, led to the opposing view that M1 sends a much more abstract signal specifying a three-dimensional vector in extrapersonal space. To distinguish between the two alternatives, movement velocity and external load were varied separately in the same experiment (Kalaska et al., 1989). The results however supported both views: velocity- and load-related signals were both present and combined additively in the activity of the same neurons. Although the directional coding hypothesis can be extended to incorporate the observed tuning for load (Kalaska et al., 1989;Georgopoulos et al., 1992;Taira et al., 1996) as well as position (Georgopoulos et al., 1984;Kettner et al., 1988), it offers no explanation as to why all these parameters should be encoded simultaneously. Furthermore, many additional correlations have been observed (see (Todorov, 2000a) that do not seem to fit in either view of M1 function. One can use regression analysis to identify the multiple parameters correlated with the activity of each neuron (Ashe and Georgopoulos, 1994;Fu et al., 1995;Taira et al., 1996); the activity of any one neuron is of course correlated with some parameters better than others, but the presence of simultaneous correlations makes it impossible to determine what that neuron really encodes (Fetz, 1992). Is it then possible to make any coherent statement regarding the role of M1 in arm movement control, at least to a first approximation? Is there one parameter of motor behavior which, if encoded in the M1 population, would explain most available results? We argue that of all the parameters previously considered, the only one that can provide the basis for a coherent model of M1 is muscle activation. The present chapter summarizes a series of recent studies (Todorov, 2000a;Todorov et al., 2000;Todorov,
meaningful quantity. It is hard to imagine how the combination of abstract directional signals encoded in extrapersonal space, and proprioceptive feedback encoded in a muscle- or joint-related frame of reference, could be useful or even interpretable. Evolutionary considerations. One would expect the general function of the spinal cord in primates, to the extent that it is used, to be related to functions it performed long before any supraspinal systems existed. The supraspinal systems most likely evolved to enhance those functions, by modulating the state of the spinal interneuronal networks in a way that helped them solve whatever problems they were already trying to solve(Loeb, 1999). The interpretation of abstract motor commands was not among these problems. A descending signal which modulates the state of the spinal networks should be compatible with the representation they use, and that representation is likely to be closely related to the pattern of muscle activation. Computational considerations. While the motor system is clearly hierarchical, the division of labor implicit in the directional coding hypothesis is unlikely. According to that hypothesis the vast supraspinal system is dedicated to the rather trivial calculation of a three-dimensional vector, while all the difficult control problems (e.g. coordinating redundant nonlinear actuators, dealing with complex interaction forces, ensuring stability in the presence of delays and uncertainty) are left to the ancient spinal circuitry.
1.2. Model overview
As explained above there are plenty of reasons to adopt a muscle-based model of M1. The question is, how can such a model account for the correlations with various behavioral parameters that seem unrelated to muscle activation? The main contribution of our work is in answering that question. Here we outline the features of the model and explain intuitively how they help us account for the many disparate results. M1 output. We focus exclusively on M1 neurons projecting to the spinal cord. Interneurons and output neurons projecting to other areas will not be modeled, since most of the arguments from Section 1.1 do not apply to them. The activity of such neurons can be systematically different, and so by focusing on the M1 output we avoid some of the response diversity that a more complete model would have to explain. M1-to-behavior mapping. We are only interested in the instantaneous activity of the M1 output neurons, and its causal relationship to the instantaneous motor behavior at some later point in time. All predictions regarding M1 will be derived by "inverting" that causal relationship (see Fig 1A). Motor planning, sensorimotor transformations, and transcortical feedback corrections are all enclosed in the black box labeled "brain" (Fig 1A) and do not concern us here. The "brain" receives online feedback,
takes into account all goals, plans, expectations, preferences, etc. and somehow generates the M1 output. This pathway will not be used to draw any predictions, and so we do not need to model it. Population analysis. Our goal is to explain average population (or subpopulation) responses rather than individual response types (Kalaska et al., 1989). Because of the redundancy in the M1-to- muscle mapping, the same pattern of muscle activity can be generated with an infinite family of M activation patterns. How the brain chooses one of them is an interesting question which we address separately in Sections 6 and 7. For the most part however we will not distinguish between M1 activity patterns that result in the same EMG pattern. Simplifications. We introduce a number of simplifications that allow us to avoid curve fitting and adhoc assumptions about poorly understood mechanisms (see Section 2). In particular, we use linear approximations to the spinal circuitry, muscle force production mechanisms, and kinematics of the multijoint arm. The effects of other descending systems are ignored, although there is evidence (Miller and Sinkjaer, 1998) that rubrospinal projections can be modeled similarly. Sensory feedback through spinal pathways is assumed to modulate the muscle activity generated by descending signals (see Section 7). While such feedback modulation is very important for achieving stability, its effect is washed out in the trial-averaged analyses used in the M1 literature, and so we do not model it explicitly here. Muscle state dependence. It is well established that for a fixed activation level, the force generated by a muscle varies with its length and velocity (Joyce et al., 1969;Zajac, 1989;Brown et al., 1999). The force-length-velocity surface at maximum (tetanic) activation is illustrated in Fig 1B, using the state-of-the-art muscle model of (Brown et al., 1999); details in Table 1 and Section 2. Note that for the same activation level and within the physiological range of muscle lengths and velocities, the force output varies from almost 0 to above maximum isometric force (defined as 1). Clearly the neural circuits controlling muscle activation have to take this state dependence into account if the resulting behavior is to resemble the desired one even remotely. The M1 signal compensating for the state dependence illustrated in Fig 1B plays a central role in our model since it correlates with position and velocity – and thus explains the observations that appear most problematic.
2. From M1 output to motor behavior
The relationship between M1 output and motor behavior is no doubt complex. If we incorporate every detail whose existence is currently suspected the resulting model will probably have enough free parameters to fit every available observation. The usefulness of a model, however, comes not from incorporating all proposed mechanisms and fitting all existing data, but from illuminating the underlying principles whose interplay gives rise to, or at least can give rise to the data. Our strategy is to avoid free parameters and the curve fitting associated with them. We do so with the help of simplifying assumptions:
descending activity may somehow be gated at the spinal level. Anatomically the M1 output is not the only descending source of activation: there are other descending systems, and the pyramidal tract itself originates from multiple cortical areas (Dum and Strick, 1991). It is not yet known, however, whether the pyramidal tract neurons in premotor and supplementary motor areas behave more similarly to the mixed population in their respective area, or to the pyramidal tract neurons in M1 (in which case the present model applies to them as well). Recordings from the red nucleus (Miller and Sinkjaer, 1998) indicate that Eq 1 may be a reasonable approximation to the effects of other descending systems. It should also be noted that lesions in M1 produce more profound deficits than lesions in any other supraspinal motor area (Johnson, 1992), further justifying our emphasis on M1. As discussed above and in Section 7, the spinal feedback contributions are assumed to average out over repeated trials.
2.2. From muscle force to endpoint force
Let f i ( ) t be the force (tension) that muscle i generates, and e i ( t )the corresponding endpoint force in the
two-dimensional workspace of the hand. While the direction of e i ( t )depends on joint configuration, that
dependence (see Fig 1A) is negligible for the small workspaces and stereotyped postures adopted in M experiments, so we will ignore it. Then
e i ( t ) = p ifi ( ) t (2)
where p i is a constant 2D vector whose length determines the scaling from muscle force to endpoint
force. Note that the present formulation ignores the effects of muscle forces causing redundant joint rotations. If we take this redundant subspace into account as outlined in (Todorov, 2000a), we could derive additional constraints on the M1 population activity (without affecting the constraints derived below). However M1 data has traditionally been analyzed in endpoint space, so these additional constraints could not be compared to existing experimental results.
2.3. Arm dynamics in endpoint space
Let M be the 2x2 matrix expressing arm inertia in endpoint space, f ext( t ) the endpoint force that the
hand applies against external objects, and x ( t )the hand position in endpoint space (the origin of the 2D
coordinate system is defined at the center of the workspace). We will assume that changes in inertia are negligible for the small movement extents and velocities used in the M1 literature, and approximate the arm dynamics with an anisotropic point mass. Making explicit the dependence of muscle force on
activation, length and velocity (the latter two being functions of x and x ^ ), the endpoint dynamics is
f ( a ( t ) ( ) t ( ) t ) ( ) t M ( ) t
i i^ i i ∑ p^ ,^ x , x ^ = f ext + x (3)
We now encounter a redundancy problem: there are many more muscle activations ai (and even
more cell activations c j ) than constraints given by Eq 3. We can proceed in two ways: a) linearize the
left hand side of Eq 3 to obtain a first-order approximation to the M1-to-behavior mapping; b) introduce further assumptions that resolve redundancy and allow us to use the detailed muscle model illustrated in Fig 1B. We will pursue the first approach (described in Section 2.4A) throughout this chapter. The second approach (described in Section 2.4B) will only be used when conclusions reached on the basis of the linearization might seem sensitive to details of the muscle model.
2.4A. Using a linearized muscle model
Assembling all a i ’s in the vector a , all c j ’s in the vector c , all f i ’s in the vector f^ , all wij ’s in the
matrix W , and all p i ’s in the columns of the matrix P , Eq 3 can be written in matrix form
P f ( a , x , x ) = f ext+ M x (3a)
Linearizing f around the baseline f ( 0 , 0 , 0 ) = 0 yields an approximation to the left hand side of Eq 3a:
P f ( a , x , x ^ ) = PS a − K x − B x . Here K and B are the endpoint stiffness and damping by definition, and
the scaling matrix
a
f
S =∂ is diagonal because each muscle is only affected by its own activation. From
Eq 1 we have PS a = PSW c , and so PSW c = f ext + M x + B x + K x.
We will now approximate PSW with FU , where the jth column vector u j of the matrix U is
a 2D unit vector pointing in the same direction as the jth column vector of PSW. Thus the j th neuron
activates a linear combination of muscles generating endpoint force in the direction u j. The 2x2 matrix
F = PSWUT^ ( UU T )− 1 , chosen to minimize the norm of the approximation error PSW − FU , can be
thought of as fitting an ellipse to the direction-dependent variation of the lengths of the column vectors in
PSW. Since the scaling from firing rates to forces is arbitrary, we can assume without loss of generality
that F = 1. Summarizing the derivation so far
FU c ( t − d ) = f ext ( ) t + M x ( ) t + B x ( ) t + K x ( ) t (4)
The 2x2 matrices M , B , K , F all reflect the geometry of the arm, in roughly the same way.
Consider a 4-degree-of-freedom arm, with the elbow pointing down and the hand constrained to a
c ( t d ) ( F ( ) t m ( ) t k ( ) t ) b Tj ( ) t
Tj
j f x x u x
− = u^ − + + +
2 1 ext (6)
This equation is essentially the dot-product of u j with the right hand side of Eq 5 (resulting in cosine
tuning). The reason we treat the damping term differently is that muscle damping is known to be
asymmetric – predominantly present during shortening (see Fig 1B). The truncated cosine term u T j x ( t )
implies that neuron j compensates for damping only when the hand is moving within 90 deg of its force
direction u j (corresponding to shortening of the muscles it activates).
Note that Eq 6 prescribes an “ideal” response which does not capture the diversity of real neuronal responses (Kalaska et al., 1989), and should be thought of as the average response of all cells
whose force directions are close to u j. To obtain an explicit model of response variability, we replace the
endpoint impedance parameters m , b , k in Eq 6 with cell-specific m j , bj , kj. These constants are
sampled independently from uniform distributions between zero and twice the average value:
m j ∈ [ 0 ; 2 ], b j ∈[ 0 ; 20 ], k j ∈[ 0 ; 100 ]. Multiplying by 2 and adding a cell-specific baseline cj
sampled uniformly from c j ∈[ 0 ; 34 ]yields an average response of 45 Hz in the preferred direction and
5Hz in the opposite direction, which is close to the range observed experimentally (Crammond and Kalaska, 1996). Eq 6 can also be used to analyse 1-degree-of-freedom movements, or to compute the average
firing of all neurons in their preferred direction ( u A ) and the direction opposite to it ( u N ). In either case
there are only two directions of interest ( u A and u N ), which simplifies Eq 6 to
c ( t d ) ( f ( ) t mx ( ) t kx ( ) t )
c t d f t mxt kxt bxt
N
A
ext
ext
2.4B. Using a detailed muscle model
We now return to Eq 3 and rederive the M1-to-behavior mapping using the detailed muscle model illustrated in Fig 1B. That model is identical to (Brown et al., 1999), except for the effects of sag,
yielding, low-pass filtering and tendon elasticity which we ignore. The muscle force f ( a , l , v )given in
Table 1 is a function of activation a , length l , and velocity v. All parameters are set to the average values for fast-twitch and slow-twitch muscles, which account for data from a wide range of experiments (Brown
et al., 1999). The activation a is normalized so that a = 0 corresponds to the recruitment threshold and
a = 1 to the frequency at which the force output saturates. Length l and velocity v are expressed in units
of l 0 and l 0 /sec, where l 0 is the length at which maximum isometric force is generated. Force output
f is in units of maximum isometric force.
The results from this section will only be applied to movement tasks where the directional asymmetries cancel each other (see Eq 5), so Eq 3 becomes
e f ( a ( t ) a l ( ( ) t ) v ( ( ) t )) m ( ) t
i i^ i i i i i max ∑^ p^ + , x , x ^ = x (8)
The vectors p i now have unit length, li and vi denote the length and velocity of muscle i , a i is the
baseline activity corresponding to posture at the center of the workspace, and e max is the maximum
isometric force each muscle can contribute in endpoint space. The baseline a i is needed here because the
muscle model is nonlinear and requires the total muscle activation. The direction of hand movement for which muscle i shortens most rapidly is very close to the
direction p i in which it produces endpoint force, and its orientation varies little over a small workspace
(see Fig 1A). We will approximate^3 this direction as being constant and equal to the (constant) line of
action p i. Then the length li is proportional to the projection p i T x of the endpoint position x on p i.
The relationship between endpoint and muscle kinematics becomes
p x
p x
i iT^
i i iT
v h
l r h
The negative sign is due to the fact that muscles shorten when the hand moves in the direction in which
they produce force. The constants ri^ specify the muscle lengths when the hand is at the center^ x =(^0 ,^0 )
of the workspace, and h is a scaling constant that relates endpoint distance to muscle length.
To resolve the redundancy in Eq 8 we will again use the assumption of cosine tuning. The total
force m x that needs to be generated will be distributed among individual muscles as
disi ( ) x ^ (^) ne^1 co 2 m p iT x max
where n is the number of muscles, the cocontraction level co is the scalar sum of the magnitudes of
endpoint forces contributed by individual muscles, and the dot-product p i T x corresponds to cosine
(^3) This approximation is very similar to the linearization used in (Mussa-Ivaldi, 1988).
model predicts that such a list should indeed accumulate, because muscle activation is correlated with all the candidate representations that have been considered.
3.1. Position
The postural activity of individual neurons predicted by the model (Eq 6) is illustrated in Fig 2A, for 8 hand positions arranged in a circle around the center of the workspace. The gradient of the response surface is oriented along the neuron's preferred direction, and the slope is determined by the amount of stiffness that needs to be compensated. Such monotonically increasing, and roughly linear postural responses have been well documented (Georgopoulos et al., 1984;Kettner et al., 1988).
3.2. Load
The predicted M1 activity in isometric force conditions (Fig 2B) is a linear function of the force being generated. In the 1-dof case (Eq 7) the response increases linearly with force magnitude for neurons
whose force vector u is within 90 deg of the isometric force f ext, and decreases otherwise. In the 2-dof
case (Eq 6) the response varies linearly with the cosine of the angle between u and f ext. Such
monotonically increasing, and roughly linear load responses have been well documented (Fetz and Cheney, 1980;Kalaska et al., 1989).
3.3. Movement
The average activity of all M1 neurons in their preferred and opposite directions (Eq 7) is illustrated in Fig 2C, for a typical 10 cm movement with 500 msec duration and bell-shaped speed profile. The M signal predicted from the linearized muscle model (left) is the sum position, velocity, and acceleration terms weighted by the direction-averaged stiffness, damping, and inertia. When the detailed muscle model is used instead (Eq 11) the individual contributions of these terms can no longer be separated. The composite M1 signal however remains similar (right). The predicted movement response closely resembles the speed profile and varies with direction – explaining why the activity of M1 neurons has often been interpreted as an encoding of movement direction (Georgopoulos et al., 1982) or velocity (Schwartz, 1994). In our model this phenomenon occurs simply because the damping that needs to be compensated dominates the other terms making up the M response. Note that the time-varying fluctuations relative to baseline are different in the preferred and opposite directions, which is a direct consequence of the fact that muscle damping is asymmetric. The asymmetry seen in Fig 2C has been observed experimentally (Kalaska et al., 1989;Crammond and Kalaska, 1996;Moran and Schwartz, 1999b). In comparison to the predicted response during an isometric
force ramp, the signal in Fig 2C has a more elaborate temporal structure – in agreement with experiments comparing movement and isometric force conditions (Sergio and Kalaska, 1998). One might wonder why the predictions in Fig 2C do not resemble the classic triphasic burst pattern of muscle activity. The reason is that this pattern is typical for large ballistic movements (Brooks, 1986), which are much faster than the movements studied in the M1 literature. If we apply our model (Eq
7) to a movement which is x times faster, the magnitude of the position component remains the same, the
velocity component scales by x , and the acceleration component scales by x^2. Thus for faster
movements the model predicts that the acceleration component of the M1 signal will become more pronounced, causing PV reversals consistent with the triphasic burst pattern. One might ask whether maintaining a position away from the center always requires a substantial increase of agonist activity (as shown in Fig 2C). This is indeed a side-effect of the linearized model where muscle stiffness is constant, but is no longer true in the detailed muscle model. As Fig 1B shows, the stiffness (i.e. partial derivative in the “length” direction) of realistic muscles is activation-dependent. Therefore the peripheral position can be maintained by substantially lower muscle activations. In the model, reducing muscle cocontration is accomplished by reducing the cocontraction parameter co in Eq 10 and 11, and the prediction is that the overall activity in M1 will also decrease. The relation between overall M1 activity and muscle cocontraction is further discussed in Section 3.6.
3.4. Movement with external load
If a constant load is added during movement, the model prediction (Eq 6) is the sum of the corresponding movement and load responses – as observed experimentally (Kalaska et al., 1989). The predicted response as a function of movement and load direction is shown in Fig 2D.
If an isotropic inertial load m ext^ is added during movement, the model prediction (Eq 6) is the
sum of the movement response and the term F −^1 m ext x ( t ). This added term accentuates the acceleration
dependent component of the signal shown in Fig 2C, especially in lateral directions where F −^1 is larger.
Thus our model predicts that when monkeys move with inertial loads of increasing magnitude, reversals of PV direction during the movement will begin to be observed – first in lateral movement directions. This is indeed the case (Kalaska et al., 1992).
3.5. Population vectors
Fig 2E shows the population vectors (Eq 5) for different postures, movements, loads, and movements in the presence of a constant load. In movement conditions the PV prediction has been averaged over the
model predicts in the absence of cocontraction. Thus increased muscle cocontraction is accompanied by a non-specific increase of activity in M1. We also conducted an fMRI experiment with human subjects asked to cocontract their forearm and wrist muscles while maintaining posture. Preliminary results were described in (Todorov et al., 2000). In agreement with the monkey data, activation over M1 was significantly elevated in the cocontraction condition compared to rest in the same posture. Since subjects were not producing movements or isometric forces in any direction, the observed change in M1 activity cannot correspond to any directional command. This converging evidence from monkeys and humans provides strong evidence against the directional coding hypothesis, and fits naturally in our framework.
4. Trajectory reconstruction
Given the correlations between M1 activity and various behavioral parameters illustrated in Fig 2, it should be possible to reconstruct the movement trajectory more or less accurately from a population of M1 responses. The potential for neuroprosthetic applications has recently generated substantial interest in this topic (Wessberg et al., 2000;Helms Tillery et al., 2001). Here we propose a simple reconstruction method that takes advantage of the multiple correlations. We also explain in the framework of our model why trajectory reconstructions based on the Population Vector hypothesis (Schwartz, 1994;Moran and Schwartz, 1999a) are systematically distorted.
4.1. Reconstruction methods
Since the M1 signal most closely resembles movement velocity (Fig 2C), the simplest viable method is to integrate the PV over time – yielding reasonably accurate trajectory reconstruction (Schwartz, 1994;Moran and Schwartz, 1999a). We applied this method to the PV predicted by our model (Eq 5) in the spiral movements used experimentally (1.5-7.5 cm spiral radius, traced in 2.5 sec according to the
2 / 3 power law). The reconstruction we obtained (Fig 3A) is very similar to (Schwartz, 1994;Moran and
Schwartz, 1999a), and displays the same distortion near the center of the spiral (the reason for the distortion is discussed later). An alternative method that has been studied is fitting a multilinear regression model (or a neural network) that maps cell firing rates to hand trajectories (Wessberg et al., 2000). These methods implicitly assume that cell firing contains information about a single movement parameter: velocity in (Schwartz, 1994;Moran and Schwartz, 1999a) and position in (Wessberg et al., 2000). Empirically this is not the case, and therefore it is possible to develop better methods taking advantage of the multiple correlations that have been documented (as well as the regularities of hand kinematics). Optimal reconstruction can be achieved by a Bayesian method which involves: i) a
generative model of how neural activity depends on arm state; ii) a dynamic model predicting the current arm state from the previous one; iii) a recursive estimator that computes the posterior distribution of the arm state given the distribution inferred previously and the current neural activity. Although steps towards building such a method are already being taken (Gao et al., 2001), this is a difficult estimation problem and the optimal solution is not likely to be found in the immediate future. Here we propose a simple reconstruction method which is related to PV integration (Schwartz,
- but also takes advantage of the correlations with hand position and acceleration (or load).
According to our model (Eq 5) the population vector U c (^ t − d )can be interpreted as the force acting on a
point mass that moves in a visco-elastic environment. Therefore we can reconstruct the hand trajectory by simulating such a point mass and driving it with the (instantaneously computed) population vector. The state of the simulation is our prediction regarding the state of the hand. As in the PV integration method, the initial state has to be known. Endpoint impedance and external loads also have to be known (or fitted to optimize performance). To assess the performance of our method, we generated a synthetic database containing 200 minimum-jerk trajectories (100 straight, 100 curved) with 300-600 msec duration and 5-10 cm extent. Each trajectory was executed in 4 dynamic conditions: normal, external load, external friction, external spring. For each execution we simulated the responses of 100 Poisson-spiking cells with mean firing rates given by Eq 6 (using the model of response variability), binned the spike data in 20 msec bins, and smoothed it with a 4th order Butterworth filter whose parameters were optimized separately for each reconstruction method. Fig 3B shows that the reconstruction error of our method (Dynamic Simulation) is lower than multilinear regression. As expected, prediction error decreases with the number of cells in the dataset. Note that our method generates simultaneous predictions about position and velocity, while the regression model is different in each case. We do not show an explicit comparison to the PV integration method because it is a special case of multilinear regression applied to velocity prediction. The reconstruction method described here may appear to be tailored to our model of M responses, but this is not quite true. It is a very generic method, applicable to any population that (for whatever reason) happens to correlate with position, velocity, and force. Such simultaneous correlations in M1 are an empirical fact, and therefore they can be taken advantage of by using similar reconstruction methods.
4.2. Variable timelag
The distortion seen in Fig 3A and (Schwartz, 1994;Moran and Schwartz, 1999a) is due to an interesting phenomenon described by (Schwartz, 1994). The timelag between the M1 population vector and hand
interpreted as encoding of direction irrespective of magnitude. After examining 6 such results (see below), we conclude that none of them is convincing evidence for direction encoding and in fact they are all consistent with our model. Note that the magnitude of both velocity (Schwartz, 1994) and force (Fetz and Cheney, 1980) is clearly being encoded in M1. So the question we are asking here is whether the encoding of direction dominates the encoding of magnitude, in a way inconsistent with treating either force or velocity as a vector.
5.1. Force magnitude versus direction
Let k denote the trial number, r ( ) k the magnitude of the 3D isometric force vector f ext( ) k produced by
the monkey, x ( ) k , y ( ) k , z ( ) k the coordinates of the unit vector which points in the same direction as
f ext ( ) k , and mfr ( ) k the rate at which a given neuron fires in trial k. To assess the relative contributions
of force magnitude and direction to the firing rate of that neuron, it might seem reasonable to compare the two regression models
mfr ( ) k b b x ( ) k b y ( ) k bz ( ) k
mfrk b brk
x y z
r
0
0
(D)
(M)
For 79% of the M1 neurons studied by (Taira et al., 1996) only regression model (D) provided a significant fit. Can we then conclude that 79% of the neurons in M1 do not encode force magnitude, and if we reach such a conclusion how can we reconcile it with earlier studies (Fetz and Cheney, 1980) showing the opposite? All we can conclude from the results of (Taira et al., 1996) is that the majority of M1 neurons do not encode force magnitude in the way prescribed by regression model (M). Note however that model (M) prescribes a rather unnatural encoding of magnitude: the response of the neuron is supposed to increase when the force magnitude increases, even when the force direction is opposite to the neuron’s preferred direction. Consider instead a more natural encoding (Eq 6) where the response of the neuron is simply the dot-product of its preferred direction and the force vector. Geometrically such a response is a plane (Fig 4A left). The regression models used by (Taira et al., 1996) attempt to fit the surfaces shown in Fig 4A (right) to M1 data. If M1 neurons encode force as we predict, regression model (D) would clearly provide a significant (although not perfect) fit, while regression model (M) would not fit at all. We have confirmed this using synthetic responses generated from Eq 6 (Todorov, 2000a).
5.2. Bias versus dynamic force
To assess the relative contributions of static and dynamic isometric force, monkeys were trained to maintain a static (bias) force in a specified direction and then produce an additional (dynamic) force pulse in the direction of a randomly chosen target (Georgopoulos et al., 1992). The finding that the preferred direction for dynamic force does not depend on the direction of bias force (and vice versa) was taken as evidence that bias and dynamic forces are somehow encoded separately in M1. This finding was illustrated with a plot almost identical to Fig 4B. Our model predicts the same effect (Fig 4B) without explicitly separating bias and dynamic force.
Consider a neuron whose firing mfr is given by the dot-product of its preferred direction u and the total
isometric force f total. Since f totalis the sum of f dynamicand f bias, we have
mfr = u T f total = u T ( f (^) dynamic+ f bias) = u T f dynamic+ u T f bias
If we now fix f biasand determine the preferred direction by varying f dynamic, the term u T^ f bias will simply
act as a baseline and the preferred direction will be u. If we change the direction or magnitude of f bias
(or remove it completely) and repeat the same procedure, the only difference will be in the baseline and
so the preferred direction will remain unchanged. The same holds for fixing f dynamicand computing the
preferred direction with respect to f bias.
5.3. Cell classification
Although it is clear that the responses of individual neurons are correlated with multiple movement parameters, one could classify neurons according to the movement parameter showing the strongest correlation. This was done in a 2D movement task (Ashe and Georgopoulos, 1994) where the time- varying mean firing rate of each neuron was correlated with target direction (which does not vary within a trial), hand position, velocity, and acceleration. The responses classified as direction-related were most common (47%) while the acceleration-related responses were least common (6%). Does this finding contradict our model which contains no explicit representation of movement direction? The answer turns out to be negative: the same classification scheme applied to Eq 6 (with response variability) yields up to 30% direction-related, and only 2% acceleration-related responses. How can a seemingly reasonable classification scheme find a large percentage of direction-related responses in a model without any explicit encoding of target direction? This is possible because a weighted sum of position, velocity and acceleration (Eq 6) does not necessarily resemble any one of these signals – as illustrated in Fig 4C. The family of synthetic response types given by Eq 6 is parameterized
