## Abstract

We used statistical methods for spherical density estimation to evaluate the distribution of preferred directions of motor cortical cells recorded from monkeys making reaching movements in 3D space. We found that this distribution, although broad enough to represent the entire 3D continuum of reaching directions, exhibited an enrichment for reaching forward from the body and, to a lesser degree, for reaching backward toward the body. The distribution of preferred directions of cells in the motor cortex may have important implications for motor cortical function and for the decoding of arm trajectories from population activity.

## INTRODUCTION

Directional tuning in reaching is a prominent feature of motor cortical activity (Caminiti et al. 1990; Georgopoulos et al. 1982, 1986; Schwartz 1994). The tuning of a motor cortical cell is characterized by its preferred direction (PD), which is a unit vector that specifies the direction of reaching at which the cell's discharge rate peaks. Although PDs have been calculated for thousands of cells, the detailed distribution of PDs (DPD) across the directional continuum has not yet been investigated.

A rigorous characterization of this distribution could influence our understanding of how the directional signal generated by populations of motor cortical neurons is decoded, or “read-out.” Such knowledge is important in the context of understanding how natural reaching behaviors are controlled by the motor system and in the context of achieving control of neuroprosthetic devices (Schwartz 2004; Taylor et al. 2002). The structure of the DPD may also have important implications for our understanding of motor cortical plasticity, as discussed in the following text.

It is important to investigate the DPD in the context of unconstrained reaching movements in three-dimensional (3D) space. Such movements occur naturally and most frequently in the everyday life of primates and thus it is reasonable to suppose that the neural mechanisms of reaching would relate to them. In contrast, two-dimensional (2D) movements (Georgopoulos et al. 1982) can only partially relate to those neural mechanisms because they are, by definition, constrained. Furthermore, correct inferences concerning the DPD can be obtained only from the 3D case. Although the basic functional property of directional tuning is present for both 3D and 2D movements, a particular PD obtained from 2D movements represents a projection into the movement plane of the cell's true, 3D PD; therefore the 2D PD is consistent with an infinite number of 3D PDs (Amirikian and Georgopoulos 2003). In this work, we considered only PDs that were estimated in the context of a 3D reaching task.

Current knowledge about the DPD for unconstrained, 3D movements derives from three studies (Caminiti et al. 1990; Mitsuda and Onorati 2002; Schwartz et al. 1988). In the first investigation of the DPD (Schwartz et al. 1988), it was observed that a sample of roughly 500 PDs extracted from two monkeys was dispersed across the 3D directional continuum. It was also observed that the DPD exhibited an enrichment for reaches made forward, away from the body; however, when applied to this sample, a Rayleigh test for uniformity against a unimodal alternative did not reject the null hypothesis of uniformity. In the study of Caminiti et al. (1990), a sample of nearly 200 PDs extracted from three monkeys again revealed that the DPD is sufficiently broad to represent all portions of the 3D continuum, although the null hypothesis of uniformity was not explicitly tested against any other alternatives. Most recently, it was reported in Mitsuda and Onorati (2002) that the hypothesis of uniformity was rejected by a Rayleigh test (*P* < 0.05) for a sample of 126 PDs obtained from a single monkey, although explicit measures of the location and strength of the implied bias in the DPD were not provided.

In this study we analyzed the structure of the DPD using a variety of statistical tests and methods for density estimation. We applied these methods to a sample that included the original data described in Schwartz et al. (1988), as well as a more recent data set that includes roughly 800 PDs extracted from widely distributed recording sites in the motor cortex of two monkeys (Naselaris et al. 2005). Our findings are consistent with previous observations that PDs are distributed broadly across the 3D directional continuum. However, we found that the DPD is not *strictly* uniform, as the result of an enrichment of the representation for forward reaching directions (as observed originally; Schwartz et al. 1988), as well as for backward reaches directed toward the body. This enrichment induces a relative increase in accuracy of the motor-cortical representation for forward and backward reaches, as measured by the neuronal population vector; however, when averaged across all directions, the effect on accuracy is quite small. Our results are consistent with previous demonstrations (Georgopoulos et al. 1986, 1988; Schwartz 1994) that the distribution of PDs in motor cortex provides an excellent substrate for a population-level code of reaching direction.

## METHODS

### Behavioral task

