Decoding Complete Reach and Grasp Actions from Local Primary Motor Cortex Populations (original) (raw)

Articles, Behavioral/Systems/Cognitive

, Gregory Shakhnarovich, Payman Yadollahpour, John M. K. Mislow, Michael J. Black and John P. Donoghue

Journal of Neuroscience 21 July 2010, 30 (29) 9659-9669; https://doi.org/10.1523/JNEUROSCI.5443-09.2010

Abstract

How the activity of populations of cortical neurons generates coordinated multijoint actions of the arm, wrist, and hand is poorly understood. This study combined multielectrode recording techniques with full arm motion capture to relate neural activity in primary motor cortex (M1) of macaques (Macaca mulatta) to arm, wrist, and hand postures during movement. We find that the firing rate of individual M1 neurons is typically modulated by the kinematics of multiple joints and that small, local ensembles of M1 neurons contain sufficient information to reconstruct 25 measured joint angles (representing an estimated 10 functionally independent degrees of freedom). Beyond showing that the spiking patterns of local M1 ensembles represent a rich set of naturalistic movements involving the entire upper limb, the results also suggest that achieving high-dimensional reach and grasp actions with neuroprosthetic devices may be possible using small intracortical arrays like those already being tested in human pilot clinical trials.

Introduction

Coordinated actions of the upper limb may engage more than a dozen cortical areas (Kalaska et al., 1997; Shadmehr and Wise, 2005), with each contributing to various components of reach and grasp (Jacob and Jeannerod, 2003). Primary motor cortex (M1) plays a central role among these areas (Lemon, 1993) and is one site in which motor plans involving proximal and distal joints are likely to merge. Despite being the focus of numerous studies, M1 activity has seldom been analyzed in the context of naturalistic multijoint movements. Most of the behavioral tasks used to examine M1 activity require the subjects to operate in a restricted movement space (such as a horizontal plane) and/or use averaged neural responses over multiple repetitions of similar movements. These constraints tend to remove trial-specific details and result in artificially high correlations between multiple joint angles, making it difficult to determine the relationship between movement variables and neuronal firing.

There is abundant evidence that neurons related to distal and proximal actions intermingle considerably within M1. Electrical stimulation and spike-triggered averaging of electromyographic (EMG) activity show that small regions and even single neurons in M1 can facilitate the movement of both proximal and distal joints (Donoghue et al., 1992; McKiernan et al., 1998; Park et al., 2001). Recent M1 recordings in humans with tetraplegia show that nearby neurons are engaged by imagined proximal and distal actions (Hochberg et al., 2006). Understanding how local M1 populations control highly flexible coordinated limb movements not only has important implications for understanding volitional movement control but also for the design of neuroprosthetic devices that attempt to reproduce reach and grasp actions from neural activity.

The goal of this study was to determine whether local neuronal populations in M1 encode information sufficient to reconstruct the joint angles of the arm, wrist, and hand across a broad range of naturalistic reach and grasp actions aimed at various objects moving through a three-dimensional (3D) workspace. Such actions elicit correlations among joints that approach those occurring during natural movements, which are much weaker than those observed during repetitive, stereotyped behaviors typically used in experimental settings. To relate motion of the entire upper limb to neural activity, we combined marker-based motion capture techniques capable of measuring 25 individual joint angles with cortical recordings from 96-microelectrode arrays chronically implanted in the arm/hand region of M1. We demonstrate the presence of information related to the entire limb by reconstructing the full set of 25 joint angles during continuous reach and grasp movements using population activity, showing that it is possible to extract information related to the large number of degrees of freedom (DoFs) engaged in naturalistic movements from a small region of M1. Our findings further suggest that current microelectrode array technology, identical to that now being evaluated in human pilot clinical trials, could potentially be used to control a realistic robotic arm and hand or even reanimate multiple muscles in a paralyzed limb using extant functional electrical stimulation techniques (Moritz et al., 2008), going far beyond the few dimensions of neural control that have been achieved in able-bodied monkeys (Carmena et al., 2003; Velliste et al., 2008) or in humans with paralysis (Hochberg et al., 2006; Kim et al., 2008).

Materials and Methods

Task.

The data presented were gathered from two male macaque monkeys (monkeys C and G) over four different experimental sessions (two for each monkey). The behavioral task they performed was designed to elicit a broad range of reach and grasp movements that would differentially engage multiple joints of the hand and arm. Monkeys were trained to single-handedly intercept and hold objects swinging toward them at the end of a string. The trajectories and speeds of the objects were varied to elicit a variety of different reach-to-grasp movements. After holding an object for ∼1 s, they received a fruit juice reward. The objects used in the grasping task (Fig. 1) were made of plastic or wood and included balls (2 and 7 cm diameter), cylinders (0.8 cm diameter × 12 cm long and 3.2 cm diameter × 15 cm long), a cube (faces 3.5 cm), a rectangular prism (9 × 1.1 × 1.1 cm), an isosceles triangular prism (8.8 × 2.3 × 3.4 cm), a disk (4.2 cm diameter, 1.1 cm thickness), and a ring (7.5 cm diameter, 0.9 cm thickness). The combination of objects used for each session is summarized in supplemental Table 1 (available at www.jneurosci.org as supplemental material).

Figure 1.

The monkeys used different grasping strategies characterized by varied patterns of grip aperture scaling and wrist motion to intercept and hold each of the objects. Monkeys were trained to intercept and hold objects of various shapes and sizes swinging toward them at the end of a string. Measurements of grip aperture, wrist pronation/supination, and wrist flexion/extension are shown for nine segments of data in which different objects were grasped. The objects are shown above the plots, and consistent scaling was used to preserve relative dimensions. All measurements were taken during a single session (C2).

Motion capture.

