HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Full Text (PDF) All Versions of this Article: 97/5/3736    most recent 01064.2006v1 Alert me when this article is cited Alert me if a correction is posted Citation Map Services Email this article to a friend Similar articles in this journal Similar articles in ISI Web of Science Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via ISI Web of Science (3) Citing Articles via Google Scholar Google Scholar Articles by Trainin, E. Articles by Karniel, A. Search for Related Content PubMed PubMed Citation Articles by Trainin, E. Articles by Karniel, A.

Explaining Patterns of Neural Activity in the Primary Motor Cortex Using Spinal Cord and Limb Biomechanics Models

¹Department of Electrical Engineering, Technion—Israel Institute of Technology, Haifa; and ²Department of Biomedical Engineering, Ben-Gurion University of the Negev, Beer-Sheva, Israel

Submitted 5 October 2006; accepted in final form 4 March 2007

	ABSTRACT

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX A APPENDIX B GRANTS ACKNOWLEDGMENTS REFERENCES

 
What determines the specific pattern of activation of primary motor cortex (M1) neurons in the context of a given motor task? We present a systems level physiological model describing the transformation from the neural activity in M1, through the muscle control signal, into joint torques and down to endpoint forces and movements. The redundancy of the system is resolved by biologically plausible optimization criteria. The model explains neural activity at both the population, and single neuron, levels. Due to the model's relative simplicity and analytic tractability, it provides intuition as to the most salient features of the system as well as a possible causal explanation of how these determine the overall behavior. Moreover, it explains a large number of recent observations, including the temporal patterns of single-neuron and population firing rates during isometric and movement tasks, narrow tuning curves, non cosine tuning curves, changes of preferred directions during a task, and changes of preferred directions due to different experimental conditions.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX A APPENDIX B GRANTS ACKNOWLEDGMENTS REFERENCES

 
The primate motor system is a highly complex system leading to sophisticated motor activities resulting from the concerted activity of many cortical, sub-cortical and skeletal modules involving multiple feedback loops (Dum and Strick 2005). One of the key components in this system is the primary motor cortex (M1), which plays a major role in voluntary limb movement. Projections from M1 influence muscles through direct synapses onto motorneurons and indirectly through spinal inter-neurons (Kandel et al. 2001). It is now well understood that M1 is highly heterogeneous (Alexander and Crutcher 1990; Ashe 2005; Crutcher and Alexander 1990; Kakei et al. 1999) and that different populations in M1 represent different motor control signals and therefore cannot be uniformly interpreted. However, physiology provides evidence that some M1 neurons are highly correlated with muscle activity (e.g., Fetz and Cheney 1980; Morrow and Miller 2003). We refer to such neurons as muscle-related cells. Our model is aimed at explaining the functional behavior of these muscle related cells in M1 in the context of arm movement.

Many statistical models (e.g., Georgopoulos et al. 1982; Paninski et al. 2004; Sanger 1996) describe how neuronal activity of single neurons changes with hand variables and play an important role in neural prosthetic applications (Schwartz 2004). However, these approaches do not provide answers to the following basic questions. How do neural activities, projected from the motor cortex toward the spinal cord, result in hand movement and force? Given the redundancy of the controlled system (Bernstein 1967), how does the brain select specific control signals to achieve a motor task? The model provided here directly addresses these questions. In particular, our study is focused on the profound influence of some of the plant's properties on the control signal's behavior, where by plant we refer to the spinal-musculo-skeletal system. Understanding the spinal cord and the limb biomechanics is essential to explaining the motor control signals descending from the brain to the spinal cord. Our systems level physiological model describes the transformation from the neural activity in M1, through the muscle control signal (MCS), into muscle forces, joint torques and down to endpoint forces and movements. The redundancy of the system is resolved by adding biologically plausible optimization criteria.

Our model aims at being the "simplest" physiological model able to predict the complex patterns of activity of muscle related neurons in M1. Despite its relative simplicity, the model captures key features of the spinal cord and of the biomechanical system in the context of the tasks studied. In particular, limb mechanics, muscle unidirectionality and the finite time response of the muscles and the spinal cord are key properties of the system, responsible for basic features of the observed neural activity. Our model makes specific predictions about novel experimental situations and is consistent with current experimental results.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX A APPENDIX B GRANTS ACKNOWLEDGMENTS REFERENCES

 
Our model is aimed at predicting the neural activity of a subset of M1 neurons, which contribute linearly and additively to the input of the motorneuron pools (whether such a contribution is achieved through direct or indirect connections to motorneurons). We refer to such neurons as muscle related neurons and assume that these cells include the cortico-motorneurons and, possibly, some of the other neurons in M1. It is expected that many of the noncortico-motorneurons in M1 do not possess such a linear and additive contribution to the muscle control signal as they either project to other parts of the spinal cord or do not project to the spinal cord at all.

We consider voluntary movement tasks, implying that the control signal is initiated from the motor cortex. Moreover we consider simple trained tasks, motivating the assumption that the control signal is reproducible to a reasonable degree of accuracy, and that it is optimized for the task. Additionally, we assume well-trained tasks, which do not require the hand to achieve high accuracy or to retain stability against significant perturbations. This implies that the brain does not use significant co-activation of antagonist muscles to increase the hand's impedance (Hogan 1984, 1989; Osu et al. 2004).

