|Home | About | Journals | Submit | Contact Us | Français|
How the activity of populations of cortical neurons generates coordinated multi-joint actions of the arm, wrist and hand is poorly understood. This study combined multielectrode recording techniques with full arm motion capture in order to relate neural activity in primary motor cortex (MI) of macaques (macaca mulatta) to arm, wrist and hand postures during movement. We find that the firing rate of individual MI neurons is typically modulated by the kinematics of multiple joints and that small, local ensembles of MI 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 MI 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.
Coordinated actions of the upper limb may engage more than a dozen cortical areas (Kalaska, 1997, Shadmehr, 2005), with each contributing to various components of reach and grasp (Jacob, 2003). Primary motor cortex (MI) plays a central role among these areas (Lemon, 1993) and is one site where motor plans involving proximal and distal joints are likely to merge. Despite being the focus of numerous studies, MI activity has seldom been analyzed in the context of naturalistic multi-joint movements. Most of the behavioral tasks used to examine MI 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 MI. Electrical stimulation and spike triggered averaging of EMG activity show that small regions, and even single neurons in MI can facilitate the movement of both proximal and distal joints (Donoghue et al., 1992, McKiernan et al., 1998, Park et al., 2001). Recent MI recordings in humans with tetraplegia show that nearby neurons are engaged by imagined proximal and distal actions (Hochberg et al., 2006). Understanding how local MI 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 MI 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 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. In order 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 MI. 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 MI. 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).
The data presented were gathered from two male macaque 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 towards them at the end of a string. The trajectories and speeds of the objects were varied in order to elicit a variety of different reach-to-grasp movements. After holding an object for ~1 sec 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 (diameter 4.2 cm, thickness 1.1 cm) and a ring (diameter 7.5 cm, thickness 0.9 cm). The combination of objects used for each session is summarized in supplementary Table 1.
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, plc.) and synchronized through hardware with neural recordings. We analyzed a comprehensive catalogue of 25 anatomically defined upper extremity joint angles (see supplemental materials for details). The system used 12 infrared cameras operating at 240 frames per second to track the position of multiple reflective markers (4mm diameter hemispheroids) with sub-mm 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.
Neural data was recorded using microelectrode arrays (Cyberkinetics Neurotechnology Systems, Inc., CKI) chronically implanted into the MI 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 Fig. 1. Each array contained a 4 × 4 mm grid with 96 silicon-based electrodes 1mm in length and spaced 400 μm apart (see (Suner et al., 2005) for details on array structure and surgical procedures). Monkeys were head-fixed during recording. Data acquisition and storage were accomplished using a Cerebus multi-channel data acquisition system (CKI). Differences in spike waveform shape and amplitude were used to identify single unit activity using custom-made software (Vargas-Irwin, 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.
Linear state-space models (Eq. 1 and 2) 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 msec. bins successively shifted by 41.67 msec. (equivalent to 10 motion capture frames). At each time step k, a new ‘hidden state’, xk, was calculated based on the previous hidden state, xk-1, and a vector of firing rates, uk (Eq. 1). The output, zk,, representing predicted motion, is a linear function of the hidden state vector (Eq. 2). The model assumes that both the output and the state are corrupted by white Gaussian noise. Note that while, at a given time instant, the relationship between the neural input and the kinematic output is linear, this relationship changes in time due to 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, Natick, MA). While the values of the hidden states are unobserved, optimization iterates between re-estimating 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 i) predicting the state, which takes d2+dc+d operations in Eq. 1, ii) predicting the kinematic value in d+1 operations in Eq. 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 minute for each decoding variable (using a MacPro with a 2.9GHz Quad-Core Intel Xeon processor). Once model building was complete, the algorithm reconstructed individual positions for each degree of freedom in ~0.1 msec. This is considerably less time than what is required to obtain spike counts for a single 100 msec. 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 6-9 trial blocks, and each block (~ two minutes in duration) consisted of 20 to 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, while 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 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.25sec intervals we carried out approximately 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 approximately 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 over-fitting, since ‘chance’ correlations are evaluated with the data used to generate the decoding models. Overall, hidden state models improved decoding accuracy by about 20% (on average) compared to 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).
During the task, monkeys reached for and grasped objects across a volume of 3D space approximately 25 × 25 × 25 cm. This covers a considerable portion of the space the monkeys can comfortably reach (given that their arms are less than 30 cm long). Hand speeds were broadly distributed (mean of 9 cm/sec, std = 12 cm/sec, with a maximum of 93 cm/sec). 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 distal-most markers on the thumb and index), wrist angles and hand locations is shown in Fig. 2C-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 Fig. 3A (mean correlations averaged over 4 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 largely decoupled three main groups of kinematic variables: those related to the upper arm, the wrist, and the hand.
While 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 where 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, Mason et al., 2004, Todorov and Ghahramani, 2004, Ingram et al., 2008).
Our analysis of the movements elicited by dynamic reaching and grasping task show that the upper limb operates in a complex multi-dimensional 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 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 MI 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 (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 msec. 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 correlation (r) between measured firing rates and predicted firing rates. Initially all lags were set to 50 msec. In order to find the optimal combination, the lag for each kinematic variable was adjusted to maximize the value of R (computed with all variables) over 5 iterations through the full set of kinematics. In order to approach the conditions of on-line decoding, only lags where neural activity preceded movement were tested. Lags ranging from -50 to -300 msec. relative to movement (in 50 msec. 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 R2 values for both monkeys is shown in Fig. 4A. All R2 values calculated were statistically significant (p<0.01). The R2 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 (sr2) for each of these sets of kinematic variables were calculated by subtracting the values of R2 obtained with the full set of kinematics from the values of R2 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 sr2 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 sr2 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 sr2 are statistically significant (p<0.01) and greater than 0.05 (representing more than 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 non-overlapping portion of the variance).
The distribution of sr2 over all recorded neurons is shown in Fig. 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 sr2 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; about 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% (C) or 26.1% (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 Fig. 5.
Although Fig. 4E shows that approximately 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). In order to determine if there were groups of neurons that were more closely associated with either 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 (sr2) for arm joints by the sum of the sr2 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 fashion (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 zero and one hundred, 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 towards 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, more than 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 roughly 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.
Multidimensional activation profiles for MI 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 Fig. 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).
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 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 msec., (at intervals of 50 msec.), were considered. The distribution of lags used for decoding is shown in supplemental Fig. 3. The training set data was also used to generate optimized decoding models. In order 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’ fashion 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 about 30 neurons were used (Fig. 8A). Root mean squared error measures displayed a similar pattern (supplemental Fig. 4). 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) 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). In monkey C, more than 1/4 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). 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 one degree (3.72deg for C and 4.79deg for G). The distribution of errors over all joint angle reconstructions (including both monkeys) had a mean of -0.26 deg (indicative of an unbiased estimator) and standard deviation 9.52deg. The median absolute error across every monkey and joint angle was 4.13deg. This is equivalent to 6.15% when normalized with respect to the measured movement range for each DoF (see supplemental Fig. 6 for details on individual DoFs). For every session, decoding accuracy for all measured DoFs was above the 95% confidence limit for chance performance (see ‘methods’ for details on statistical testing). Fig. 9A shows a sample of measured and decoded postures from a single trial (see supplemental Movie 2 for a sample of fully reconstructed movement sequences). Fig. 9B shows 1D reconstructions of two individual kinematic parameters in detail. The Pearson correlation coefficient (r) between measured and decoded variables for each of the four datasets is summarized in Fig. 9A. The mean r across sessions over all decoded joint angles was 0.72 (standard deviation 0.094); more than 90% of the r values were above 0.6, and more than 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 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[arm vs. wrist]>0.60, p[arm vs. hand]>0.28, p[wrist vs. hand]>0.93). For monkey C, there was a trend towards more accurate decoding of arm joints (median 0.80) than wrist (median 0.7) or hand (median 0.74) joints (two-sample Kolmogorov-Smirnov test, p[arm vs. wrist]>0.032, p[arm vs. hand]>0.049, p[wrist vs. hand]>0.33).
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 MI neurons. These results demonstrate that local neighborhoods of MI neurons carry diverse information simultaneously related to both proximal and distal joints used in reaching and grasping motions of the upper limb.
We demonstrated that it is possible to reconstruct a broad range of naturalistic reach-to-grasp movements from local ensembles of MI 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 roughly comparable to what has been reported in previous studies decoding simpler movements (supplementary Table 3). At least 10 cells were needed in order to achieve near-optimal decoding accuracy for any one DoF (Fig. 8A, supplemental Fig. 4). At the same time, groups of neurons selected for hand, wrist and arm joint angle decoders formed largely 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 (even though 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) and could generalize well across different types of movements. For example, decoded grip aperture values could be consistently decoded even though six to nine different objects (engaging apertures between 10 and 70mm) were grasped during each session (Fig. 1 and supplemental Fig. 8). Based on these results we can conclude that local neuronal populations in the MI 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 about 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 multi-dimensional 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 MI neurons. Aflalo and Graziano (Aflalo and Graziano, 2006, Aflalo and Graziano, 2007) measured 8 DoFs of the upper limb while recording MI neural activity using single electrodes. This study also used semipartial correlations in order to evaluate the relationship between neural activity and individual DOFs. They concluded that neurons were, on average, significantly related to 4 out of the 8 parameters measured, suggesting that MI neurons simultaneously represent various aspects of limb motion, with different neurons emphasizing different combinations. Jackson and colleagues (Jackson et al., 2007) reached a similar conclusion based on the relationship between neural activity and electromyographic (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 MI employs a distributed control scheme where multiple DoFs are represented simultaneously even at the level of individual neurons, often including both proximal and distal joints (Fig. 4, Fig. 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 MI neurons: any correlate of joint angles may be the ‘true’ coded variable. Our results cannot conclusively determine how MI represents or transforms information for reach and grasp actions, but they do mean that theories of MI 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 multi-joint kinematics. However, simple anatomical explanations cannot account for the combinations of anatomically disparate, weakly correlated parameters that we observe.
It is well established that the output of MI neurons diverges to multiple muscles (Shinoda et al., 1981, Buys et al., 1986, Fetz, 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, 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). While these findings support the concept of a holistic control scheme where individual neurons influence many DoFs, no prior study has explored how MI 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 MI, 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 MI, 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 MI neurons with direct cortico-motoneuronal connections to distal musculature Rathelot and Strick, 2009) could display higher selectively to particular upper extremity muscles.
Among all motor areas, MI has the most direct influence on alpha 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. Even though 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 3 or 4 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, since it could allow direct and simultaneous control of arm posture (elbow and shoulder), hand orientation in space (wrist) and grip aperture from a local MI 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 MI population providing reach and grasp information. Further experimentation will be required to determine if it is possible to decode the more complex movements of the fingers independently in order 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 in order 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), further studies will be necessary to determine if 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 on-line 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 regimes. If MI 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 neuron's full activation profile. 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 multi-joint decoding.
While 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 MI 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.
The authors 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. This research project was supported by the National Institute of Health (NIH/NINDS F32 NS061483, NIH/NINDS NS25074, RO1EB007401-01) the National Center for Research Resources (NCRR C06-16549-01A1), National Science Foundation (IIS-0534858), the Office of Naval Research (N00014-07-1-0803, award N00014-16-1-0692) the Katie Samson Foundation and the Office of Research and Development, Rehabilitation R&D Service, Department of Veterans Affairs.
Conflict of interest disclosure for John P. Donoghue: prior Chief Scientific Officer and director, stock holdings, Cyberkinetics Inc.