Upper limb motion from shoulder to fingers was recorded using reflective markers tracked with an optical motion capture system (Vicon Motion Systems; Oxford Metrics Group) and synchronized through hardware with neural recordings. We analyzed a comprehensive catalog of 25 anatomically defined upper extremity joint angles (for details, see supplemental material, available at www.jneurosci.org). The system used 12 infrared cameras operating at 240 frames/s to track the position of multiple reflective markers (4-mm-diameter hemispheroids) with submillimeter accuracy. A total of 29 markers were attached to each monkey using mild water-soluble adhesive (Fig. 2a). A kinematic model of the arm and hand (Fig. 2b) was specified and fit to observed marker data, by least-squares error minimization using a commercially available software package (Vicon iQ). The model was rooted at the shoulder, extending to the distal interphalangeal joints. Parameters of the articulated model (segment lengths, marker offsets) were calibrated in the beginning of each session. Values for each DoF were derived from an articulated model fit to the measured positions of markers attached to the arm and hand. In particular, Euler angles were used to represent relative joint rotations. The representation for shoulder rotations was converted to polar coordinates (azimuth, elevation, and internal/external rotation). Supplemental Movie 1 shows a sample of the captured arm and hand movements (available at www.jneurosci.org as supplemental material).

Figure 2.

The dynamic grasping task elicits a wide range of upper limb movements. a, Movement was recorded by tracking 29 reflective markers attached to the monkey with mild water-soluble adhesive. b, Model of the hand and arm fit to the 3D position of the markers to calculate the joint angles for each frame. c, Distribution of measured grip apertures (top) and wrist angles (bottom). d, Projection of the position of the arm endpoint (proximal wrist marker) onto the coronal, sagittal, and transverse planes (dataset C1). Color saturation indicates the density of the points. Marginal distributions of density are shown along each axis. Crosses denote the mean and quartiles. The highest density of points (highest color saturation) was noted at the location where the monkeys typically rested their arm after each reach and grasp action.

Neural recording.

Neural data were recorded using microelectrode arrays (Cyberkinetics Neurotechnology Systems) chronically implanted into the M1 upper limb area of two macaque monkeys (just anterior to the central sulcus at the level of the genu of the arcuate sulcus). The position of each array on the cortical surface is shown in supplemental Figure 1 (available at www.jneurosci.org as supplemental material). Each array contained a 4 × 4 mm grid with 96 silicon-based electrodes 1 mm in length and spaced 400 μm apart (for details on array structure and surgical procedures, see Suner et al., 2005). Monkeys were head fixed during recording. Data acquisition and storage were accomplished using a Cerebus multichannel data acquisition system (Cyberkinetics Neurotechnology Systems). Differences in spike waveform shape and amplitude were used to identify single-unit activity using custom-made software (Vargas-Irwin and Donoghue, 2007). We recorded 122 and 88 neurons in two sessions in monkey C and 35 and 30 neurons in two sessions in monkey G. Within-animal recordings may have included overlapping populations. The degree of uniqueness of each population is impossible to establish definitively. Consequently, we treated each dataset as an independent population for analysis.

Decoding 25-dimensional movement.

Linear state-space models were used to reconstruct the time-varying value of each DoF using neural ensemble firing rates (for an introduction to state space models, see Ljung, 1999). Firing rates were calculated in partially overlapping 100 ms bins successively shifted by 41.67 ms (equivalent to 10 motion capture frames). At each time step k, a new “hidden state,” xk, was calculated based on the previous hidden state, x_k_−1, and a vector of firing rates, uk: Embedded Image The matrix F represents the evolution of the hidden state vector x based on the previous state, L represents the influence of firing rates u on the state vector x, and wk represents Gaussian zero-mean noise in the state dynamics. The output, zk, representing predicted motion, is a linear function of the hidden state vector (Eq. 2) corrupted by zero-mean Gaussian noise (Vk): Embedded Image Note that although, at a given time instant, the relationship between the neural input and the kinematic output is linear, this relationship changes in time as a result of the evolution of the hidden state. The size of the hidden state vectors was empirically set to 3. The model was trained (optimized) using an iterative algorithm implemented in the system identification toolbox of MatLab (MathWorks). Although the values of the hidden states are unobserved, optimization of the matrices F, L, and C was accomplished by iterating between reestimating the most likely values of the hidden states and minimizing, via gradient descent, the output prediction error under these values.

The computational complexity of the decoding algorithm can be described as follows. Suppose the hidden state x has dimension d, the dimension of neural signal u is c, and the kinematic value predicted is a scalar (i.e., we decode a single DoF). After the model has been trained, the computations for each time step require (1) predicting the state, which takes _d_2 + dc + d operations in Equation 1, and (2) predicting the kinematic value in d + 1 operations in Equation 2. Thus, for n time steps, the complexity of decoding is O(nd(d + c + 2)), which is linear in the number of time steps as well as in the dimension of the neural signal. For the datasets described here, model building using training data was performed in under 1 min for each decoding variable (using a MacPro with a 2.9 GHz Quad-Core Intel Xeon processor). Once model building was complete, the algorithm reconstructed individual positions for each degree of freedom in ∼0.1 ms. This is considerably less time than what is required to obtain spike counts for a single 100 ms time bin, demonstrating that the algorithm is fast enough for real-time online decoding.

