Subthalamic Neurons Encode Both Single- and Multi-Limb Movements in Parkinson’s Disease Patients

The subthalamic nucleus (STN) is the main target for neurosurgical treatment of motor signs of Parkinson’s disease (PD). Despite the therapeutic effect on both upper and lower extremities, its role in motor control and coordination and its changes in Parkinson’s disease are not fully clear. We intraoperatively recorded single unit activity in ten patients with PD who performed repetitive feet or hand movements while undergoing implantation of a deep brain stimulator. We found both distinct and overlapping representations of upper and lower extremity movement kinematics in subthalamic units and observed evidence for re-routing to a multi-limb representation that participates in limb coordination. The well-known subthalamic somatotopy showed a large overlap of feet and hand representations in the PD patients. This overlap and excessive amounts of kinematics or coordination units may reflect pathophysiology or compensatory mechanisms. Our findings thus explain, at the single neuron level, the important subthalamic role in motor control and coordination and indicate the effect of PD on the neuronal representation of movement.


Results
We recorded the activity of 89 units in the STN of ten Parkinson's disease patients intra-operatively (Table 1).

Direct Neuronal Encoding of Kinematics.
To directly evaluate the neuronal code, we modeled the firing rate function by linear models based solely on kinematics: orientation, angular velocity, and acceleration (see Methods). Figure 1 displays the firing rates of three example neurons along with their linear estimation models. For a neuron to be considered "directly encoding kinematics", we required two criteria: the correlation between its model and firing rate is significant (see Methods) and the model explains a considerable percentage of the variability of the firing rate (coefficient of determination: r 2 > 0. 30). Based on these two criteria, 43% (38/89) of the recorded units directly encoded movement kinematics. This percentage is comparable to the one reported in the literature for passive movements (49%) 20 .
The time lags between the neuronal firing and movement kinematics were almost uniformly distributed in the examined range between (− 1)s and 1s (Fig. 2a). The optimal time lags of models generated for each type of movement (unipedal/bipedal/unimanual/bimanual) also expand the whole range ( Fig. 2(b-i)). For all types of movement, there were units whose firing preceded the movement itself, whereas for others, it succeeded the movement. Similarly, when examining the different types of kinematic models (i.e., models based on orientation, angular velocity, acceleration or the combination of the three), time lags showed large variations. As discussed below, this lack of specificity may be due to Parkinson's disease.
Neurons Containing Kinematic Information. Do subthalamic neurons also utilize other, indirect, maybe non-linear, encoding schemes to represent kinematic information? We investigated this question using the entropy correlation coefficient, an information-based measure of dependence 21 , and found significant relations between firing rate and kinematics in 93% (83/89) of recorded units in at least one condition (Bonferroni corrected). All 83 responsive units were significantly related to all three kinematic parameters: orientation, angular velocity, and acceleration in a manner that was independent of the other two parameters (conditional entropy correlation coefficient) in at least one condition. This indicates that almost all of the recorded units encoded kinematic information in their firing rate.