Two essential properties related to the plant model are the number of degrees of freedoms (DOF) of the joints, and the number of muscles spanning the arm's joints. Arm models of varying complexity have been introduced in the literature, starting from 2 DOF with two muscles (e.g., Harris and Wolpert 1998) and at most 13 DOF with 42 muscle bundles (Garner and Pandy 2001). Our solution is general in the sense that it assumes neither a specific number of DOF nor a specific number of muscles. We use general formulae, in which D stands for the number of DOF and M stands for the number of muscles. Additionally we illustrate the model and its solution using a relatively simple arm model with 2 DOF (D = 2) and six muscles (M = 6).

The model addresses the following issues: 1) the relation between endpoint force and joint torques. 2) The relation between joint torques and joint movement. 3) The relation between equivalent muscle forces and joint torques. 4) The relation between MCSs and muscular force. 5) Resolving the redundancy of MCSs by an optimality criterion. 6) The relation between neuron activities and the MCS. 7) Resolving the redundancy of neuron activities by an optimality criterion.

Next, we present a detailed description of the model and refer the reader to Fig. 1A for a graphical explanation.

View larger version (19K): [in this window] [in a new window]   FIG. 1. Overview. A: general model block diagram. B: solution flow chart and the general model's parameters.

We assume that joint angle trajectories (or some equivalent data) are given by either direct measurements or by some model. The transformation between endpoint position and joint angles is given by

(1)

where P is a vector designating hand position and is a vector of joint angles. The function f is nonlinear and uses the lengths of the arm's links as parameters. It can be calculated based on multiplication of transformation matrices (Asada and Slotine 2000; Spong et al. 2006).

The kinematic equations of a bijoint planar arm (see Fig. 2A) are given in APPENDIX B.

View larger version (12K): [in this window] [in a new window]   FIG. 2. Arm model with 2 joints and 6 muscles. A: arm mechanics. L1 and L2 are the upper arm and forearm lengths, respectively. 1 and 2 are the shoulder and elbow angles, respectively. 1 and 2 are the shoulder and elbow torques, respectively. P = (x,y) are hand coordinates relative to the shoulder at (0,0). F is the force vector exerted by the hand. B: equivalent muscles. Sh, shoulder; Bi, bijoint; El, elbow; F, flexor; E, extensor.

Once joint angles and external force trajectories are given, it is possible to calculate the total torque exerted by the muscles at each joint. The joint torque is given by

(2)

where _S is the joint torque due to the external force and _D is the joint torque due to arm movement. These torques are derived using the following transformations and relations, which are standard in the robotics literature (e.g., Asada and Slotine 2000; Spong et al. 2006).

The transformation between endpoint force and joint torques (statics) is given by

(3)

where F^(ex) is the external force vector and J is the manipulator Jacobian. The manipulator Jacobian J is a N x D matrix, where N is the external space dimensionality and D is the number of DOF. It is derived from Eq. 1 according to

The Jacobian of a bijoint planar arm is given in APPENDIX B.

The relation between torques and movement (dynamics) is given by

(4)

where is a vector of joint angles, M accounts for inertial forces, C accounts for velocity-dependent forces, and G represents the influence of gravity. The derivation of the dynamical equation for a given arm model can be done using Euler-Lagrange equations (Asada and Slotine 2000; Spong et al. 2006). The explicit form of Eq. 4 for a bijoint arm without gravity (see Fig. 2A) is given in APPENDIX B.

Musculoskeletal geometry

An equivalent muscle describes an anatomical muscle, or a synergetic group of anatomical muscles, along with the spinal cord circuits related to these muscles. The moment arm describes the relationship between the equivalent muscle force and the torque produced at each joint. In our model, it is assumed that the moment arm is constant for all joint angles so that _ij(t) = R_ijF_i(t), where _ij(t) is the moment produced by equivalent muscle i on joint j and R_ij is the (constant) moment arm of equivalent muscle i acting on joint j.

The total joint moment at joint j is given by _j(t) = j = 1, 2, ... , D or

(5)

where R is the M x D matrix composed of the arm moments R_ij, and F(t) is a M x 1 vector composed of the muscle forces. Specific moment arm values for a 2-DOF arm model are given in APPENDIX B.

Equivalent muscle dynamics

In our model it is assumed, that the dynamic relation between the control signal u_i(t) and the muscular force F_i(t) is given by a simple first-order low-pass filter

(6)

where U(t) is a M x 1 vector of the MCSs u_i(t).

It should be noted that in our model, MCSs describe the population activity of cortical neurons contributing to the input of motorneuron pools. This population activity is neither an electromyographic (EMG) signal nor does it represent motorneuron activity although it is closely related to them. In particular, the relation between muscular force and its control signal is not identical with the dynamics of an anatomical muscle (or with a combined action of several anatomical muscles, in case of a simplified arm model) because the dynamics of equivalent muscles includes the influence of the spinal cord as well. Thereby we have used a relatively large value of the decay coefficient . The value of the decay coefficient was estimated to be 200 ms, leading to a good fit with experimental results. This relatively large value accounts for the overall response time of the spinal cord, muscle excitation, and muscle activation. Additionally, we've tested the sensitivity of the model's predictions to the value of and show that the general patterns predicted by the model are robust over a wide range of values of —see RESULTS.