All decoding was performed offline on test data not used for training the decoder. Although all kinematic variables and neural activity were recorded simultaneously, each DoF was decoded individually. Each session consisted of six to nine trial blocks, and each block (∼2 min in duration) consisted of 20–40 reaching and grasping actions using a single object. Each trial block was subdivided into training, validation, and testing segments. The training segment spanned the first 60% of the data for each block, whereas the validation and testing segments each included 20%. The training and validation sets were used to generate and optimize the decoding models. The final decoding results were obtained by running the optimized models on the “test” dataset (20% of the data collected for each object). The accuracy of the decoding results was evaluated using Pearson's correlations (r) between measured and reconstructed kinematics. For each parameter, we evaluated the statistical significance of the decoding results using a Monte Carlo approach. We calculated the correlation between kinematic variables reconstructed using neural test data and sets of kinematics of equivalent length taken from the training and validation segments. This strategy estimates the “chance” performance by comparing decoded movements with kinematics collected from the same monkey at a different point in the session. By shifting the training/validation data sampling window by 0.25 s intervals, we performed ∼2000 comparisons for each dataset (with all comparisons made between data segments collected the same day from the same monkey). Overall, this procedure allowed us to sample ∼8000 chance correlations for each kinematic variable. We used the empirical 95% confidence limit over these distributions of values as the bound for chance performance. Note that this procedure relies on the large variation in movement parameters present in our data. Applying this method to highly stereotyped repetitive movements, chance levels of performance would easily match decoding results. Additionally, this evaluation method guards against overfitting, because chance correlations are evaluated with the data used to generate the decoding models. Overall, hidden state models improved decoding accuracy by ∼20% (on average) compared with a simple least-squares linear filter (results not shown). All decoders were constructed using 0th-order kinematics (for example, joint angles). Explicitly decoding velocities and accelerations using the same model did not improve performance (results not shown).

Results

Kinematic properties of the dynamic reaching and grasping task

During the task, monkeys reached for and grasped objects across a volume of 3D space ∼25 × 25 × 25 cm. This covers a considerable portion of the space the monkeys can comfortably reach (given that their arms are <30 cm long). Hand speeds were broadly distributed (mean ± SD of 9 ± 12 cm/s, with a maximum of 93 cm/s). The presentation of objects of various sizes and shapes swinging in different directions both close to and away from the monkey elicited a broad range of reaches and grasps. The distribution of grip apertures (the distance between the distalmost markers on the thumb and index), wrist angles, and hand locations is shown in Figure 2, c and d. Different objects elicited varying patterns of grip aperture and wrist motion associated with a range of possible one-handed reaching and grasping strategies (Fig. 1).

Although the behavioral task produced a rich variety of different reaches and grasps, the inherent natural correlations across joint angles were not completely dissociated. Kinematic correlations across a representative set of joints are shown in Figure 3a (mean correlations averaged over four sessions). Stronger correlations were most evident among actions about the same joint (such as the wrist or shoulder). Digit joint angles also tended to be highly correlated with each other because the monkeys nearly always grasped objects while moving all fingers in concert (in Fig. 3a, these variables are summarized using grip aperture). Despite these correlations, the task primarily decoupled three main groups of kinematic variables: those related to the upper arm, the wrist, and the hand.

Figure 3.

The dynamic grasping task essentially decouples the DoFs of the hand, wrist, and arm. a, Correlations display the relative coupling between DoFs, which are weaker as distance between joints increases. Values for eight representative parameters averaged across monkeys are shown. s., Shoulder; e., elbow; w., wrist; int./ext. rot., internal/external rotation; pron./sup., pronation/supination; uln./rad., ulnar/radial deviation; flex./ext., flexion/extension. b, Dimensionality estimation for reach and grasp actions. This plot shows the cumulative percentage of variance accounted for as a function of the number of PCs used to represent the joint angles for each of four sessions. c, The firing rate of individual M1 neurons was only moderately correlated with any given joint angle. Correlation coefficients (CC) between the firing rate of each neuron and each measured joint angle are divided into three groups: those related to the arm (shoulder and elbow, 4 DoFs), the wrist (3 DoFs), and the hand (18 DoFs).

Although we measured 25 joint angles, their correlations demonstrate they are not all used independently (as is to be expected in a linked and constrained system). We used principal component analysis (PCA) to estimate the number of linearly uncorrelated DoFs that were actually engaged in the task. PCA generates a linear transform that maps the space of measured joint angles onto a new coordinate system in which axes are aligned with the directions of maximum angular variation. Our results show that an average of 10 principal components (PCs) accounted for >95% of the variance across all datasets (Fig. 3b), suggesting that the recorded actions could be closely reproduced by a 10-dimensional control signal. Previous studies using PCA to estimate the underlying dimensionality of arm and hand movements in both monkeys and humans have reported both lower and higher values, depending on the number of DoFs measured and the relative complexity of the behavior analyzed (Mason et al., 2001, 2004; Todorov and Ghahramani, 2004; Ingram et al., 2008).

Neuron tuning properties

Our analysis of the movements elicited by dynamic reaching and grasping task show that the upper limb operates in a complex multidimensional space. Our next goal was to determine what types of motion were most closely related to changes in the activity of individual neurons. We first evaluated the relationship between neural activity and motor behavior by calculating Pearson's correlations (r) between single-neuron firing rates and individual kinematic variables. Rather than being highly correlated with small subsets of tightly coupled joint angles, the firing rates of M1 neurons tended to be moderately correlated with multiple relatively independent joint angles (often including both proximal and distal joints). Of these correlation coefficients, 95% had an absolute value <0.38 and 50% were <0.16 (Fig. 3c).

