## Abstract

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

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

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. 1*A* for a graphical explanation.

### Hand position and joint angles

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. 2*A*) are given in appendix b.

### Joints torques

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* × *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. 2*A*) 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 _{ij}*

*F*

_{i}(

*t*), where τ

_{ij}(

*t*) is the moment produced by equivalent muscle

*i*on joint

*j*and

*R*is the (constant) moment arm of equivalent muscle

_{ij}*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* × *D* matrix composed of the arm moments *R _{ij}*, and

*F*(

*t*) is a

*M*× 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*× 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

*i*th equivalent muscle,

*n*(

_{ij}*t*) is the control signal (e.g., firing rate) corresponding to the

*j*th neuron contributing to the

*i*th equivalent muscle,

*w*is its weight, and

_{ij}*d*is the latency of the neuron. We further assume that the latencies

_{ij}*d*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 {

_{ij}*n*(

_{ij}*t*−

*d*)}

_{ij}_{j=1}

^{Ni}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*_{0} 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. 1*B*. 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.

### Experimental setup

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

### 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. 3*A*. 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.

#### PDs.

While different MCSs have tuning functions with similar shapes, they are distinguished by their PDs. For example, the PDs in the 2-DOF arm model are 355° for the shoulder flexor, 30°for the bijoint flexor, 87° for the elbow flexor, 175° for the shoulder extensor, 210° for the bijoint extensor, and 267° for the elbow extensor—see Fig. 3*B*. In fact, the PD of each MCS is in one of two opposite directions. For example, the PD of *u*_{1}(*t*) is 355 or 175° depending on the sign of the expression *F*(*t*) + α*F′*(*t*). Thus the PD of the control signal is *reversed* during the task whenever the expression *F*(*t*) + α*F′*(*t*) changes its sign. We will address this phenomenon in the movement task section. Elsewhere, we refer by PD to the PD given that *F*(*t*) + α*F′*(*t*) > 0.

#### 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 3*C* 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. 4*C*. The value of α only determines the relative height of the pulse in the pulse-step response.

To obtain an averaged response, the preceding profile is multiplied by a truncated cosine as described by *Eq. 12*. As explained in the preceding text, the averaged response is smoothed to account for the effects of the simple averaging. Recall that *Eq. 12* is independent of the arm model, and therefore the comparison with Sergio et al. (2005) is appropriate. The averaged response in Sergio et al. (2005) was generated from those neurons that were directionally tuned during target hold time using the PD from that epoch. The good qualitative match between the experiment and model can be clearly observed in Fig. 5.

###### Interpretation.

The population vector is defined as the sum of neuronal firing rates (relative to each cell baseline) multiplied by each neuron's PD unit vector (Georgopoulos et al. 1983). The population vector has been used to estimate output variables based on neural activity. In the experiment of Sergio et al. (2005), the population vector significantly deviates from the force direction as shown in Fig. 11 in Sergio et al. (2005). The description of the population activity as a control signal has two advantages over the description of the population activity by a population vector of hand force. First, the control signal can be approximated by a simple average. Second, the agreement between the MCS and the averaged neural activity is better than the agreement of the population vector with the external force. Indeed, modern extraction algorithms (e.g., Salinas and Abbott 1994; Shpigelman et al. 2005; Taylor et al. 2002; Wu et al. 2004) significantly improve the accuracy of the simple population vector. However, such techniques are much more complex than the simple average we have used. This does not reduce the importance of such algorithms for BMI applications. However, in the preceding case, “the population codes a MCS,” is the preferable interpretation of the population activity because it is the simplest interpretation that matches the experimental results.

###### 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.

### Movement task

##### 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. 7*C*.

Figure 7, *D* and *E*, show that the triphasic response is robust as it occurs for all values of α. The value of α only determines certain details of the shape of the curve. The first pulse gets narrower with the increase of α, the second pulse shifts to the left and the third pulse broadens with the increase of α. To understand the influence of α, notice that the force profile (Fig. 7*A*) is quite similar to a sine wave. Therefore the “phase shift” of *F*(*t*) + α*F′*(*t*) (Fig. 7*B*) becomes larger with the increase of α.

##### 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.

##### SINGLE-CELL ACTIVITY.