From single neuron activities to MCS

In our model, each cortical muscle-related neuron contributes to one of the MCSs. We assume that the control signal is given by a sum of delayed neural activities. Specifically

(7)

where u_i(t) is the control signal corresponding to the ith equivalent muscle, n_ij(t) is the control signal (e.g., firing rate) corresponding to the jth neuron contributing to the ith equivalent muscle, w_ij is its weight, and d_ij is the latency of the neuron. We further assume that the latencies d_ij are distributed uniformly in the range of [50 ms, 100 ms], as the typical delay between M1 neural activity and muscle activity is 50–100 ms (Morrow and Miller 2003). These delays are long compared with the conduction times (10 ms) due to inward currents in the motorneurons (Morrow and Miller 2003). The neural signals {n_ij(t – d_ij)}_j=1^N_i are computed in APPENDIX A for any fixed weights and delays and are given by Eq. 11.

Redundancy resolution

The motor system is highly redundant (Bernstein 1967). Yet arm movements are highly stereotypical under a large variety of experimental conditions and generally vary little within and between subjects. Furthermore, variability tends to decrease with practice. Thus the redundancy of the motor system is resolved in a principled manner.

One may derive a criterion for redundancy resolution in several ways. One possible approach is based on the engineering principle of optimal control (Chow and Jacobson 1971; Engelbrecht 2001; Pandy 2001), which proposes that the system operates by optimizing a given cost function under the appropriate constraints. A second possible approach is to simply postulate an empirical redundancy resolution rule, verifying that it indeed explains and predicts observed phenomena. Finally, one may use self-organizing neural network models whereby the redundancy is resolved by the dynamics of these systems. Although each of these approaches possesses its merit, we have chosen to base our derivation on optimal control theory, which fits in nicely with an evolutionary view of the nervous system.

However, such an approach may pose some difficulties. The fact that the brain tries to reduce some costs and converges to some solution does not necessarily imply that it always reaches the optimal solution. The question arises as to whether biological neural networks actually solve optimization problems (for discussions, see Carpenter et al. 1987). The organization of the brain may impose further constraints on the motor plan. The studies of Morasso (1981), Shadmehr and Mussa-Ivaldi (1994), and Torres and Andersen (2006) suggest that the brain employs an independent path plan and thus restricts the possible muscle control signals. Similarly, motor primitives (Bizzi et al. 1991) may impose restrictions on the possible control signals. Yet the motor planning restrictions may be included in the constraints of the optimal control problem. For example, the stochastic optimal feedback control model of Burdet and Milner (1998) incorporates motor primitives as constraints.

We believe that the three redundancy resolution approaches alluded to in the preceding text are complementary rather than contradictory. The empirical approach focuses on the question "what," the optimal control theory approach focuses on the question "why," and the neural network approach deals with the question "how."

In the present study, we use optimal control theory based on simple optimality criteria. The redundancy of the system is resolved in two stages. First, we compute the MCSs according to a muscle optimization criterion and then compute how each MCS is divided between different neurons, according to a neuron-based optimization criterion.

Optimization criterion for the MCSs

A natural question that arises pertains to the cost function used. We can get a qualitative general answer to this question from observations of motor learning processes. For example, the studies of Burdet et al. (2001) and Franklin et al. (2003) show that when a human subject learns to reach a target through an unstable force field, two phenomena occur during the learning period. First, the percentage of successful trials increases and second, the EMG, and concomitantly, energy consumption, decreases. Such observations are supported by the clear evolutionary advantage of improving both task performance and saving energy. Furthermore, each of these two tendencies cannot be explained as a result of the other one. For example, unless a trained task requires high accuracy, the brain will not employ significant co-activation to achieve higher accuracy (Hogan 1984; Osu et al. 2004). The conclusion is that the cost function should combine both task performance and energy consumption (Hogan 1984; Miyamoto et al. 2004).

To resolve the MCS redundancy, we have used a quadratic optimization criterion, given by

(8)

where t₀ and t_f are the initial and final trajectory times. The simplifying factor of 1/2 has no influence on the solution.

Our choice of a quadratic optimality criterion is based on the following reasoning.

First, a squared neural control input is a simplified model of muscle energy consumption. Under isometric force conditions and less than 30% of maximum contraction, the metabolic power consumption of the muscle is proportional to the squared neural control input (Hogan 1984).

Second, some studies have demonstrated the predictive capability of such a criterion. For example, the study of Bolhuis and Gielen (1999) demonstrated that a quadratic cost function successfully predicted EMG activity of arm muscles under conditions of isometric static force.

Third, as mentioned in the preceding text, in general, the cost function combines both task performance and energy consumption. In the present case, we assume well-trained tasks, which do not require the hand to achieve high accuracy or to retain stability against significant perturbations. In this case, using a minimum energy criterion is reasonable.

Optimization criterion for neural activities

We assume that neural activity at the level of the single cell is selected in a way that minimizes the randomness of the muscle control signal. Our motivation here is that single neurons are noisy elements, and noise reduction through averaging may be one reason for using a large number of neurons for controlling the muscles. We assume here an optimization criterion, which is related to a measure of MCS stochasticity.

The criterion is based on a very general assumption regarding the way in which a single neuron contributes to the stochasticity measure of the muscle control signal. In particular, we assume that