Correlations between a neuron and multiple kinematic parameters will necessarily reflect the correlation among the kinematics themselves to some degree. Multiple regression analysis using the conventions described by Cohen and Cohen (1983) was used to determine the extent to which neural correlations could be accounted for by kinematic correlations. In the first step of the analysis, coefficients of multiple correlation (R) between the full set of joint angles (and their velocities) and the firing rates of individual neurons (in 100 ms bins) were computed as follows: linear regression was used to obtain a least-squares prediction of the firing rate for each neuron given the full set of joint angles. R values were then obtained by computing the Pearson's correlation (r) between measured firing rates and predicted firing rates. Initially all lags were set to 50 ms. To find the optimal combination, the lag for each kinematic variable was adjusted to maximize the value of R (computed with all variables) over five iterations through the full set of kinematics. To approach the conditions of online decoding, only lags in which neural activity preceded movement were tested. Lags ranging from −50 to −300 ms relative to movement (in 50 ms steps) were considered. For each neuron, only the regression results with optimized lags were analyzed further. Squaring the value of R yields an estimate of the variance in firing rate explained by the kinematic parameters. The distributions of _R_2 values for both monkeys is shown in Figure 4a. All _R_2 values calculated were statistically significant (p < 0.01). The _R_2 values obtained in the first step of the analysis were then used to calculate semipartial correlations between firing rates and three sets of joint angles: those associated with the arm (shoulder and elbow), the wrist, and the hand. Squared semipartial correlation coefficients (_sr_2) for each of these sets of kinematic variables were calculated by subtracting the values of _R_2 obtained with the full set of kinematics from the values of _R_2 obtained with a reduced set excluding the variables of interest. This procedure excludes the fraction of variance redundantly accounted for by multiple sets of kinematics. Values of _sr_2 can therefore be interpreted as the amount of firing rate variance accounted for by one set of kinematics (arm, wrist, or hand) beyond what is accounted for by the other two. In this sense, the values of _sr_2 reflect the variance in firing rates exclusively related to one set of kinematic variables. We shall refer to neurons as having a set-specific relationship (i.e., with a specific subset of kinematic variables) when the values of _sr_2 are statistically significant (_p_ < 0.01) and >0.05 (representing >5% of the variance in firing rate). Note that a cell can have multiple set-specific relationships if the variance in firing rates can be most completely explained by linear regression when multiple sets of kinematics are used (with each set explaining a nonoverlapping portion of the variance).

Figure 4.

Neural correlations with multiple sets of kinematic parameters cannot be accounted for by correlations among the kinematics themselves. In all histograms, results are shown in solid gray for monkey C and in black outlines for monkey G. a, Coefficients of multiple correlation (_R_2) for each of the recorded neurons. These values can be interpreted as the fraction of variance in firing rates accounted for by linear relationships with the full set of kinematics (joint angles and joint angle velocities). Note that in all cases neural activity preceded movement, because the analysis was limited to predictive lags to approximate online decoding conditions. b–d, Semipartial correlation coefficients (sr2) for three sets of kinematic variables. These values can be interpreted as the amount of firing rate variance accounted for by one set of kinematic variables (arm, wrist, or hand) beyond that accounted for by the other two. Thus, b shows the variance in neuronal firing rates explained exclusively by arm kinematics after removing the variance explained by hand and wrist variables. c and d show similar distributions for wrist and hand joint angles. e, Overlap between groups of neurons that show statistically significant semipartial correlations >0.5 (representing >5% of the variance in firing rate) to arm, wrist, or hand kinematics. Values for monkeys C and G are shown separately in each partition. Overall, 86.7% of the neurons recorded in monkey C and 75.4% of the neurons recorded in monkey G displayed a set-specific relationship with at least one set group of kinematic variables. Note that the areas in the Venn diagram are not to scale. MAD, Median absolute deviation.

The distribution of _sr_2 over all recorded neurons is shown in Figure 4b–d. For monkey C, 84.2% of the neurons had a set-specific relationship to hand kinematics, demonstrating that neural correlations to hand joints exceeded what could be accounted for by kinematic correlations for most of the neurons examined. A relatively smaller fraction of neurons, 46.6%, displayed a set-specific relationship with arm kinematics, and only 3.8% of the neurons displayed a set-specific relationship with wrist kinematics. Similar results were obtained for monkey G, although the fraction of neurons with a significant relationship to the wrist were not as markedly reduced (with 62.9, 39.9, and 12.1% for the hand, arm, and wrist respectively). The groups of neurons with significant _sr_2 values for arm and hand joints overlapped considerably (Fig. 4e), indicating that many neurons were involved in controlling movement about multiple joints, often including proximal and distal; approximately one-third of the neurons in each monkey showed a set-specific relationship for both hand and arm kinematics (44.7% in monkey C and 30.7% in monkey G). It is important to note that these overlaps, by definition, cannot be accounted for by the linear correlations in the kinematic variables. Other combinations comprised relatively small fractions of the neuronal populations with one exception: neurons with a set-specific relationship only to the hand accounted for either 37.6% (monkey C) or 26.1% (monkey G) of the total number of neurons. Raster plots showing the movement-related activity of various kinds of neurons classified according to their semipartial correlations are shown in Figure 5.

Figure 5.

Single-unit activity in M1 during naturalistic reaching and grasping movements. Raster plot showing the action potentials fired by 25 individual neurons over a span of 5.5 s. Raster plots are colored according to set specific relationships with kinematics (significant semipartial correlations >0.05; for details, see Results, Neuron tuning properties): blue, hand only; green, wrist only; red, hand only; orange, hand + wrist; magenta, hand + arm; black, hand + arm + wrist. The data shown corresponds to session G2. Of the 30 neurons recorded in M1, only the 25 shown had set-specific relationships to at least one set of kinematic parameters. During the time span shown, the monkey performed three separate reach-to-grasp movements targeting a small ball. Grip aperture measures are overlaid over the raster plot to highlight the moments when prehension occurs (troughs). The full posture of the arm is shown for four time points in the first movement.