Figure 9 displays a comparison between a subset of neurons in the experiment and the predicted response. For the simulation, *Equation 12* was used without smoothing. In the experiment, 87.5% of the cells, displaying step or pulse-step response in the isometric task (56% of all the cells) displayed the burst-pause-burst response in their PD (Fig. 9*A*) as predicted by the model (Fig. 9*B*). Most of these cells (78 or 44% of all cells), as well as part of the cells in the remaining group, displayed the predicted response of pause-burst-pause in the opposite direction (Fig. 9, *A* and *B*); 45.4% of the cells, exhibiting a pulse response in the isometric task (13% of all the cells), displayed the burst-pause-burst response in their PD. Many of these cells displayed the response of pause-burst-pause in the opposite direction. Other neurons in the movement task displayed pulse response (20%), step response (13%), pulse-step response (4%), and the rest were unclassified. Similarly to the single-cell dynamic response for the isometric task, the agreement with our predictions is partial though significant. We elaborate on this issue in the discussion.

##### PD REVERSAL.

According to a preceding section, the PD of the control signal is reversed each time the expression *F*(*t*) + α*F′*(*t*) changes its sign. Therefore in movement tasks, the PD of the control signals is reversed during the second phase. According to the model, we expect to see such a phenomenon at the level of the single cell as well. Figure 10 shows that this prediction matches the experimental data. The graph of Sergio et al. (2005) was generated from the cell activity shown in Fig. 9, which represents 44% of the cells. Note that the PD goes through two reversals during the movement task, as expected from the model.

##### NON COSINE BEHAVIOR.

The control signal in our model can be decomposed into two components. One component is related to arm movement, and the second is related to the load movement. The isometric force approximation is based on the assumption that the torque needed to move the arm is small compared with the torque needed to move the load. This assumption holds only when the expression *F*(*t*) + α*F′*(*t*) is significant. However, when the load is insignificant and/or during periods around the transition times, the isometric force approximation is poor. The component related to arm movement is significant during periods around the transition times, because *F*(*t*) + α*F′*(*t*) ≈ 0. In the experimental setup, we've examined the load was significant. However, we still expect the isometric force approximation to be poor around the transition times. During these periods, the nonlinear arm dynamics becomes more significant. Consequently, half cosine tuning is not expected anymore.

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. 10*C*. The noncosine prediction is in agreement with the experimental results shown in Fig. 10*A*. During periods around the transition times, cell activity was not directionally tuned (time windows represented by circles).

## DISCUSSION

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 Cisek et al. (1998) does not use optimal control theory to explain why the brain produces a certain control signal. On the other hand, it introduces a computational scheme for M1. Similarly, the model of Baraduc et al. (2001) focuses on how the motor cortex learns and computes control signals, using a neural network model. The neural network in this model transforms position in extrinsic coordinates into position in joint angle coordinates.

The model of Todorov (2000, 2002) describes the average neural activity of neurons with a certain PD by a population vector related to the external endpoint force, endpoint position, velocity, and acceleration. We share with Todorov (2000, 2002) the general philosophy, according to which, understanding the spinal cord and the limb biomechanics is essential to interpreting the motor control signals descending from the brain to the spinal cord. However, our view and understanding of the spinal-muscular-skeletal system is very different from Todorov (2000, 2002). In fact, our model and Todorov's (2000, 2002) are based on very different assumptions. Most importantly, the model used in Todorov (2000, 2002) does not contain muscle unidirectionality, muscle and spinal cord finite time response, and realistic multi-joint mechanics (see, for example, Moran and Schwartz 2000)—features that we find crucial for achieving qualitatively correct predictions. Due to the conceptual difference between the models, the predictions of the two approaches contradict each other, even in the simplest scenarios. For example, in an isometric force task the average neural activity predicted by Todorov (2000, 2002) differs from our prediction. First, the predicted tuning function in Todorov (2000, 2002) is a full cosine, whereas our model predicts a truncated cosine tuning function. Recall that both narrow and wide tuning functions are present in M1; however, muscle-related cells have narrow tuning functions. Second, the response at the PD predicted by Todorov is a step, whereas our model predicts pulse-step response. The predictions of Todorov (2000, 2002) also differ from ours with respect to the movement task. For example, comparing Fig. 2C in Todorov (2002), and Fig. 9*B* in this work, with the results in Sergio and Kalaska (1998; Fig. 1b) shows the current model provides a significantly better qualitative description—see comparison in Fig. 9.

Finally, the model of Neilson and Neilson (2005) uses a quadratic optimality criterion to resolve the redundancy of the muscle control signals. The model assumes known flexion and extension torques, which are used in linear equality constraints. Given that the flexion and extension torques are known, one can use this assumption with our general solution for the MCS as this assumption simply doubles the number of linear equality constraints. Unlike our MCS model, the model of Neilson and Neilson (2005) is limited to static scenarios. Another important difference between this model and our model is that it does not include the inequality constraints due to muscle unidirectionality. In the case where all the predicted control signals are positive, the solution of this model is similar. In any other case, this model cannot be applied due to negative control signals. Thus our solution is more general.