(9)

where f(x) is a strictly convex function. Examples for functions that obey the preceding restriction are f(x) = x^r where r > 1. A function that does not obey this restriction is f(x) = x. The restriction on f(x) ensures that the optimality criterion will impose a single solution. The exact shape of f(x) does not influence the solution—see APPENDIX A. As a specific example, we may assume that f(x) is the variance of x.

Next we present the full optimization problem leading to the computation of the MCS and the neural activity.

Complete optimal control problem

(10A)

(10B)

(10C)

Then for each i = 1, 2, ... , M

(10D)

(10E)

Solution

The preceding optimal control problem is solved analytically in APPENDIX A, leading to the optimal neural activity

(11A)

where the constants _ij are computed in APPENDIX A. The optimal MCS is given by

(11B)

where the constant matrix R^(k) contains a subset of the rows of R as explained in APPENDIX A. For an isometric task (see APPENDIX A), the MCSs are given by

(11C)

Under the assumption of generalized symmetry (see APPENDIX A for details) a simplified expression results

The solution flow chart is summarized by Fig. 1B. The variables appearing in the model are summarized in Table 1. A detailed solution with specific parameters and arm configuration for the 2-DOF arm model is given in APPENDIX B.

We refer to a variety of experimental results published in the literature and present comparative model results. In particular, we provide a detailed comparison of our model's predictions with the experimental results of Sergio and Kalaska (1998) and Sergio et al. (2005). Our comparisons to the experimental results used only published data and figures.

We describe briefly the experimental setup of Sergio and Kalaska (1998) and Sergio et al. (2005). Two juvenile rhesus monkeys (a.k.a. Macaca mulatta) were trained to perform an isometric task and a movement task. The activity of single cells in the caudal part of M1 was recorded during both tasks. A neuron was selected if it was related to movements of the shoulder and/or the elbow but not to more distal joints and displayed directional tuning in at least one of the tasks. Most neurons were related to the shoulder and shoulder girdle, with a smaller number related to the elbow. The recorded neurons were in the caudal part of M1 of three hemispheres from the two monkeys (71 and 17 cells from the left and right hemispheres of monkey A, 44 cells from the right hemisphere of monkey B). Attempts were made to record neurons from all cortical layers, but the need for stable isolation over long periods of time led to a bias toward neurons in intermediate cortical layers. The set of recorded neurons included 132 cells. A particularly interesting aspect of the experiment is that the activity of the same neurons was recorded for the two different tasks (isometric and movement), enabling a detailed study of the (considerable!) change of activity of the same cell between tasks.

In an isometric task, the subject retains a fixed endpoint position in the face of an external force field. The dynamic profile of absolute force is the same in all trials, while force direction changes between trials. The force in this case was assumed to be F^(ex)(t) = F(t)[cos() sin() 0]^T = F(t)D(). The monkey held a static handle during 1–3 s (center hold time) and was then required to exert a ramp force of 1.5 N in one of eight directions. Force directions were spaced at 45° intervals, starting at 0°. We've estimated the shape of the force ramp, according to graph 1A in Sergio and Kalaska (1998) as a second-order spline with rise duration of 150 ms. The monkey was required to produce a force trajectory, which was confined to a horizontal plane (a significant vertical force resulted in an error). The handle was positioned in front of the monkey while hand location was at the midline, 20 cm in front of the sternum.

In the movement task, the monkey was required to push a load of 1.3 kg by 8 cm. Movement duration was 0.6 s, and the directions of movement were spaced at 45° intervals, starting at 0°. Because we did not have experimental data about the hand trajectory in the movement task, we assumed a minimum jerk trajectory (see APPENDIX B). The inertial force during movement was calculated according to the relation F = ma. After movement, the monkey was required to exert a force of 1 N against the pendulum. We've assumed that the force profile following the movement period is described by F(t) = 1 – ae^–bt, where a and b are chosen in a way that retains the continuity of the force and its first derivative at the end of movement.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX A APPENDIX B GRANTS ACKNOWLEDGMENTS REFERENCES

 
Intuition

We provide here a brief summary of the model's predictions, and their relation to the model assumptions described in METHODS. We discuss both the directional behavior and the response through time.

A temporal tuning function maps a time-dependent external force or movement direction into a related signal, e.g., the temporal firing rate of an M1 cell. A weighted sum of several cosine functions is a cosine function [i.e. a cos( – ) + b cos( – ) = c cos( – )]. In the isometric force task, each one of the external force components changes with direction as a cosine function, and therefore a linear relationship between endpoint force and the neural activity implies cosine tuning. Under isometric conditions, the relation between endpoint force and joint torques is linear (Eq. 3). However, due to the muscle unidirectionality (Eq. 6), the relation between the MCS and joint torques is linear for positive values and zero otherwise (Eq. 11B). Therefore the tuning function, predicted by the model under isometric conditions, is a truncated cosine rather than a full cosine.

The preferred direction (PD) is the direction corresponding to the peak of the tuning function. Since the linear relation between endpoint force and joint torques depends on the Jacobian (Eq. 3), the PD changes with arm posture.

The truncated cosine was obtained due to the linear relation between endpoint force and joint torques (Eq. 3) under isometric conditions. Yet in the movement task, the nonlinearity of arm dynamics (Eq. 4) leads to noncosine tuning at certain periods of the movement task.