Although Figure 4e shows that ∼40% of the sampled neurons have significant semipartial correlations to multiple sets of kinematic parameters, it does not show whether these neurons emphasize one particular set over others (that is to say, if one particular set tends to contribute to most of the variance). To determine whether there were groups of neurons that were more closely associated with arm, wrist, or hand kinematics, we defined a descriptive statistic referred to as the index of kinematic selectivity. The index of kinematic selectivity for arm joints was calculated for each neuron by dividing the squared semipartial correlation (_sr_2) for arm joints by the sum of the _sr_2 values for arm, wrist, and hand joints and expressing the result as a percentage (Fig. 6a). A value of zero for arm selectivity means that whatever variance in firing rate can be explained by arm kinematics is redundantly explained by either wrist or hand kinematics. A value of 100% means that all set-specific variance (the sum of all semipartial correlations) is accounted for by arm kinematics, meaning that, aside from variance explained by arm kinematics, no additional variance is explained by the hand or wrist. Indices of kinematic selectivity for the wrist and hand were computed in a similar manner (Fig. 6b,c). As expected from the semipartial correlations, wrist kinematics tended to contribute a relatively small portion of independent variance in firing rates. Most of the variance could be attributed to hand and arm joint angles. If neurons mainly reflected either proximal or distal kinematics, we would expect to see a bimodal distribution of hand and arm kinematic selectivity indices with peaks close to 0 and 100, indicating that only one set of parameters accounted for most of the variance of a given neuron. The observed pattern did not conform to this hypothesis: there was not a sharp division between neurons more tightly coupled to either hand or arm. There was a bias toward higher kinematic selectivity for the hand, with median values of 60.40 and 46.83 for monkeys C and G, respectively (compared with 28.35 and 32.85% for arm selectivity). However, >65% of neurons (66.19 in monkey C and 73.85 in monkey G) displayed hand selectivity values between 30 and 70. Likewise, between 46.19% (monkey C) and 60% (monkey G) of the neurons analyzed displayed arm selectivity indexes between 30 and 70. This pattern suggests that the properties of individual neurons span a broad spectrum of combinations between proximal and distal kinematics, often emphasizing both to a generally comparable extent. These findings contradict the hypothesis that neurons are directly correlated to a single kinematic parameter and that any correlation to others is merely an incidental reflection of the relationship among the kinematics themselves. Our results support the alternative hypothesis that neural encoding operates in high dimensional spaces and that these spaces can include combinations of proximal and distal kinematic parameters.

Figure 6.

Most neurons were not preferentially related to specific subsets of kinematic parameters, even after removing the influence of correlations between kinematics. In all panels, results are shown in solid gray for monkey C and in black outlines for monkey G. a, Index of kinematic selectivity for arm joints (for details, see Results). These values are a ratio of the variance in firing rates exclusively represented by arm joints (semipartial correlation for arm kinematics) over the total set-specific variance (the sum of semipartial correlations for the arm, wrist, and hand). b, c, Index of kinematic selectivity for wrist and hand joints, respectively. MAD, Median absolute deviation.

Multidimensional activation profiles for M1 neurons can illustrate the relationship between firing and multiple behavioral variables (akin to a tuning function or receptive field). Although it is not possible to visualize the full 25-dimensional joint movement space, it is possible to display multiple 2D slices representing the relationship between pairs of variables and neuronal spiking. Examples of these activation profiles for neurons simultaneously tuned to grip aperture, shoulder, and wrist kinematics are shown in Figure 7. The intermingling of kinematic representations is also evident in the spatial organization of neurons, which can be derived from the regular 10 × 10 organization of the electrode array. Neurons with similar primary features (highest joint angle correlates) showed no clear spatial pattern of aggregation across the regularly sampled region of the electrode array (supplemental Fig. 2, available at www.jneurosci.org as supplemental material).

Figure 7.

M1 neurons are simultaneously tuned to proximal and distal DoFs. Pairs of measured DoFs were discretized into 100 equally sized 2D bins spanning 99% of the data (after removing outliers at each extreme). Colored circles show the average firing rate for each bin (color was interpolated for the space between the circles). Each row shows the activity of a single M1 neuron. The first column shows firing rate as a function of shoulder elevation and grip aperture. The second and third columns have similar displays, using grip aperture versus wrist pronation/supination and shoulder elevation versus wrist pronation/supination.

Arm and hand movement reconstruction