Although we have been able to demonstrate good agreement with experimental results, many open questions remain for future research of which we mention a few. *1*) Cortical neurons projecting to the spinal cord have further functionalities beyond providing input to the motorneuron pools. For example, projections to Ia-inhibitory interneurons enable the brain to regulate the gain of the reflex loop and thus change arm impedance. A more comprehensive model should take such connections into account. *2*) To extend the model's applicability and predictive power, a more realistic spinal-biomechanical model is needed. For example, to apply the model to tasks involving impedance control, it should take into account neural noise, control of the reflex loop gains, influence of the MCS amplitude on muscle stiffness and viscosity and the dependence of arm moment on joint angles.

We should emphasize that this study provides but one step in a long journey toward the formidable task of constructing a unified controller-plant based model of the biological motor system. In general, the plant is very complex, the control policy is task dependent, and the diversity of tasks is large. Furthermore, the biomechanical system (Loeb et al. 2002) and the spinal cord (Burke 2004) are not fully understood. To combine the plant model with a controller model, while incorporating appropriate optimality criteria, one needs solve the associated optimal control problem. In general, this computational problem is very difficult due to its high dimensionality and the inherent nonlinearity and stochasticity. Finally, to better test such models, the connectivity strengths of recorded neurons to different motorneurons and interneurons should be identified. Clearly the goal of achieving a full computational understanding of voluntary motor control is a long distance away. Nevertheless we believe that the encouraging results of the simplified model we developed show that the task is not insurmountable.

## APPENDIX A

### General solution

From *Eqs. 10a* and *10b*, we obtain Thus we have replaced *Eqs. 10a* and *10b* with equality constraints that depend only on *U*(*t*). Now because the integrand of the cost function and the constraints depend only on *U*(*t*), the solution at a given time depends only on the integrand and the constraint at that particular time. Therefore the optimal control problem can be replaced by the following optimization problem This problem can be solved using Lagrange multipliers—see, for example, Bryson (1999).

We denote Lagrange multipliers for the equality constraints by λ_{j}(*t*), *j* = 1, 2, … , *D*, and the Lagrange multipliers for the inequality constraints by λ′_{i}(*t*), *i* = 1, 2, … , *M*.

At the minimum the equality constraints *R*^{T}*U*(*t*) = τ(*t*) + ατ′(*t*) are satisfied.

In addition, for each *i* = 1, 2, … , *M*, the solution satisfies one of the following two conditions: either the minimum is on the boundary of the admissible region (so-called active constraints) or the minimum is inside the admissible region (so-called inactive constraints), and then or Thus either We conclude that

Therefore the *D* equality constraints become *U*(*t*) = [*R*λ(*t)*]_{+} and therefore the set of equations *R*λ(*t) = 0* divides the space Λ, spanned by λ(*t*), into regions Λ^{(k)}, *k* = 1, 2, … , *K*. For each region Λ^{(k)}, there is a unique group of active muscle indices *m _{k}*, for which

*u*(

_{i}*t*) > 0,

*i*∈

*m*. More precisely Because

_{k}*R*

^{T}[

*R*λ(

*t*)]

_{+}= τ(

*t*) + ατ′(

*t*) is a continuous transformation from Λ space to space, each region Λ

^{(k)}in the Λ space corresponds to a region in the space.

Thus the solution is computed in the following way. Given that belongs to the *k*th region of the space, then where *R*^{(k)} is a matrix containing the rows of R, corresponding to the index group *m _{k}*. Hence

*Equation (*)*provides the general solution to the optimization problem. This solution can be further simplified under the assumption of

*generalized symmetry*, which is defined as independence of the expression

*R*

^{(k)T}

*R*

^{(k)}on the specific region

*k*. Under this assumption, the expression

*R*

^{(k)T}

*R*

^{(k)}can be computed once, and there is no need to find the region to which τ(

*t*) + ατ′(

*t*) belongs. In particular, the condition of generalized symmetry holds in the case of

*strict symmetry*, implying that muscles are composed from symmetric antagonist pairs. Namely, for each pair of antagonist muscles

*i*and

*i′*,

*R*= −

_{i′j}*R*,

_{ij}*j*= 1, 2, … ,

*D*. For a given region

*k*, either

*i*∈

*m*or

_{k}*i′*∈

*m*, and because

_{k}*R*

_{ij}*R*= (−

_{il}*R*)(−

_{ij}*R*) =

_{il}*R*

_{i′j}*R*,