Let us assume that the time response of the control signal simply follows the joint torques. In the isometric task, we would expect a positive or negative step response, and in the movement task, we would expect a positive step followed by negative one or a negative step followed by positive one (depending on movement direction). Because muscles are unidirectional (Eq. 6), negative pulses are impossible and negative torques are produced by the antagonist muscles alone. Taking this consideration into account, we would expect a nonsymmetric response. Thus in the isometric task we would expect a step response or no response, depending on force direction. In the movement task, we would expect a pulse-pause response or a pause-pulse response, depending on movement direction. We denote the latter pattern as biphasic. However, we should take into account the finite time response of the equivalent muscle, expressed by the decay coefficient (Eq. 6). Due to the finite time response, we obtain a pulse-step rather than a step response in the isometric task and a triphasic rather than a biphasic response in the movement task. We elaborate on these issues in the following sections.

Isometric task

DIRECTIONAL BEHAVIOR. In this section we test different properties of the predicted neural command signal, given by Eq. 11.

Tuning Function.    Equation 11C shows that the shape of the resulting tuning function of all MCSs is half cosine. If gravity is considered, the tuning function is a truncated cosine with a width, which is wider or narrower than 180°—see APPENDIX A. An example of the tuning functions for the 2-DOF arm model is given in Fig. 3A. From Eq. 11A, we expect that the tuning function of a single neuron should also be a truncated cosine function. Such a tuning function is narrower than a full cosine. How does the prediction of a narrow tuning function match the experimental data? Figure 3D shows an example of a representative neuron from Sergio et al. (Sergio and Kalaska 1998; Sergio et al. 2005) for which the width of the tuning function is 180°. The tuning functions of 51% of the cells in Amirikian and Georgopoulos (2000) were considerably narrower than 360°. According to Paninski et al. (2004), about a third of the neurons have narrow tuning function.

View larger version (33K): [in this window] [in a new window]   FIG. 3. Directional tuning of the muscle control signals. A: tuning functions for the 2 degrees of freedom (DOF) arm model. Right hand position is (–5 cm, 20 cm). B: preferred directions for the 2-DOF arm model. Parameters and notations are the same. C: preferred direction (PD) at different locations of the right hand for a neuron related to the bijoint flexor according to the 2-DOF arm model (APPENDIX B). The position of each arrow corresponds to the relative location of the hand on the planar work surface with the top arrow corresponding to the most distal hand location. The central location is (–5 cm, 20 cm). The remaining 8 hand locations were at a distance of 8 cm, spaced at 45° intervals, starting at 0°. The arrow represents the PD with respect to the isometric force. D: cell directional response in the isometric task, based on Fig. 1a from Sergio and Kalaska (1998) with permission. Discharge pattern of a shoulder-related M1 cell during the isometric force task. Each raster illustrates cell activity during 5 trials, and the raster location corresponds to the direction of the force. Data are aligned on the 1st significant force change, denoted by a solid vertical line (M). For each trial, the heavy tick mark to the left of the cursor movement onset line shows the time of target onset and the heavy tick mark to the right shows the time at which the final static level of force within the peripheral target was attained. It can be seen that after the force change (M) the response is about zero in 4 of 8 directions. Therefore the tuning function width is 180°.