Population decoding methods provide a means to test whether information is contained in the activity of a set of neurons. Here we tested the hypothesis that small neuronal populations contained sufficient information to reconstruct entire reach and grasp actions. Specifically, we used linear state space models (see Materials and Methods) to reconstruct all 25 measured joint angles based on neural firing rates. Four additional parameters, including grip aperture and the x, y, and z position of the arm endpoint, were also decoded in an identical way. The training set (60% of the data) was used to select the lag maximizing the linear correlation coefficient (r) between the firing rate of each cell and each DoF. A set of lags from 50 to 300 ms (at intervals of 50 ms) were considered. The distribution of lags used for decoding is shown in supplemental Figure 3 (available at www.jneurosci.org as supplemental material). The training set data was also used to generate optimized decoding models. To avoid overfitting, the subset of cells selected to reconstruct each DoF was chosen by optimizing decoding accuracy on the validation set (20% of the data) in the following procedure. First, a decoder was built using the single cell with the highest r value (raw firing rate vs measured DoF). Next, a new decoder was built using the two cells with the highest r value (the same cell used in the first decoder plus a new one). This process continued in a “greedy” manner by adding one cell at a time and rebuilding the decoder. The mean r value (measured vs decoded) over all variables tended to peak or asymptote when ∼30 neurons were used (Fig. 8a). Root mean squared error measures displayed a similar pattern (supplemental Fig. 4, available at www.jneurosci.org as supplemental material). We therefore used 30 neurons per DoF to build the final set of decoders used to reconstruct the kinematics from the testing data segment (not used in parameter optimization steps). All decoding results shown [except for Fig. 8a and supplemental Fig. 4 (available at www.jneurosci.org as supplemental material)] were taken from this testing segment. Note that this method may yield a different set of cells in the decoder for each DoF. However, neuron subgroups chosen to decode individual DoFs tended to overlap extensively, even when a large number of neurons (80–100) were available (supplemental Fig. 5, available at www.jneurosci.org as supplemental material). In monkey C, more than one-fourth of the neurons recorded in each session were simultaneously used in decoders for at least one DoF from the hand, wrist, and upper arm (supplemental Table 2, available at www.jneurosci.org as supplemental material). Because 30 or 35 units (session 1 or session 2) were recorded from monkey G, nearly all were used to decode each DoF. Despite the limited number of neurons for monkey G, decoding accuracy was similar between the two animals (Fig. 8c). Although the magnitude of the angular errors was significantly higher for monkey G (two-sample Kolmogorov–Smirnov test, p < 1 × 10−6), the difference between the median absolute error measures for both monkeys was only slightly more than 1° (3.72° for monkey C and 4.79° for monkey G). The distribution of errors over all joint angle reconstructions (including both monkeys) had a mean ± SD of −0.26 ± 9.52° (indicative of an unbiased estimator). The median absolute error across every monkey and joint angle was 4.13°. This is equivalent to 6.15% when normalized with respect to the measured movement range for each DoF (for details on individual DoFs, see supplemental Fig. 6, available at www.jneurosci.org as supplemental material). For every session, decoding accuracy for all measured DoFs was above the 95% confidence limit for chance performance (for details on statistical testing, see Materials and Methods). Figure 9_a_ shows a sample of measured and decoded postures from a single trial (for a sample of fully reconstructed movement sequences, see supplemental Movie 2, available at www.jneurosci.org as supplemental material). Figure 9_b_ shows 1D reconstructions of two individual kinematic parameters in detail. The Pearson's correlation coefficient (_r_) between measured and decoded variables for each of the four datasets is summarized in Figure 9_a_. The mean ± SD _r_ across sessions over all decoded joint angles was 0.72 ± 0.094; >90% of the r values were above 0.6, and >60% were above 0.7. Although decoding accuracy was significantly higher for monkey C (two-sample Kolmogorov–Smirnov test, p < 0.02), the difference between the distributions of _r_ values was small, with means of 0.74 and 0.70 for monkeys C and G, respectively. The _r_ values obtained for hand, wrist, and arm joints (with medians of 0.72, 0.73, and 0.73, respectively) were not significantly different in monkey G (two-sample Kolmogorov–Smirnov test, _p_ > 0.60 for arm vs wrist, p > 0.28 for arm vs hand, p > 0.93 for wrist vs hand). For monkey C, there was a trend toward more accurate decoding of arm joints (median of 0.80) than wrist (median of 0.7) or hand (median of 0.74) joints (two-sample Kolmogorov–Smirnov test, p > 0.032 for arm vs wrist, p > 0.049 for arm vs hand, p > 0.33 for wrist vs hand).

Figure 8.