_{i′l}*j,l*∈ (1, 2, … ,

*D*), the expression does not depend on the specific region

*k*. Thus under the assumption of strict symmetry, the expression R

^{(k)T}

*R*

^{(k)}can be computed once, while the index group

*m*can be any group that contains a single representative from each pair of antagonist muscles. The condition of generalized symmetry is less restrictive. For example, it suffices that for two groups of flexors and extensors surrounding the same joint, while not requiring each of the muscles in the group to have a symmetric pair.

_{k}### Neural activity

Having computed the MCSs, {*u _{i}*(

*t*)}

_{i=1}

^{M}, we now ask how the neural activity of each control signal is divided between the neurons contributing to this control signal. Because the integrand of the cost function and the constraints depend only on

*n*(

_{i,j}*t*−

*d*), the solution at a given time depends only on the integrand and the constraint at that time. Therefore the optimal control problem can be replaced by the following optimization problem Assume initially that the inequality constraints are absent. Using Lagrange multipliers we get that where (

_{i,j}*f′*)

^{−1}exists due to the strictly convexity of

*f*. Because this solution obeys the inequality constraints, it serves as a solution of the full problem (including the inequality constraints). Thus the neural activity of a single neuron

*n*(

_{i.j}*t*−

*d*) is proportional to the total neural activity of the control signal

_{i,j}*u*(

_{i}*t*).

### Isometric task

Consider an isometric task, i.e., the subject retains a fixed endpoint position in the face of an external force field. Assume that the force is applied at a certain direction φ on the horizontal plane, namely According to the relation between endpoint force and joint moments, given by *Eq. 3*, we get that and therefore where *C*^{(k)}(θ) = *R*(*R*^{(k)T}*R*^{(k)})^{−1} *J*^{T}(θ) is a *M* × 3 matrix.

Under the assumption of generalized symmetry mentioned in the preceding text Hence Using trigonometry, we get that the tuning function is

### Influence of gravity

In the case where the arm is subject to gravity, some constant term is added to each torque, according to *Eq. 4*. Therefore where *rg _{i}*(θ) is a constant bias due to gravity.

As before, the tuning function is a truncated cosine. However, the width of the tuning function is not exactly 180° but wider or narrower, depending on the sign of the constant bias *rg _{i}*(θ).

## APPENDIX B

### Solution for the bijoint 6-muscles model

The kinematic equations of a bijoint planar arm (see Fig. 2*A*) are given by where *L*_{1}, *L*_{2}, θ_{1}, θ_{2} are explained in Fig. 2*A*.

The joint angles can be derived in this case solely from the hand position relative to the shoulder position, according to Using these equations, we have calculated joint angle trajectories. For the isometric task, it was assumed that the hand position is (−5 cm, 20 cm) in shoulder coordinates. For the movement task. it was assumed that the hand trajectory is a minimum jerk trajectory (Flash and Hogan 1985), i.e. where ; *x*_{0}, *y*_{0} and *x*_{1}, *y*_{1} are hand coordinates at *t* = 0 and at *t* = *t*_{f} respectively.

The static torque was calculated according to *Eq. 3*, where the Jacobian in this case is given by The dynamic torque was calculated according to *Eq. 4*, which is given in explicit form for a bijoint arm without gravity (see Fig. 2*A*) by where *m _{i}* is the mass of link

*i*,

*r*is the distance between the joint and the center of link

_{i}*i*, and

*I*is the inertia of link

_{i}*i*. We used the following parameter values, applicable to a Rhesus monkey, based on Cheng and Scott (2000), Melis et al. (2002), and Graham and Scott (2003) Given the static and the dynamic torque, we have calculated the total torque according to

*Eq. 2*.

The musculoskeletal geometry is given by six equivalent muscles, corresponding to the six classes of anatomical arm muscles: shoulder flexor and extensor, bijoint flexor and extensor, and elbow flexor and extensor (Fig. 2*B*). For this model we've used the moment arm values shown in Table 2. These values are rough averages of the moment arms reported in Graham and Scott (2003).

The arm moments matrix, under the symmetry conditions discussed in appendix a, is given by Then the MCSs were computed according to *Eq. 11.*

## GRANTS

This work was partially supported by the Ollendorff Foundation at the electrical engineering department at the Technion and by National Institute for Psychobiology (Israel) Grant to Dr. Karniel.

## Acknowledgments

We are grateful to S. Shoham for very helpful comments on the manuscript and to L. Sergio and J. Kalaska for providing details about their experimental setup. The helpful and instructive comments by the anonymous reviewers are also gratefully acknowledged.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2007 by the American Physiological Society