Hand location and arm posture dependence of directional tuning.    According to Eq. 11C, the PD depends on hand position and hand location through the Jacobian. Figure 3C shows an example of PD dependence on hand location in the 2-DOF arm. Systematic changes in directional tuning due to hand location or arm posture have been observed in several studies (Ajemian et al. 2000; Caminiti et al. 1991; Hocherman and Wise 1991; Lacquaniti et al. 1995; Scott and Kalaska 1997; Sergio and Kalaska 1997, 2003). According to Eq. 11B of our model, neural directionality is defined in [(t) + '(t)] space. We elaborate on this issue in the DISCUSSION.

DYNAMIC BEHAVIOR.    Average response.
In the isometric task, different tuning functions differ only in their PD. Therefore according to Eq. 11C, the neural activity, aligned to the neuron's PD, of all cortical muscle related neurons is given by

(12)

where Q(t) is the aligned neural response and _CM is the force direction relative to the neuron's PD. This result is very robust because it depends only on a single parameter—. Furthermore, this result does not depend on the specific arm model used.

Equation 12 can also be used to describe the average of aligned neural responses, where the average is taken over many neurons. In practice, such averaging procedures suffer from two problems. First, different neurons have different latencies (Eq. 7). A possible way to resolve this problem is to estimate the latency of each neuron and align its response with respect to time before averaging. A neuron's latency can be estimated by computing the cross-correlation between its activity and the EMG signals (Morrow and Miller 2003). However, because the experimental data we've examined is based on a simple average, we had to compensate for the effect of different latencies by smoothing the predicted averaged response. We've smoothed the predicted averaged response by assuming that the latencies are distributed uniformly in the range of [50 ms, 100 ms].

Second, in the calculation of the experimental averaged response, each neuron's response is not exactly aligned to its PD because in practice the PDs are given within a resolution of 45°. This low resolution creates a smoothness effect in the direction axis as well. To compensate for this, we've smoothed the predicted averaged response by assuming that the directionality alignment is distributed uniformly in the range of [–22.5°, 22.5°].

   Predicted control signal.
If the control signal was defined merely relative to the muscle force, we would get a signal with similar shape to the muscle force profile. However, the control signal is also related to the first derivative of the muscle force (see Eq. 6). Therefore instead of obtaining a step response, we get a pulse step response—see Fig. 4. The pulse-step result is robust with respect to a wide range of values of the force derivative coefficient () as shown in Fig. 4C. The value of only determines the relative height of the pulse in the pulse-step response.

View larger version (9K): [in this window] [in a new window]   FIG. 4. Pulse step response in the isometric task for the general model. The ordinate is normalized to arbitrary units. A: step of muscle force. B: 1st derivative of muscle force. C: combined pulse step response. Solid bold line, = 100 ms; solid line, = 200 ms; dotted line, = 300 ms.

View larger version (16K): [in this window] [in a new window]   FIG. 5. Population activity for the isometric task for the general model. A: experiment, based on Fig. 9 of Sergio et al. (2005) with permission. Mean population response as a function of time and force direction, where the direction is relative to the PD of each cell. All data were aligned to the time of force onset (time 0) and the PD of each neuron was arbitrarily rotated to the right. B: simulation of the average neural response in our model. The ordinate of the simulation is normalized to arbitrary units.

   Single-cell activity.
The predicted single-cell response is proportional to the MCS response (see Eq. 11), i.e., a pulse step response multiplied by a half cosine tuning function. Notice that for the single cell, there is no need to apply a smoothing procedure and therefore Eq. 12 was used without further processing. In the experiment of Sergio and Kalaska (1998), 28% of the cells demonstrated a dynamic response similar to the response predicted by our model. See Fig. 6 for a comparison between this subset of neurons and the predicted response. The time response of most other cells was similar to part of the predicted response: 36% of the cells displayed a step response at their PD, 29.3% displayed a pulse response, and 6.7% of the cells were unclassifiable. There are several possible reasons why the predictions of our model with respect to single-cell activity apply only to a subset of the neurons reported on in the experiments. We elaborate on these in the DISCUSSION.

View larger version (16K): [in this window] [in a new window]   FIG. 6. Isometric task: single-cell response for the general model. A: experiment, based on Fig. 1a from Sergio and Kalaska (1998) with permission. Discharge pattern at the PD (top) and at the opposite direction (bottom) of a shoulder-related M1 cell in histogram format (10-ms bins). Data are aligned on the 1st significant force change, denoted by a solid vertical line (M). B: simulation of our model at the PD (top) and at the opposite direction (bottom). We assumed that the control signal is delayed by 100 ms. The ordinate of the simulation is normalized to the same arbitrary units as in Fig. 5B.

ISOMETRIC FORCE APPROXIMATION. Let us assume an isometric force experiment, in which the force profile is given by F(t) = ma(t) where m is the mass of the external load and a(t) is the acceleration in the movement experiment. We refer to this approximation as the isometric force approximation. There are two differences between this approximation and an accurate description of arm mechanics: 1) hand position, and consequently the Jacobian, changes during movement. Yet our simulation with the 2-DOF arm model and 8 cm movement distance shows that the assumption of fixed hand position has a minor effect on the tuning function. 2) Besides moving the load, the muscles also need to move the arm itself. This effect is calculated using the dynamics (Eq. 4). We discuss the effect of this component later. In the meantime we assume that this effect is insignificant.

Under these simplifying assumptions we can calculate the control signal's joint torques in a similar way to the isometric task.

TRI-PHASIC RESPONSE. Figure 7, A and B, shows graphs of F(t) and of F(t) + F'(t), for the isometric force approximation of the movement task. As shown, there is an important difference between the movement task and the isometric task with force ramp profile. Unlike the force ramp profile scenario, the expression F(t) + F'(t) is not always positive. The profile of F(t) + F'(t) is triphasic, where the first phase is positive, the second is negative, and the third is positive again. We refer to this type of response as triphasic as it relates to a known phenomenon in which the activation period of one muscle is followed by an activation period of the antagonistic muscle, followed by a further activation period of the first muscle. We refer to the transition times between the first and the second phases and between the second and third phases as transition times. If we examine two antagonistic muscles in the 2-DOF arm model, for example, the shoulder flexor and the shoulder extensor, it is clear that the predicted pattern is triphasic. The triphasic response is a pattern of flexor-extensor-flexor, or alternatively extensor-flexor-extensor, depending on the movement direction. Note that the brain would not employ a triphasic pattern if not for the finite time response of the equivalent muscles. When the neural input to a muscle ceases, the muscular force does not vanish immediately because it has some decay period. To change the muscular force rapidly enough, the antagonistic muscle should be activated, while the agonist muscle is still active—see Fig. 7C.

View larger version (19K): [in this window] [in a new window]   FIG. 7. Tri-phasic response in the movement task for the general model. A: hand inertial force during a minimum jerk trajectory, given by , where (the units are arbitrary). B: . The positive phases (1st and 3rd) represent the profile of 1 muscle control signal (MCS). The negative phase (2nd) represents the profile of the antagonist MCS. C: shoulder muscle torques in the 2-DOF arm model for the movement task. Movement direction is 0. Dashed line, shoulder flexor torque; dotted line, shoulder extensor torque; solid line, total torque; vertical lines, transitions between antagonist control signals. D: influence of the value of (force derivative coefficient) on cell response at the PD. Solid bold line, = 100 ms; solid line, = 200 ms; dotted line, = 300 ms. E: same as D at the opposite direction.