The neural data used in these studies were collected from four rhesus monkeys (NI, IO, A, and B) engaged in the 3D center-out task (Georgopoulos et al. 1986). Light-cued reaches were made from a fixed, central starting position toward eight targets located near the corners of a cube. Data from all four monkeys refer to reaches of the left arm. The starting position was level with the shoulder of the monkey's reaching arm. The monkey was instructed to hold the starting position for a variable amount of time (0.5–1.5 s), after which a button located on a randomly selected target was illuminated. The monkeys reached toward the target and depressed the button for 0.08 s to obtain a juice reward. For each recording site, monkeys performed 40 reaches toward the targets (five blocks of eight reaches; targets were presented in random order within each block).

### Neural recordings

All recordings were obtained from the right hemisphere of the motor cortex. A detailed account of the recording protocols for monkeys NI and IO was previously given elsewhere (Schwartz et al. 1988). For monkeys A and B, recording sites were located at various depths beneath a patch of cortex that extended 3–4 mm along the central sulcus and 7–12 mm in the direction perpendicular to the central sulcus. This region was centered about 15 mm from the midline along the medial–lateral axis. Microstimulation at 5–20 μA within the boundaries of this region evoked contractions of distal and proximal arm muscles. Recording sites along each electrode penetration were spaced 150 μm apart. At each site, raw extracellular membrane potentials were sampled at 60 kHz. Single-unit activity was extracted from these records using the Plexon off-line sorter (Plexon, Dallas, TX).

### Calculation of PDs

PDs were calculated using a multiple linear regression analysis (Georgopoulos et al. 1986). For each reach, an average firing rate was obtained for the interval spanning from 0 to 60 ms after the onset of the target to the end of the movement. The average firing rate *r* was treated as the dependent variable in a regression relating the direction of reach to the firing rate where *d*_{x}, *d*_{y}, and *d*_{z} are the direction cosines of the target and *b*_{0}, *b*_{x}, *b*_{y}, and *b*_{z} are regression coefficients. The components of the cell's preferred direction (PD) vector **p** are given by where All cells included in these analyses demonstrated a significant directional tuning effect (*P* < 0.05), as indicated by a bootstrapping procedure (Lurito et al. 1991).

### Nonparametric density estimation

We used a kernel-based method to obtain a nonparametric estimate of the DPD as the probability density *g*(**x**), of reaching direction **x**, given a set of *n* experimentally determined PDs where **p**_{i} is the PD of cell *i*. The kernel function *f* is the von Mises–Fisher (vMF) probability density function (Fisher et al. 1987) In this context, κ > 0 acts as a smoothing parameter.

The value of κ was determined using a maximum likelihood with cross-validation approach suggested in Fisher et al. (1987). In this approach, we estimate the density at **p**_{j} using a sample of *n* − 1 PDs created by omitting **p**_{j} We then maximized the log-likelihood function *L*(κ) = ∑_{i} log [*g*_{i}(**p**_{i})], with respect to κ over all *n* PDs in the sample.

### Parametric density estimation

The parametric model of the DPD was specified as where α_{i} ∈ **α** is a mixing coefficient, κ_{i} ∈ **κ** is the concentration parameter for the *i*th component of the model, and *m* = 2 is the number of components. The density associated with the *i*th component peaks at **μ**_{i} ∈ **M**, the mode location. The size of the peak at **μ**_{i} is determined by κ_{i} > 0 and α_{i} > 0; peak size increases as κ_{i} and/or α_{i} increases. The mixing coefficients satisfy ∑ α_{i} = 1 and are interpreted as the probability that a given PD belongs to component *i*. **α**, **κ**, and **M** were estimated using the expectation maximization (EM) procedure (Banerjee et al. 2003). To avoid solutions corresponding to local maxima, the EM algorithm was run 100 times from different, randomly chosen initial conditions (0 < κ < 1 for each initial condition). Solutions obtained from these randomly varying initial conditions were highly consistent. The solution that generated the highest likelihood under the model is reported in Table 1.

### Test for significance of density peaks in the DPD

We used a bootstrapping procedure to test the significance of the deviation from uniformity of individual peaks in the nonparametric distribution. In this procedure, 10^{3} samples of *n* PDs were drawn from a uniform distribution on the sphere. For each such sample, a nonparametric density estimate was constructed using the procedure described above and the maximum density values were stored. The resulting distribution was used to test the significance of local peaks in the density values obtained from the nonparametric model of the DPD.