Optimal decoding was accomplished using populations of 30 neurons to reconstruct each DoF. a, Mean r (over all variables) between measured and decoded kinematics as a function of the number of neurons used to decode each DoF. Neurons were added to the models in order from most to least correlated with the DoF being decoded (following the “greedy”' selection method described in Results, Arm and hand movement reconstruction). Based on this graph, the final set of decoding models were constructed using 30 neurons for each DoF. The total number of neurons used for decoding in each session (pooled across all DoFs) is shown next to each plot as a fraction of the total number of neurons recorded. For details on the overlap of populations used to decode arm, wrist, and hand movements, see supplemental Figure 7 (available at www.jneurosci.org as supplemental material). b, Error distribution across all joint angles reconstructed across all test datasets using sets of 30 neurons (solid gray histogram for monkey C and outlined black histogram for monkey G). Triangles along the _x_-axis denote the range that encompasses 95% of the data for each distribution (green for monkey G and red for monkey C). More than 50% of the reconstructed joint angles were within 5° of the measured values. MAD, Median absolute deviation.

Figure 9.

Decoding continuous 25-dimensional movement. a, Examples of measured (ghost) and decoded (solid) arm postures from a reach and grasp trial (with each of the 25 joint angles decoded independently). b, Detailed view of measured (blue) and reconstructed (black) values for grip aperture and shoulder azimuth. c, Correlation coefficients between measured and decoded variables. Colored dots represent the values for each experimental session, and solid bars mark the mean over all sessions. Black asterisks represent chance levels of performance (95% confidence limit for correlations between reconstructed kinematics and temporally shifted kinematics; for details, see Results). In addition to joint angles, grip aperture, as well as the x, y, and z position of the endpoint of the arm were directly decoded. Decoding accuracy was above chance for every degree of freedom examined. MAE, Mean absolute error; In./Ex. Rot., internal/external rotation; Flex./Ext., flexion/extension; Ul./Rad., ulnar/radial deviation; Pron./Sup., pronation/supination; MCP, metacarpophalangeal; Ante./Retro., anteposition/retroposition; Rad. Ab./Ad., radial abduction/adduction; Palm. Ab./Ad., palmar abduction/adduction; PIP, proximal interphalangeal.

The lack of significant differences in decoding accuracy across joints shows that the neuronal population in monkey G contained similar amounts of information about all measured DoFs. Although movement reconstructions in monkey C tended to be more accurate for arm joint angles, wrist and hand kinematics were also decoded above chance levels. Our findings show that it is possible to reconstruct complex naturalistic multidimensional reaching and grasping movements using the firing rates of relatively small populations of M1 neurons. These results demonstrate that local neighborhoods of M1 neurons carry diverse information simultaneously related to both proximal and distal joints used in reaching and grasping motions of the upper limb.

Discussion

We demonstrated that it is possible to reconstruct a broad range of naturalistic reach-to-grasp movements from local ensembles of M1 arm/hand area neurons using standard decoding algorithms. In every session, the decoding accuracy of all 25 measured DoFs exceeded chance performance. Decoding accuracy of each DoF was approximately comparable with what has been reported in previous studies decoding simpler movements (supplemental Table 3, available at www.jneurosci.org as supplemental material). At least 10 cells were needed to achieve near-optimal decoding accuracy for any one DoF (Fig. 8a) (supplemental Fig. 4, available at www.jneurosci.org as supplemental material). At the same time, groups of neurons selected for hand, wrist, and arm joint angle decoders formed primarily overlapping sets (Fig. 8b), and small populations of as few as 30 neurons were capable of providing information related to the entire upper limb. The fact that decoding results were very similar between the two animals (although the number of neurons recorded in monkey G was much smaller) suggests that neuron selection is not critical. Our decoding models faithfully replicated the complex high-dimensional characteristics observed in the behavioral task (supplemental Fig. 7, available at www.jneurosci.org as supplemental material) and could generalize well across different types of movements. For example, decoded grip aperture values could be consistently decoded, although six to nine different objects (engaging apertures between 10 and 70 mm) were grasped during each session (Fig. 1) (supplemental Fig. 8, available at www.jneurosci.org as supplemental material). Based on these results, we can conclude that local neuronal populations in the M1 arm/hand area contain information related to a large set of non-redundant DoFs of the upper limb. Our examination of neural activity supports and expands this hypothesis by showing that even individual neurons are capable of representing unrelated kinematic parameters such as shoulder elevation and grip aperture (which are not only anatomically distant but are also shown to be only weakly correlated in the data we collected). Multiple regression analysis showed that approximately one-third of the neurons in each monkey displayed a direct relationship with both proximal and distal joint angles exceeding what can be accounted for based solely on kinematic correlations (Fig. 4). Classifying these cells as “hand” or “arm” neurons would ignore the richness of their tuning properties. Instead of pigeonholing neurons into fixed classes, we developed the index of kinematic selectivity to represent the strength of the relationship between firing rates and kinematics along a continuous scale. Rather than clustering at extreme values (as would be expected if proximal and distal movements were represented separately at the level of individual neurons), the range of selectivity indices observed for arm and hand kinematics suggests that the properties of individual neurons span a broad spectrum of combinations between proximal and distal movement parameters (Fig. 6). The rich variety of information represented at the level of single neurons would help explain why relatively small ensembles are sufficient to perform high-dimensional movement reconstructions.

Our results may seem at odds with previous studies that identify neurons with strong correlations to particular movement parameters. However, most of these studies operate in restricted movement spaces, which tend to enhance correlations across kinematic variables. This type of experimental paradigm could effectively collapse complex multidimensional tuning functions into simpler forms, resulting in artificially higher correlations between kinematics and neural activity.

The results of previous studies analyzing complex movements representative of the natural behavioral repertoire of primates support our conclusions about the information content of individual M1 neurons. Aflalo and Graziano (2006, 2007) measured eight DoFs of the upper limb while recording M1 neural activity using single electrodes. This study also used semipartial correlations to evaluate the relationship between neural activity and individual DOFs. They concluded that neurons were, on average, significantly related to four of the eight parameters measured, suggesting that M1 neurons simultaneously represent various aspects of limb motion, with different neurons emphasizing different combinations. Jackson et al. (2007) reached a similar conclusion based on the relationship between neural activity and EMG data obtained from the forearm, suggesting that “precise information about activation of a particular arm muscle can be readily obtained from most cells in the arm area of motor cortex.” In agreement with these findings, our results suggest that M1 uses a distributed control scheme in which multiple DoFs are represented simultaneously even at the level of individual neurons, often including both proximal and distal joints (Figs. 4, 7). Our findings expand on previous work by quantifying the independent contribution of different groups of joints to single-neuron activity and demonstrating that information reflecting high-dimensional movement of the entire upper limb can be extracted from simultaneously recorded local neuronal ensembles. Simultaneous multielectrode recordings also allowed us to evaluate all neurons under precisely the same set of behavioral conditions, which was not the case for previous single-electrode studies.

Although we showed that joint angles could be directly decoded from neural activity, this does not mean that these variables are directly encoded by M1 neurons: any correlate of joint angles may be the “true” coded variable. Our results cannot conclusively determine how M1 represents or transforms information for reach and grasp actions, but they do mean that theories of M1 function must incorporate the presence of integrated reach and grasp information at the single-cell and local population level. Our data does not include EMG or other records of movement dynamics. It is possible that synergies relating the activity of some individual neurons to many joints could be explained by a representation of multi-articulate muscles instead of multijoint kinematics. However, simple anatomical explanations cannot account for the combinations of anatomically disparate, weakly correlated parameters that we observe.

Implications for cortical organization

It is well established that the output of M1 neurons diverges to multiple muscles (Shinoda et al., 1981; Buys et al., 1986; Fetz et al., 1989), sometimes including groups acting on both proximal and distal joints (McKiernan et al., 1998). Conversely, the input to a given muscle appears to converge from a wide area of cortex (typically several square millimeters) that overlaps extensively with the neuronal pools of other muscles (Landgren et al., 1962; Andersen et al., 1975; Jankowska et al., 1975; Sato and Tanji, 1989; Donoghue et al., 1992; Schieber and Hibbard, 1993; Schieber, 1996; Sanes and Schieber, 2001; Smith and Fetz, 2009). Although these findings support the concept of a holistic control scheme in which individual neurons influence many DoFs, no previous study has explored how M1 ensembles relate to multidimensional naturalistic movements engaging the entire limb. The fact that we can extract information related to the arm, wrist, and hand from a relatively small area of cortex is consistent with a highly integrated representation that controls reach and grasp as a functional collective. Although the area spanned by our microelectrode arrays does not constitute a comprehensive survey of M1, our findings suggest that at least some small regions are capable of representing entire limb actions. Our recordings were limited to a 4 × 4 mm area on the surface of M1, and it is possible that the region located deeper within the bank of the central sulcus [which has been shown to contain most of the M1 neurons with direct cortico-motoneuronal connections to distal musculature (Rathelot and Strick, 2009)] could display higher selectively to particular upper extremity muscles.

Implications for neuroprosthetics

Among all motor areas, M1 has the most direct influence on α motor neuron pools driving upper limb muscles (Dum and Strick, 1991) and is therefore a promising location for the detection of complex multidimensional control signals for neuroprosthetic devices. Our findings suggest that signals from small regions of motor cortex of persons with tetraplegia could potentially provide control signals, allowing high-dimensional control of robotic limbs or functional electrical stimulation systems capable of reanimating multiple paralyzed muscles. Previous studies have accomplished robotic limb control but with many dimension-reducing constraints that limit the range and flexibility of possible motions. Hochberg et al. (2006) showed that a person with tetraplegia was able to operate a 1D proportional control robotic hand and could use a 2D interface to perform reach and grasp with a simple robot arm. Velliste et al. (2008) demonstrated neural control of a robotic arm by able-bodied monkeys. This was accomplished by decoding velocity information for a 3D endpoint and converting it to 4D joint angle coordinates by setting a fixed elbow swivel angle. Although the arm was equipped with a gripper under proportional control, the task only required binary control (open or closed). In a similar experiment, Carmena et al. (2003) demonstrated explicit proportional control of grip force, although arm endpoint trajectories were restricted to a 2D plane. In both of these studies, correlations induced by the constraints of the task mean that the intrinsic dimensionality may have been lower than the three or four decoded DoFs. In comparison, we successfully measured and reconstructed 25 joint angles. Although correlations among joints reduced the intrinsic DoFs necessary to describe the task to ∼10, this is substantially more than addressed by any previous study. This marked dimensionality increase in the control signal is of great relevance for neural prosthetic systems, because it could allow direct and simultaneous control of arm posture (elbow and shoulder), hand orientation in space (wrist), and grip aperture from a local M1 population using signals that may be used in the CNS to accomplish the same function. Because our task did not elicit significant individuated finger movements, our results cannot address the ability to reconstruct them from the same M1 population providing reach and grasp information. Additional experimentation will be required to determine whether it is possible to decode the more complex movements of the fingers independently to achieve direct neural control of an even higher number of functional DoFs.

Significant challenges still remain before high-dimensional neuroprosthetic technology can be fully developed for clinical applications. Although we restricted our decoding analysis to rely exclusively on neural activity preceding movement to minimize the influence of proprioceptive inputs, the results presented here were obtained in primates with intact somatic sensation. Although low-dimensional cursor control has been demonstrated for paralyzed humans with limited somatosensory feedback (Hochberg et al., 2006; Kim et al., 2008), additional studies will be necessary to determine whether high-dimensional control can be achieved under similar conditions. Despite the difficulties posed by implementing real-time closed-loop systems, it is also possible that the presence of natural or artificial feedback could improve the accuracy of high-dimensional neuroprosthetic control through online error correction. However, means to correlate the desired movement of multiple joints with neural activity related to intended or imagined actions for those who cannot actually move (and may have compromised sensory feedback) is not obvious. Our data suggests that taking into account the integrated nature of neural representations may be important for the development of neuroprosthetic training regimens. If M1 neurons do indeed operate in high-dimensional spaces that integrate diverse aspects of movement kinematics, characterizing the activity of a given neuron while one joint moves and the rest remain fixed would only yield a one-dimensional projection of the full activation profile of a neuron. A proper characterization of neural activity would require a wide sampling of possible kinematics, which would be impossible to achieve one joint at a time. Therefore, training data including naturalistic movement such as the ones elicited by our behavioral task may be necessary for effective multijoint decoding.

Although our decoding strategy produced accurate results for a wide range of complex actions, improved control might be obtained with more advanced mathematical methods. In our current decoding scheme, each DoF is reconstructed independently without taking into account interactions between multiple DoFs. However, our observations suggest that the activity of M1 neurons reflects these interactions, potentially providing an even greater source of information for decoding algorithms to exploit. It is possible to modify our state space model to represent multiple kinematic parameters simultaneously. However, as the dimensionality of the decoding space grows, the amount of training data required to provide a representative sample increases exponentially (a problem commonly referred to as the “curse of dimensionality”). When dealing with a limited set of training data, this generally results in decreased decoding accuracy. Algorithms that efficiently represent the integrated motion of the limb using dimensionality reduction techniques could potentially circumvent this problem and improve decoding accuracy.

Footnotes

This research project was supported by National Institutes of Health (NIH)/National Institute of Neurological Disorders and Stroke Grants F32 NS061483 and NS25074, NIH Grant RO1EB007401-01, National Center for Research Resources Grant NCRR C06-16549-01A1, National Science Foundation Grant IIS-0534858, Office of Naval Research Grant N00014-07-1-0803 and Award N00014-16-1-0692, the Katie Samson Foundation, and the Office of Research and Development, Rehabilitation Research and Development Service, Department of Veterans Affairs. We thank Beth Travers, Allan Rydberg, and Loretta Reiss for their assistance with animal care and instrumentation design, as well as Lachlan Franquemont, Kaivon Paroo, and Jessica Feldman for their help during data collection.
Correspondence should be addressed to Carlos E. Vargas-Irwin,Department of Neuroscience, Brown University, Box G-LN, Sidney E. Frank Hall for Life Sciences, 185 Meeting Street, Providence, RI 02912. Carlos_Vargas_Irwin{at}brown.edu