POPULATION ACTIVITY. Figure 8 presents a comparison of the experimental results with the model's prediction. The simulation is based on Eq. 12, as we used the isometric force approximation. Recall that Eq. 12 is independent of the arm model, and therefore the comparison with Sergio et al. (2005) is appropriate. The average response in Sergio et al. (2005) was generated for those neurons that were directionally tuned during target hold time, using the PD from that epoch. As can be seen in Fig. 8, the temporal behavior of the neural signal predicted by the model corresponds to the experimental results. Particularly important is the fact that the results for both the isometric (Fig. 5) and the movement (Fig. 8) tasks reproduce the qualitative shift in neural activity between the two tasks.

View larger version (17K): [in this window] [in a new window]   FIG. 8. Population activity in the movement task for the general model. A: experiment, based on Fig. 9 of Sergio et al. (2005), with permission. Mean population response as a function of time and force direction, where the direction is relative to the PD of each cell. All data were aligned to the time of force onset (time 0) and the PD of each neuron, calculated at target hold time, was arbitrarily rotated to the right. B: simulation of a MCS in our model. The ordinate of the simulation of our control signal is normalized to the same arbitrary units as in Fig. 5B.

View larger version (23K): [in this window] [in a new window]   FIG. 9. Movement task: single-cell response for the general model. A: experiment, based on Fig. 1b from Sergio and Kalaska (1998) with permission. Discharge pattern at the PD (top) and at the opposite direction (bottom) of a shoulder-related M1 cell in histogram format (10-ms bins). Data are aligned on the 1st significant force change, denoted by a solid vertical line (M). B: simulation of the single neuron activity at the PD (top) and at the opposite direction (bottom). We assume a delay of 100 ms. The ordinate of simulation is normalized to the same arbitrary units as in Fig. 5B. C: reconstruction of Fig. 2C in Todorov (2002). Activity at the PD (top) and at the opposite direction (bottom).

View larger version (27K): [in this window] [in a new window]   FIG. 10. Reversals of the PD during movement task. A: based on Fig. 2b of Sergio and Kalaska (1998) with permission. The temporal trajectory of the PD of a representative cell during movement tasks. The trajectory was determined by a 50-ms sliding-window analysis. Time windows within which the cell was significantly related to direction are shown by an asterisk. Time windows within which the cell was not directionally related are shown by circle. Notice that around the transition times (characterized by reversal of PD) the cell was not directionally tuned. Large thick concentric circles denote movement onset and offset. B: simulation of a single neuron PD in the 2-DOF arm model. The PD was calculated as the angle at which the tuning function attains its maximum. The PD is relative to its value during target hold time. C: temporal tuning function in the 2-DOF arm model for the movement task at time = 80 ms.

To simulate this phenomenon, we cannot use the general equation for the isometric task (Eq. 12), as the calculations in this case depend on the arm model. We need to calculate the control signals according to Eq. 11B. However, the arm in this experiment was allowed to move in three-dimensional space, whereas we do not have data about the monkey's arm moments for a 4-DOF model. Therefore we have studied this phenomenon using the 2-DOF arm model as described in APPENDIX B, assuming that qualitatively, the kind of effect of dynamic nonlinearity around transition times would be similar for the 2- and the 4-DOF models. Our simulations (with the 2-DOF model) show that during these periods the tuning function becomes wider and looses its cosine shape—an example of a bimodal tuning function can be seen in Fig. 10C. The noncosine prediction is in agreement with the experimental results shown in Fig. 10A. During periods around the transition times, cell activity was not directionally tuned (time windows represented by circles).

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX A APPENDIX B GRANTS ACKNOWLEDGMENTS REFERENCES

 
We have presented a mechanistically explicit modeling approach based on the following principles. 1) It is possible to attribute functional meaning to neural activity in the motor areas in the brain. 2) Understanding the spinal cord and the limb biomechanics is essential to interpreting the motor control signals descending from the brain to the spinal cord. In particular, the motor cortex control signals are related to endpoint variables through dynamic models. 3) The redundancy of the control signals can be resolved using optimal control theory. 4) Neural activity characteristics, such as tuning functions, PDs and coordinate systems, can be derived within a general theoretical framework.

Based on this modeling approach, we have presented a model that explains muscle-related neural activity in the M1. The model has demonstrated good agreement with a variety of experimental results.

The model successfully predicted the temporal population activity in the isometric and movement tasks (pulse-step, tri-phasic) and for different force or movement directions (truncated cosine) as well as the qualitative nature of the transition between the two tasks. There are a few differences between the population predicted responses and the experimental results. The peaks in the experiment are wider and continue beyond the force ramp/movement duration, and the tuning function in the experiment is somewhat wider. Additionally, there is a 10 spike/s ground level in the experiment.

The predictions of the time responses of the single cell partially agree with the experimental results. A subset of the cells (28% in the isometric task, 44% in the movement task) displayed close similarity with the predicted response. Most other cells displayed similarity only to part of the response. Comparing population and single-cell predictions, a subset of the single cells displayed better similarity than the population response (narrow peaks, no ground level), whereas others displayed less similarity than the population response (e.g., single-cell pulse response in the movement task).