Motor Coordination.
Bimanual coordination dysfunction is a sign of Parkinson's disease 22,23 . We therefore first compared whether patients performed movements of two limbs and of a single limb at a similar pace. The duration of the two types of movement (coordinated vs. single-limb) performed by the same limb did not differ significantly (paired-sample t-test, p > 0.06 in all 12 tests: all paces (slow, normal, fast), examined for each of the 4 limbs).
Most of the recorded units (78%; 69/89) were related to motor control of both limbs, where movement of the other (non-encoded) limb either activated the neuron (i.e., units responsive during bipedal or bimanual movements, but not during unipedal or unimanual ones), or deactivated it (i.e., units responsive during unipedal or unimanual movements, but not during bipedal or bimanual ones). We refer to the former as "pure bipedal" or "pure bimanual" units, and to the latter as "pure unipedal" or "pure unimanual" units. An example of a pure bipedal neuron appears in Fig. 1b, where significant relationships between the firing rate and foot orientation can be established only during bipedal movements, but not during unipedal ones. Figure 1c exhibits a pure unimanual neuron, whose firing and kinematics are strongly correlated during the unimanual movements, but not during the bimanual ones. Only a minority of 16% (14/89) of the recorded units were not affected by the other limb. They represented kinematics of the same limb independently of whether or not the other limb participated in the movement (i.e., during both uni-and bi-pedal movements or during both uni-and bi-manual ones; see Fig. 1a P1  M  45  R  15  3  19*  27*  60*   P2  M  41  R  8  2  Not Available   P3  M  70  R  17  3  3  0  15  6  19  7   P4  M  66  R  9  2.5  15  8  33  26  50  36   P5  M  61  R  6  3  9  1  15  3  28  8   P6  M  59  R  8  2  17  14  39  23  59  40   P7  M  57  R  10  2  10  2  17  15  27   Each graph shows the smoothed firing rate (blue) and the corresponding linear regression kinematic model (red) at a certain condition, along with the coefficient of determination (r 2 ) between them (see Methods). In parentheses above each graph are the kinematic parameters composing the model, and the pace of movements or performing limb employed during this condition. (a) The neuron is responsive to left foot movements during both the unipedal and bipedal conditions. It is also related to orientation of the right foot during bipedal movements, but this may be due to the high correlation between orientation of the left and right feet. The correlation is much lower for unipedal right foot movements. (b) Bipedal control neuron. The firing rate is related to normal-pace orientation of either foot during bipedal movements, but not during unipedal movements. (c) Unimanual control neuron. A combined model of orientation, angular velocity, and acceleration during slow movements explains high percentages of the variability in firing rate during unimanual, but not bimanual, movements.
in this population was significant during bipedal or bimanual tapping, but not during any of the two unipedal or two unimanual conditions. As these findings suggest, in a subpopulation of 7% (6/89) of the recorded units, each unit was related to both bipedal and bimanual movements, but not to any movement of one limb by itself, indicating involvement in left-right limb coordination in both upper and lower extremities. The proportion of this population agreed with random mixture of the pure bipedal with the pure bimanual properties (see Methods; Chi-square test of homogeneity, p = 0.17, χ 2 = 1.90, 1 degree of freedom). More than half the recorded neurons (51%) encoded the kinematics of single-limb movements only, but not during bipedal or bimanual movements. Table 2 presents the segmentation of this population according to the laterality of the related limb: contralateral, ipsilateral or units related to both ipsi-and contra-lateral single-limb movements, and to the type of extremity: upper, lower, or related to single-limb movements of both upper and lower extremities. The percentage of units in the latter subpopulation (Table 2, row 3) matches a random mixture of encoding of feet and hand movements (Chi-square test of homogeneity, p = 0.53, χ 2 = 0.39, 1 degree of freedom; see Methods). Thus, our recordings do not lend support to a significant subthalamic abstraction of whether the acting limb is upper or lower. In contrast, the percentage of units related to both contra-and ipsi-lateral movements (Table 2, Column 3) was significantly higher than would be expected by a random mixture of properties in the population (Chi-square test of homogeneity, p = 0.04, χ 2 = 4.44, 1 degree of freedom; see Methods). Our findings thus indicate that subthalamic units encode movement kinematics independently of the laterality of the performing limb, and similarly for single-limb hand movements. Some units represented 1-foot movements (independently of the acting foot), and similarly, some units encoded 1-hand kinematics. Localization of Neurons Representing Kinematics and Somatotopy. To localize sub-areas of the STN with high concentration of units directly encoding kinematics, we examined for each electrode the percentage of units it recorded which directly encode movement kinematics (Fig. 3a). We found that the more inferior electrodes, at least 3.6 mm below the AC-PC line, recorded significantly higher percentages of kinematics-encoding neurons (mean = 64%, SE = 11%) in comparison with those in the superior part of the recorded area (less than 3.6 mm below the AC-PC line; mean = 15%, SE = 9%; two-sample right-tailed t-test: p = 0.0027). Note, that due to operating room constraints which prevented an exhaustive search, the most inferior electrode sampled was 5.4 mm below the AC-PC line.
Over a quarter of the recorded units (27%; 24/89) held kinematic information related to feet movements only (i.e., not to hand movements), whereas only 4% (4/89) were solely related to hand movements (i.e., not feet movements). The proportions of the two populations were significantly different (Chi-square test of homogeneity, p < 0.0001, χ 2 = 17.0, 1 degree of freedom). These populations were organized somatotopically within the superior 40% of the STN, with feet movements at the ankle represented more superior-mesially, and hand movements at the wrist, more inferior-laterally (Fig. 3b). Notwithstanding, almost two thirds (62%; 55/89) of the recorded units were related to both feet movements and hand movements. This large overlap in representation was also located somatotopically (Fig. 3b).

Discussion
We identified units in the STN of patients with Parkinson's disease representing movement kinematics in a complex manner that depends on the limb performing the movement. The vast majority of recorded units were affected by the movement, or lack of movement, of both limbs. We identified a population of STN units related to the control of bipedal and bimanual movements which was not involved in controlling the corresponding unipedal or unimanual movements. A complementary subthalamic population participated in the control of unipedal or unimanual movements, but not of bipedal or bimanual ones. We observed somatotopic organization of the STN, with the ankle represented more superior-mesially, and the wrist, inferior-laterally, and with a large overlap between the two.
In the neurologically-intact STN of the monkey, the responses of 91% of the responding cells were elicited by manipulations around a single joint only 24 . That study indicates a high degree of limb specificity in the responsiveness of normal STN neurons. When comparing this with our findings of a large neuronal population responsive to multi-limb movements, assuming similarity between humans and monkeys, we infer that the excessive amount of multi-limb-related units may be a result of either pathophysiology due to Parkinson's disease, or a compensatory mechanism that re-routes motor coordination function to the STN. The decrease in specificity of the neuronal response in the STN may also be related to the decrease in accuracy of movements of patients with Parkinson's disease (i.e., higher variability of movement endpoints) 25 .
In patients with Parkinson's disease, previous work by Theodosopoulos et al. demonstrated that the majority of the movement-related units responded to passive movements of a single joint or multiple joints of the same extremity, with only less than 5% of the units responsive to movement of multiple limbs 20 . This percentage perfectly matches our finding of 4.5% of the units related to both feet and hand movements, each performed during single-limb movement (Table 2, 3 rd data row, last column). Our results thus extend these findings from passive, non-rhythmic movements 20 to also include voluntary, rhythmic active ones. More importantly, our research expands the type of movements to include ipsilateral movements and two-limb coordinated movements, actions not studied by Theodosopoulos et al. 20 .
The present results demonstrate that in Parkinson's disease the STN combines information about the acting limbs with kinematic information about their movement. We found a population of units whose encoding of single-limb movements is apparently independent of the laterality of the performing limb: feet movements, independent of the actual performing foot, or hand movements, independent of the actual hand. The high proportion of this population within the recorded neuronal population cannot be explained by a random mixture of properties, indicating a representational abstraction of laterality in the STN. However, it seems that the STN may not "bind together" movements of the feet and hand of the same body side, as the proportion of neurons encoding this type of movements does match a random mixture of properties in the population. The proportion of the population of units that are both pure bipedal and pure bimanual also matches random mixture of these properties, lending further support to the lack of inherent "binding" of upper and lower extremities in the STN. The separate  populations of pure bimanual and pure unimanual units explain why stimulation of the STN has a different effect on bimanual and unimanual movements 11 . Our finding of a somatotopic organization of the STN expands the reported somatotopy in the monkey 24,26,27 and human 28,29 to the human wrist and ankle joints. We describe, in addition, a large sub-area with an overlap between feet and hand representations, where the same unit participates in control of both types of movement. This overlap zone may explain why no clear somatotopic map of upper-and lower-extremity receptive fields could be discerned in the sagittal plane 30 .
The large overlap may reflect changes in the somatotopy due to Parkinson's disease. A similar finding has been reported in the monkey: neurons in the GPe/GPi and pallidonigral thalamus of the parkinsonian monkey lost their specificity to single contralateral joint and single preferred direction, and responded to multiple joints of the upper or lower extremities 31,32 . The reduced neuronal specificity also matches the responsiveness to multiple kinematic parameters: acceleration, angular velocity and orientation, in our results. Similarly, the high percentage Higher percentages of such units were recorded by electrodes inferior to the yellow dashed line, which is 3.6 mm below the AC-PC line. (b) Subthalamic somatotopic organization with overlap in the representation of feet and hand movements. The amount of units related to feet movements only (blue), hand movements only (red), and both hand movements and feet movements (green) is displayed as a stacked histogram plotted at recording locations. Units related solely to feet movements were found in the superior-mesial part only, whereas those related merely to hand movements lied in the more inferior-lateral part of the recording zone. The figure shows only 8 recording sites, because it depicts the results for the 9 patients implanted in the Right STN, and for one of these subjects this recording site did not yield any responsive units. of units holding kinematic information (93%) may also reflect an increase due to Parkinson's disease, as was demonstrated in the parkinsonian monkey, where the proportion of responsive neurons quadrupled in GPi and nearly doubled in GPe with respect to the intact one 31 . In the monkey STN, treatment with MPTP did not significantly change the proportion of responding neurons or the latencies of neuronal responses of nonoscillatory STN neurons 33 whereas in our results, large variations in time lags were observed in the human STN. However, the same study found that the average duration and magnitude of phasic responses to the application of flexion and extension torque pulses to the elbow tended to be increased in the parkinsonian monkey with respect to the intact one. These phenomena may be due to a compensation mechanism. The apparent discrepancies between the monkey and human STN may stem from inter-species differences, the different type of movement employed in each study (kinesthetic torque pulses vs. rhythmic voluntary movements), or differences in the compensatory mechanisms. In humans, somatosensory receptive fields were widened in patients with Parkinson's disease compared with those with dystonia 34 .
The overlap is also in concordance with a certain interaction between the subthalamic representations of upper and lower extremities indirectly implied by the increased fMRI activation in patients with freezing of gait (lower extremity) during bilateral finger movements (upper extremity), compared to Parkinson's disease patients who had no freezing and to healthy controls 35 .
Intra-operative single cell recordings provide a unique setup for studying the STN. To compare with the neurologically-normal STN, the inapplicability of the method to healthy humans limits comparisons to studies of other primates or lower species. The time constraints in the operating room also limit the possibility to record movement of additional body parts. Thus, some of the recorded neuronal activity may possibly be related to movement of unmonitored parts of the body.
Our results are an important step towards a more complete understanding of the pathophysiology of Parkinson's disease and its compensatory mechanisms at the single neuron level. The findings may help facilitate more focused implantation of electrodes for deep brain stimulation or subthalamic lesioning therapies, targeting specific symptoms related to kinematics, and at more accurate locations within the STN. Thus, future implantation procedures may include intra-operative mapping the STN for representation of upper and lower limb movements for implantation in a specific subregion. For example, based on this mapping, tremor dominant patients exhibiting no gait disorders may be implanted in the more inferior-lateral part of the STN to avoid the subarea containing representation of gait. Whereas gripping force was shown to be decodable from the activity of subthalamic units in patients with Parkinson's disease 1 , our findings indicate that care must be taken when attempting to generalize such reconstructions to individuals without Parkinson's disease, for example for brain-machine interfaces for paralyzed persons. Interestingly, the decoded signal explained only 68% of the variability in hand gripping force 1 , which may possibly result from reduced specificity due to Parkinson's disease. Nevertheless, decoding single unit activity 36 in the STN has the potential to help restore functional control of hand and feet movements for new types of brain-machine interfaces designed to restore mobility, hand movements and/or their coordinated actions 37 .

Methods
Subjects and Electrophysiology. Participants in this study were ten patients with Parkinson's disease who were undergoing implantation of electrodes for deep brain stimulation (3389, Medtronic Inc., Minneapolis, MN) when off medications for clinical reasons only. As part of the clinical procedure, two microelectrodes (Neuroprobe, Alpha-Omega, Nazareth, Israel), central and anterior, were temporarily implanted during surgery to locate the target brain area for implantation, based on single unit activity patterns. Bandpass filtered signals (0.3-3.0 kHz) from these microwires were recorded at 48 kHz using a 2-channel acquisition system (Neuroguide, Alpha-Omega, Nazareth, Israel). Spikes of individual cells were isolated based on their wavelet coefficients, the distribution of interspike intervals and the presence of a refractory period for the single units (that is, less than 1% of the inter-spike intervals are shorter than 3 ms) using superparamagnetic clustering ("WaveClus") 38,39 . The study was approved by the Medical Institutional Review Board at the Tel Aviv Sourasky Medical Center and all patients provided written informed consent. All experiments were performed in accordance with relevant guidelines and regulations.

Localization of Electrodes and Target Selection. MRI images co-registered to pre-operative CT scans
were used with the BrainLab navigation planning software to locate the target and planned trajectory. The recording site, in Talairach coordinates 40 , was identified based on its distance from target along the implantation trajectory. These coordinates were then converted into the standard Montreal Neurological Institute (MNI) coordinates using the template adopted by the International Consortium for Brain Mapping (ICBM152; Statistical Parametric Mapping (SPM) toolbox version 12 for MATLAB, Mathworks, Inc.).
STN target area is initially identified based on standard AC-PC coordinates (3 mm posterior to MCP, 11 mm lateral to the midline, and 4 mm below AC-PC plane) and then further refined using indirect localization based on the red nucleus as an internal landmark. In most cases, the STN could also be directly identified on proton-density weighted images.
Experimental Paradigms. Experiments were conducted intra-operatively (and thus, time was an important constraint), with the awake subject lying in a supine position, when the microelectrodes were inside the target area, the STN. Patients performed hand and feet tapping movements (defined as palmar flexion and extension at the wrist and dorsiflexion and plantar flexion at the ankle, respectively), one limb at a time, bipedally, or bimanually in an alternating pattern, each at three paces: slow, normal (i.e., self-selected) and fast. Subjects were instructed orally to start and stop each type of movement. All patients performed the task in the same order of movements. To keep movements voluntary, no external cuing was provided during the task. During the tasks, patients wore small, light, wireless movement-measurement devices (Opal monitors; APDM, Inc.) on the dorsa of their feet and hands, synchronized with the neuronal recording. The Opal monitors consist of accelerometers, gyroscopes and magnetometers. Accelerometer range is ± 6 g, with a typical noise density of 1.3 mm/s 2 /√ Hz. The gyroscopes have a range of: ± 34.9 rad/s (typical noise density: 0.81 mrad/s/√ Hz) in the X and Y axis, and ± 26.8 rad/s (noise: 2.2 mrad/s/√ Hz) in the Z axis. Magnetometers have a range of ± 6 Gauss, with a typical noise density of 160 nT/√ Hz. The sensors sample at 128 Hz.
Statistical Analysis. The kinematic parameters and firing rate function were computed in the same 10 ms bins. The latter was then smoothed by convolution with a Gaussian function (SD = 50 ms) to obtain a continuous function. We examined time lags between the neuronal firing and actual movement in the range of (− 1000) ms to 1000 ms, because the supplementary motor area, which sends inputs to the STN, can be active as early as 1200 ms before a new movement sequence is executed 41,42 . We modeled the smoothed firing rate function by a 12-fold cross-validated linear regression of each triaxial kinematic parameter: orientation, angular velocity and acceleration. The separate model for each parameter allows us to attribute the neuronal activity to a specific parameter. Orientation-based models also employed the sine and cosine of orientation in each axis. In addition, we constructed models employing all 3 kinematic parameters together. A model was constructed for each time lag independently. The time lag resulting in the maximal absolute value of the Pearson correlation coefficient between the model and the actual firing rate on a training set was selected 42 . The training set consisted of part of the tapping cycles performed by the patient. Reported correlation coefficients were calculated on an independent test set, composed of the tapping cycles that were not employed for constructing the linear model. Their p-values were Bonferroni-corrected (see "Significance of Correlation and Bonferroni Correction" below).
The entropy correlation coefficient 21 is a measure of dependence between two variables, whose properties are direct consequences of the properties of the expected mutual information measure. However, the values of the entropy correlation coefficient range between 0 and 1, with 0 indicating full independence and 1 complete dependence. This information-based measure can quantify the dependence of two variables even if their relations are non-linear. To quantify the dependence between the firing rate and kinematics, their entropy correlation coefficient was computed for each tapping cycle. All entropy computations were bias-corrected using the Panzeri-Treves term 43 . Bonferroni-corrected p-values were computed using a nonparametric two-sided sign test for the null hypothesis that the distribution of the entropy correlation coefficients comes from a distribution whose median is zero. The median entropy correlation coefficient (over all tapping cycles) and the p-value for its being 0 were computed for each time lag (between (− 1000) ms and 1000 ms, in 40 ms shifts). Among those time lags with a significant p-value (< 0.05; adjusted by the false discovery rate), the one with the highest median entropy correlation coefficient is defined as the optimal time lag. Because the distribution of entropy correlation coefficients is not normal, the Mack-Skillings test 44 served for nonparametric two-way unbalanced ANOVA.

Significance of Correlation and Bonferroni Correction.
We define that the correlation between the firing rate and its kinematics-based model is significant, when the probability p of getting a correlation as large as the observed value by random chance, when the true correlation is zero, is: p < 0.05. Reported p-values are Bonferroni-corrected with correction factor: 3 × 8 × 4 × 51 = 4,896 (3 paces (slow/normal/fast); 8 limb signals: (left foot, right foot, bipedal-left, bipedal-right, left hand, right hand, bimanual-left, bimanual-right); 4 types of kinematic models (acceleration, angular velocity, orientation, all parameters together); 51 time lags (between (−1000) ms and 1000 ms, in 40 ms shifts). Thus, the criterion is: p < 0.05/4,896 = 1.0E-5. A neuron was considered directly (linearly) encoding kinematics if the criteria for significance held for any of these tests, and in addition: r 2 > 0.30 held, where r 2 is the coefficient of determination between the model and actual firing rate.
Random Mixture of Properties. To evaluate whether or not in a population of size n, a certain combination of properties may be due to a random mixture of properties in the population, we compared the actually observed proportion of units with the combined properties (n obser /n) with the proportion of units expected theoretically under the assumption of a random mixture of properties. Under this assumption, the expected proportion can be estimated by the product of the observed proportion of units with each property in the population: n 1 /n, n 2 /n. We therefore compare (n obser /n) with: (n 1 /n) * (n 2 /n). We tested whether the actually observed population and the theoretical one have the same proportions of observations using the Chi-square test of homogeneity between a population of size n with n obser observations and a theoretical population of size n with n * (n 1 /n) * (n 2 /n) observations. Thus, Table 2 (bottom row) suggests that the probability of encoding contralateral movements of feet or hands can be estimated by: P contra = n contra /n = 37/89 = 0.42, and of encoding ipsilateral movements: P ipsi = n ipsi /n = 35/89 = 0.39. The probability of observing units encoding both ipsi-and contra-lateral movements under the assumption of random mixture of properties is therefore: P ipsi contra = P ipsi * P contra = 0.42 * 0.39 = 0.16. For a population of 89 recorded units, that should yield, on average, 89 * 0.16 = 14.24 units. Comparing the theoretical proportion with the observed 26 units, the Chi-square test of homogeneity yielded: p = 0.04, χ 2 = 4.44, 1 degree of freedom, indicating that the observed proportion is significantly different from that drawn from a population with random mixture of properties.