### Decoding of simulated neural activity

The firing rate *r*(**x**) for a simulated cell when a reach is made in direction **x** was defined as *r*(**x**) ≡ *N*(**x**)/Δ*t*. The temporal interval Δ*t* was set to 700 ms, corresponding to the average trial duration. The spike count *N*(**x**) was drawn from a Poisson distribution with mean given by (*b*_{0} + *k***p**^{T}**x**)Δ*t*, where *b*_{0} and *k* were drawn at random from the population of real motor-cortical cells and **p** was drawn either from this same population, a uniform distribution on the sphere, or a two-component mixture of vMF distributions. Movement direction was decoded from a population of such cells using the neuronal population vector (NPV; Georgopoulos et al. 1986). The normalized firing rate for each cell, in the form (*r* − *b*_{0})/*k*, was used as the weighting factor for calculating the NPV. Decoding error was measured as the angle between the actual reaching direction and the direction estimated by the decoding algorithm.

## RESULTS

We analyzed a sample of 1,139 PDs recorded from four monkeys performing a 3D center-out task. We will refer to this as the *combined sample*. Recordings were obtained from sites across the arm area of the primary motor cortex, including sites in the anterior bank of the central sulcus, and the exposed surface of the precentral gyrus. An equal-area projection plot of the PDs in the combined sample (Fig. 1*D*; see Fig. 1, *A*–*C* for an explanation of the projected coordinate system) illustrates the two most salient features of the DPD. First, PDs cover the entire 3D directional continuum. Second, there is a local increase in PD density for reaches that are almost straight forward (φ ≈ 180°) and slightly above the equator (θ < 90°) and another, smaller increase in PD density for reaches that are nearly backward (φ ≈ 360°) and below the equator (θ > 90°).

We estimated PD density directly by constructing a nonparametric model of the DPD (Figs. 1, *E* and *F* and 2). The nonparametric model revealed local peaks in PD density for forward and backward reaches; it also confirmed that the magnitude of the backward peak was clearly smaller than that of the forward peak. A two-component mixture model was then used to capture locations of the peaks (black dots, Figs. 1*E* and 2; see Table 1). The ratio of PD densities at the forward and backward peak locations identified by the mixture model is 1.3, indicating that the enrichment of the representation for backward reaches is indeed weaker than that for forward reaches.

To determine whether the forward and backward peaks represent true departures from uniformity, we tested their significance using a bootstrap procedure. This test revealed that both the forward and backward peaks are significantly (*P* < 10^{−3}) larger than the random fluctuations in density expected from a uniform DPD (left and right arrows in Fig. 1*E* show expected density and bound at which *P* = 10^{−3}, respectively). Consistent with this result, the Bingham test of uniformity against a bimodal alternative rejected the null hypothesis of uniformity (*P* < 0.01; see Table 1).

To test the consistency of the properties of the DPD, we applied the same statistical tests and density estimation procedures to two data sets separately (Fig. 2). Enrichments of forward and backward reaching directions are evident in both data sets. The results of the tests of significance for the separate data sets, shown in Table 1, leave little doubt that the properties of the DPD are robust.

The nonuniformity of the DPD could affect the accuracy with which movement direction can be decoded from neuronal population activity, enhancing decoding accuracy for the forward and backward reaching directions, while introducing a global bias that increases the average level of decoding error. A widely used method for decoding arm directions is the neuronal population vector (NPV) (Caminiti et al. 1991; Georgopoulos et al. 1986, 1988; Schwartz 1994; Steinberg et al. 2002; Taylor et al. 2002). The NPV is a weighted sum of PDs and is thus sensitive to first-order features of the underlying PD distribution (Georgopoulos et al. 1988; Mussa-Ivaldi 1988; Salinas and Abbott 1994; Sanger 1994; Scott et al. 2001; Steinberg et al. 2002). To assess the effect of the DPD on global decoding accuracy, we calculated population vectors from the activity of simulated cells whose PDs were drawn at random from our experimental sample. The asymptotic level of error, averaged over 100 uniformly and randomly distributed reaching directions, is roughly 12° (Fig. 3*A*). Thus the statistically significant bias in the DPD has a fairly small effect on average decoding accuracy. However, the error of the NPV is smallest for reaching directions that are close (or opposite) to the forward peak in the DPD (Fig. 3*B*). Furthermore, the length of the NPV—which is an indication of its reliability as an estimate of reaching direction—is greatest for reaching directions near (or opposite to) the forward peak. Thus the quality of the directional command produced by the NPV is greatest for reaching directions near the forward/backward peaks in the DPD (Fig. 3*C*).

As shown in Fig. 3*B*, there is a local minimum in the error for reaches at exactly 90 ° to the forward peak. For reaches in this direction, cells with preferred directions near the forward peak do not fire at rates significantly different from baseline and thus make no contribution to the population vector. As a consequence, the population vector length reaches a global minimum for these directions.

## DISCUSSION

Our findings are consistent with previous observations that preferred directions in the motor cortex are sufficiently broadly dispersed to cover the full continuum of 3D reaching directions (Caminiti et al. 1990; Schwartz et al. 1988). We have also shown that embedded within this broad dispersal of PDs is an enrichment for reaches directed forward (Schwartz et al. 1988) and an enrichment for reaches directed backward. These enrichments are manifested as local peaks in the DPD, and represent a significant departure from strict uniformity. The structure of the DPD permits a highly accurate representation of the direction of reach, as measured by the application of the NPV. The accuracy of this representation is greatest for forward and backward reaching directions.

What is the explanation for the enrichment of the motor cortical representation for forward and backward reaches? One possible explanation has to do with the motor cortex's capacity for massive, behaviorally dependent reorganization (Sanes and Donoghue 2000). Motor-cortical plasticity has been demonstrated using a variety of techniques, including intracortical microstimulation (Nudo et al. 1996), slice preparations (Hess and Donoghue 1994), functional magnetic resonance imaging (Karni et al. 1995), and transcranial magnetic stimulation (TCMS) (Classen et al. 1998; Pascual-Leone et al. 1995). Collectively, work on motor-cortical plasticity has demonstrated that the cortical representation of movements and movement sequences associated with the learning and production of a stereotyped motor task are enlarged as a result of the learning. This well-replicated finding has been interpreted as evidence that “motor practice induces the recruitment of additional M1 units into a network specifically representing the trained motor sequence” (Karni et al. 1995).

One study that supports this interpretation is particularly important in the context of our findings. In the study of Classen et al. (1998), focal TCMS of the motor cortex of humans consistently evoked thumb movements in a single direction. Subjects were then instructed to repeatedly produce movements of the thumb in the opposite of the evoked direction for 30 min. Subsequently, the direction of the TCMS-evoked thumb movement changed to match the “trained” direction. This finding suggests that the repetition of a movement in a single direction can bias the motor-cortical representation of direction toward the repeated one.

Given these results on motor-cortical plasticity, we conjecture that the enrichment of the motor cortical representation for forward and backward reaching directions may simply reflect an increase in the incidence of reaches in these directions, relative to other directions, in the everyday life of the monkey. A direct test of this conjecture would require a characterization of the statistics of the natural reaching behavior of the monkey and a thorough sampling of PDs taken before and after intentional manipulation of these statistics. To our knowledge, such an experiment has not yet been attempted.

Whether the enrichment for specific reaching directions in the DPD has any functional effect on arm movement depends on the manner in which motor cortical commands are “read out,” or translated into action by circuits in the spinal cord. If read-out is implemented by a weighted sum of cortical activities (an NPV-like mechanism), the nonuniformity of the DPD can induce a relative increase in accuracy (hyperacuity) for forward and backward reaches (see Fig. 3*C*). However, if read-out is implemented by an optimal method (Salinas and Abbott 1994), the direction-dependent errors will be eliminated. Our results thus provide a basis for testing whether an NPV-like read-out mechanism is implemented by the motor system. Although some studies have shown directionally dependent variations in reaching accuracy (Lamotte and Acuna 1978; Smyrnis et al. 2000), no measurements of reaching accuracy have been made for 3D movements performed by monkeys with a known DPD.

A final important issue related to the observed enrichment of forward and backward directions in the DPD is the spatial scale at which this enrichment exists. PDs in this study were extracted from cells dispersed across a large (about 30 mm^{2}) cortical region. The large sampling area underlying our database thus admits of at least two possibilities for the spatial organization of PDs. One possibility is that spatially separate motor-cortical regions may generate the distinct peaks apparent in the DPD. Alternatively, the structure of the DPD may be replicated locally in all regions of the motor cortex. We address this issue directly in the companion paper.

## 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 © 2006 by the American Physiological Society