The gaps between model and experiment are possibly due to the model's simplifications and/or because only part of the neurons examined in Sergio and Kalaska (2005) were muscle-related cells. It is difficult to determine which of these two reasons is more important as we do not know which of the cells were muscle related. To better examine models of neural activity in M1, one should use recorded activity of neurons with similar connectivity, i.e., the connectivity of the measured neurons should be identified.

The prediction of single-cell tuning functions is supported by experimental results. The model successfully predicts the half cosine tuning function in Sergio and Kalaska (1998). This prediction also agrees with Amirikian and Georgopoulos (2000) and with Paninski et al. (2004) showing that a significant fraction of the M1 neurons have narrow tuning functions. Interestingly, the analysis of Paninski et al. (2004) shows that narrowness of the tuning function correlates with the information about hand trajectory conveyed by the cell. The model also successfully predicted noncosine tuning-functions around the transitions times of the movement task.

Although tuning functions and PDs do not possess an intrinsic functional meaning in the context of our model, we've calculated them to compare our results with published experimental data. The model predicts the dependency of the PD on arm posture and on hand location as well as the two reversals of the PD during the movement task.

In addition to explaining known results, one may try to test the model's predictions in new situations. For example, our simulations show that the phenomenon of unusual tuning functions around transition times becomes more apparent as the load becomes lighter. The reason for this is the nonlinear nature of arm dynamics, which becomes more significant at small loads, where the isometric force approximation is poor. This prediction has not yet been tested experimentally.

The model's relative simplicity and analytic tractability support a causal explanation of how each of the plant's properties influences the overall behavior. Such an explanation was given in RESULTS.

Our model is aimed at predicting the behavior of muscle-related cells. Yet many M1 neurons are not muscle related cells (Alexander and Crutcher 1990; Ashe 2005; Crutcher and Alexander 1990; Kakei et al. 1999). Thus it is possible that the neural activity of other neurons in M1 have very different interpretations, such as movement direction. One possible interpretation of cosine tuning (Georgopoulos et al. 1982) and of the population vector (Georgopoulos et al. 1983) is that M1 employs a coding scheme based on external space directions. However, correlations with high-level parameters are not sufficient to prove high-level representation. Cosine tuning and the existence of a population vector are determined by very general necessary and sufficient conditions rather than by any intrinsically coded coordinate system (Mussa-Ivaldi 1988; Sanger 1994). However, the analysis of Mussa-Ivaldi (1988) and of Sanger (1994) does not exclude the possibility of direction coding. It just shows that the cosine tuning and the population vector are not sufficient to determine whether neurons code directions in external coordinates or in some intrinsic coordinate system. To better model the neural activity, further information is needed. Indeed, dissociation experiments (e.g., Crutcher and Alexander 1990; Evarts 1968; Scott and Kalaska 1997; Thach 1978) provide further information, which enables one to distinguish between different interpretations of neural activity. Interestingly, the early dissociation experiment of Evarts (1968) showed that the discharge frequency of the majority of pyramidal tract neurons was related primarily to the force and to the force derivative and was only secondarily related to the direction of displacement.

Currently it is known that neural coding in the motor cortex is much more complex than described by cosine tuning and the population vector (Johnson et al. 2001). Neurons in M1 do not possess a single fixed PD. The PD depends on external variables (such as force or velocity) with which the neuron is correlated. The correlations change for different tasks, and also as a function of time within a single task. Many of the tuning functions in M1 are narrow, asymmetric, and sometimes bimodal (Amirikian and Georgopoulos 2000), and the population vector may deviate from extrinsic variables (Scott et al. 2001). Yet the question of distinguishing between different interpretations of neural activity is always relevant. Interpretations of neural activity are in fact models. Therefore it is possible to choose between different interpretations on the basis of the agreement with the experimental results, the simplicity of the model, and the relative weight given to each of these two factors.

The study of Scott et al. (2001) shows that the population vector may deviate from extrinsic variables. Our study further shows in RESULTS that the MCS provides a better interpretation of the population activity reported in the experiments of Sergio et al. (2005) compared with the interpretation of endpoint force.

Observed PD changes, due to arm posture, can be explained by intrinsic coordinates but not by external coordinates (Ajemian et al. 2000; Scott and Kalaska 1997). Our model further suggests that such intrinsic coordinates do not necessarily follow simple mechanical expressions, such as joint angular velocities (Ajemian et al. 2000; Scott and Kalaska 1997) or torques (Scott and Kalaska 1997). In fact, the coordinates might be hard to guess without a physiological model as even in our simplified model, the derived coordinates are composed of both joint torques and joint torque derivatives.

Several previous studies have presented models of neural activity in M1 and its functional interpretation. The early work of Humphrey (1972) described the relation between neural activity in the motor cortex and muscular torque by a dynamic linear model. We have followed the same approach. However, Humphrey (1972) addressed the forward control problem, whereas our model addresses the inverse control problem.

The model of Bullock and Grossberg (1988) and of Cisek et al. (1998) assumes that automatic processes convert the target position command and the movement's overall speed command into an arm trajectory. The focus of this model is quite different from ours. The model of Bullock and Grossberg (1988) and of