|Home | About | Journals | Submit | Contact Us | Français|
Brain machine interface (BMI) systems hold the potential to return lost functions to patients with motor disorders. Currently most efforts in BMI have concentrated on decoding neural activity from forearm areas of cortex to operate a robotic arm or perform other manipulation tasks. Efforts have neglected the locomotion functions of hindlimb / trunk cortex. However, the role of cortex in hindlimb locomotion of intact rats, which are often model systems for BMI testing, is usually considered to be small. Thus, the quality of representations of locomotion available in this area was uncertain. We designed a new rodent BMI system, and tested decoding of the kinematics of trunk and hindlimbs during locomotion using linear regression. Recordings were made from the motor cortex of the hindlimb / trunk area in rats using arrays of 6 tetrodes (24 channels total). We found that multiple movement related variables could be decoded simultaneously during locomotion, ranging from the proximal robot / pelvis attachment point, and the distal toe position, through hindlimb joint angles and limb endpoint in a polar coordinate system. Remarkably, the best reconstructed motion parameters were the more proximal kinematics, which might relate to global task variables. The pelvis motion was significantly better reconstructed than any other motion features.
IT has been suggested that brain machine interfaces (BMI) may potentially be used in order to bypass the injured site of the spinal cord. The aim of such a bypass is to replace the lost functions and to help rehabilitate the residual functions. To date, most efforts in BMI have concentrated on recording from different areas of cortex in order to perform upper limb manipulation tasks. With a single exception  all of the BMI research has largely neglected the hindlimb / trunk cortex area for restoring locomotion -. Such research is essential because neural circuits for control of forelimbs differ in various ways from those that control the trunk / hindlimb, and both neural recording types may figure in various therapies and restoration approaches , . We sought to discover how motor cortical neurons encode the different hindlimb / trunk kinematic variables in locomotion under normal and load conditions, and which variables could best be decoded. Some of these are clearly suitable for the control of functional electrical stimulation or artificial devices or for epidural stimulation in the rodent models that are a staple of spinal cord injury research. To this end, this letter presents a first analysis of cortex recordings in rats walking on a treadmill obtained with our new hindlimb targeted BMI system.
Four adult Sprague-Dawley rats, weighing 250~300g, were used in these experiments. Care and treatment of the animals conformed to protocols approved by the ULAR and IACUC of the College of Medicine of Drexel University. After anesthesia using KXA (63 mg/Kg Ketamine, 6 mg/Kg Xylazine, 0.05 mg/Kg Acepromazine), rats were implanted with a pelvic orthosis, which allowed a PHANTOM robot (Sensable Devices) to apply forces directly to the skeleton, as shown in Fig. 1(a). Tetrode arrays, which consisted of 6 tetrodes positioned around the neck of an etched spear-shaped tungsten rod, were implanted stereotactically in the hindlimb / trunk area of the motor cortex (2.0mm lateral of the midline, 1.6mm caudal to bregma and 1.0 -1.5mm in depth). Intracortical microstimulation of rat M1 in this recorded region generates hindlimb/trunk movements in all rats tested in our hands. At the end of the final recording session, strong electrolytic lesions were made and brain slices were cut in 40um sections [see Fig. 1(c)]. Neural data were recorded using a Cerebus system (Cyberkinetics, Inc. / Blackrock systems) after animals recovered from surgery. Multiple single units were isolated offline using the Plexon system implementation of the Shoham algorithm for tetrode recordings . Up to 24 channels of neural activity at a time in a single array could be recorded, and 1 or 2 individual cells from most wires of the tetrodes could be recorded, as seen in the Fig. 1(d), around 15 to 38 well isolated single units were used in each trial.
The robotic system employed has been described in detail previously . Briefly, different force fields could be applied to the pelvis of the rat, using the robot, updated at a frequency of 1 KHz, and coupled with neural and video recordings. Rats were trained to walk on a treadmill at a speed ranging from 0.1m/s to 0.15 m/s. Unidirectional (vertical forces only, i.e., radially anisotropic) elastic fields were used in experiments here. Elastic loading (downward) fields were applied by setting an equilibrium point on the horizontal equilibrium plane of the elastic field 12.5 mm under the rat’s normal pelvic height with a field stiffness of 44 N/m. Reflective markers were attached to the toe, iliac crest, hip joint, and ankle joint. Marker movements were recorded with a digital camera (PULNiX America Inc.) at a sampling rate of 100 Hz. The side of the recording hindlimb was contralateral to the side of the electrode placements recording cortical cells. Given the coordinates of hip marker and ankle marker and the lengths from the knee joint to the two markers, the knee joint positions could be calculated. The joint angles were calculated as
where Li, Li+1 are vectors formed by adjacent pairs of the three consecutive markers, as shown in Fig. 1(a). In addition to deriving a joint coordinate system, we also calculated the toe marker relative to the hip marker in a polar coordinate system (ρ, θ), as in . The movement related variables were smoothed with low pass filter of zero-phase distortion (<25Hz, 6th order Butterworth, ‘filtfilt’ in MATLAB). Robot data were down sampled to 100 Hz and smoothed with the same low pass filter. All collected data were synchronized by a trigger event from the start of the experiment. A total of 22 normal locomotion trials and 20 loaded locomotion trials in which an elastic load was applied at the pelvis were obtained on different recording days from the 4 rats. From these data statistically near-stationary intervals of locomotion were selected. There were about 50-100 steps in each trial interval used.
Stride to stride kinematics statistics were calculated by normalizing individual stride durations, and dividing the stride interval into 100 segments. The mean cycle value of each spatial or angular variable in each segment in the interval and its standard deviation were calculated in the standard manner. The 100 standard deviations were then averaged to obtain a measure of average deviations around the mean cycle pattern for the entire interval. We next expressed this average standard deviation across the cycle as a percentage of the range of the mean values obtained in the interval. This percentage measure captured how much cycle to cycle variation occurred in relation to the depth of modulation of the parameter in the repeatable gait cycle pattern. This allowed us to compare the degree of variation of different angles and of different spatial variables and their relations to the R-squared values obtained in their regressions to neural firing patterns.
In order to make a preliminary test of the feasibility of decoding movement related variables from neurons in primary motor cortex of hindlimb/trunk area in the rat, simple offline decoders were constructed by regressing individual movement related variables on firing rate as
where Ni(t+ τ) is a firing rate vector of cell i at time t with time delay τ (negative corresponding to the past); Y(t) is a movement related variable representing robot data or hindlimb kinematics at time t; Wi(τ) is a weight vector for neuron i at time lag τ; C is a constant intercept, and ε(t) is a residual error. The continuous firing rates were calculated by averaging the spikes in a preceding window, and expressed at 10 ms intervals. After testing with different window lengths, we fixed window length to 250 ms as most kinematic parameters had better performance in this duration window and ensemble average regression R-squared values were highest. It should be noted that for some slow kinematic parameters longer windows could also have good performance. The rate time series was then subject to the same low pass filter used for kinematics smoothing. In the regression, the linear weight filter Wi(τ) obtained was fully causal, with L=25. That is, we only used information from the past 250 ms (L=25) of the filtered rate time series to predict the current kinematic states. Regression parameters were obtained using the MATLAB function ‘regress.’ The R-squared value of the linear regression, which was calculated from the observed and the fitted data, represents the proportion of explained variance in the movement related variables predicted by the neural population activities, and was used as a criterion for performance evaluation. We then tested prediction quality for a subset of the test datasets (9 normal and 11 load field). We used data sets where cells’ firing rates were higher than 1 Hz and where such active cell numbers were greater than 15. The predictive models were fitted with 90 percent of the data in each session and then tested for prediction quality with the remaining excerpted 10 percent. The variance of prediction quality was then estimated by a 9-fold cross-validation, repeating this process with successive 10 percent excerpts. R-squared values indicate a strong potential for prediction and strong correlation between the movement related variable and cell firing rate.
We recorded neural activities in the motor cortex with the tetrode arrays, and simultaneously multiple movement related variables: robot attachment-point position, velocity, and the hindlimb state. These were expressed as marker positions and velocities (in world or body centered coordinates), and were also expressed in joint-based or hip-centered polar coordinate systems. These movement related data were estimated simultaneously from the neural data off-line, with a simple linear decoder, as animals walked under normal and load conditions.
The R-sqr values of each variable were higher for fitted datasets than for prediction testing, while the trend for variance accounted for in the prediction testing of different variables was consistent with the simple fitting. The decoding of kinematic velocities for both robot / pelvis attachment and the hindlimb were not as good as decoding of robot / pelvis position and hindlimb posture (joint angles, positions) under both normal and load locomotion conditions (paired t-test, p<0.05) for most of the parameters. The reasons might be that rats in the locomotor task were more attentive to hindlimb stance postures, than to hindlimb swing velocity. This result was similar to some forearm manipulation tasks in primates . Our analysis here focuses on the position coding.
We found that proximal variables were best represented and predicted in rat trunk/hindlimb cortex. There was a session to session variance for the prediction R-squared values, as there might be different neurons recorded in different sessions recorded in different days. On all days, however, the R-squared values were significantly higher for decoding the lateral position of pelvis-robot attachment point than for decoding distal hindlimb joint angles in joint-based coordinate system and polar radius / leg length in limb-based coordinate system. The stride-to-stride average percentage standard deviations of all the robot positions (proximal trunk) were significantly higher than any of the more distal parameter’s standard deviations (paired t-test, p<0.05), though lower than the more imprecise video estimates of the proximal hip and iliac crest position markers. Thus the higher R-squared values of these pelvic motions when related to neural activity likely did not simply represent more repeatable pelvic behaviors. This proximal bias of neural decoding held in both normal and load conditions (paired t-test, p<0.05). Fig. 2 shows sample reconstructions of a rat and clearly shows that the estimation of pelvic motions at the robot end-point was more accurate than the estimation of joint angles or foot position. In polar coordinates the polar orientation of the hindlimb is better represented than polar radius / leg length and the individual joint angles. In joint coordinates the proximal joint angles were significantly better represented that the distal joint angles (paired t-test, p<0.05). The decoding of hip joint angle was the best among all the three hindlimb joint angles in joint-based coordinates. In the case of hip angle and limb orientation, in contrast to pelvic translation, these more proximal parameters did show lower absolute standard deviations and lower standard deviations as a percentage of range, compared to the more distal joint angles. The prediction of limb endpoint in the polar (limb-based) coordinate system was better than that represented in the joint-based coordinate system (paired t-test, p<<0.05). This result was in accordance with reports of decoding from sensory neurons in the dorsal root ganglion . The presence of loads appeared to slightly improve kinematic encoding, although the average stride-to-stride standard deviation was increased and was bigger in the load condition than in the normal condition.
To summarize, limb joint kinematics were on average better predicted for proximal joints, and this was consistent with lower variance around the gait cycle motions at the proximal joints. In contrast, pelvic translation parameters, which showed significantly less repeatable motions, and higher standard variation as a percentage of mean range, were also significantly better related to recorded neural firings than were distal markers, or the limb joint angles with more repeatable motions. These data suggest this cortical area best encodes proximal task-related variables in this task set.
We have tested if multiple movement related variables can be fairly accurately decoded by using simultaneous activity of populations of cells in the hindlimb / trunk area of the rat motor cortex during locomotion. Our data show that in the areas we recorded the best predictions were obtained for proximal variables, whether pelvic translation, or limb joint angles. Proximal variables are extremely important in determining limb endpoint behavior, and trunk translation is likely to be a highly significant task variable. The high R-squared value with neural firing obtained for proximal kinematics was not likely due solely to the smaller variability in the kinematics of some of these variables during locomotion. The best reconstructed of these proximal variables often showed larger standard deviations than more distal variables, including when deviations were expressed as a percentage of mean range of variables across a stride. Further, load increased standard deviations of all variables, while improving proximal R-squared values. Consistent with these load effects, we speculate that during the step cycle the rat hindlimb / trunk cortex may be most interested in representing and controlling stance, not the faster but more kinematically stereotyped swing phase, and this may account for the differences in quality of position and velocity regressions we observed. Rodent hindlimb / trunk cortex in our recorded area may be most engaged in representing such task related variables during locomotion.
By using only task level open-loop controllers, bio-inspired hexapod robots can exhibit stable dynamic locomotion over unstable terrain . It could be possible to employ rodent trunk cortex activity to classify each step cycle into its different phases or estimate desired proximal states and then to decode these states for open loop task level control of a prosthetic. A significant and related role of rodent cortex in recovered weight support after neonatal spinal transaction is suggested by our published data, and lesion of this same area alters trunk roll in locomotion .
For prediction, in contrast to fitting, we only tested data with higher firing rates and larger numbers of such recorded cells. Clearly, the prediction quality is increased when such cell numbers are increased, while low firing rate cells will contribute less to the decoding performance. Still higher numbers of neurons would be needed for good prediction, and we observed significant shrinkage in the cross-validation. However, our qualitative differences in decoding among variables, on which we focus here, remained similar for fitting and prediction. Compared with the performance of decoding reported for recordings from sensory DRG neurons , the R-squared value of decoding limb motion from cortical neurons here was lower. However, firing rates were much lower in our cortical neurons. The afferent cells directly encode the skin or muscle receptors, while the projection to cortical neurons from the receptors is less direct and more abstracted. We speculate that the reason why cells were more related to pelvic motion than to the hindlimb joint angles is that the pelvic motion is of crucial importance for locomotion as well as coordination of the four limbs; these variables might represent a task level representation and control strategy. Cortex might closely monitor or play a role in control of these task variables. Indeed, we observed that loading slightly improved the prediction quality. Cortex might not be as strongly involved in moment-to-moment control during normal locomotion. However, under load conditions walking initially may need more attention, and cortical engagement in the activity might increase. This might account for the higher R-squared value in the load condition compared with normal locomotion. The patterns of activity and representation here might be specific to quadrupedal locomotor patterns. The biomechanical and control requirements of quadrupedal gaits may have significant differences from bipedal gaits and this difference may account for differences in our findings in relation to findings in . Further, these may represent differences between evolutionarily older default patterns (in rats here) and newly learned and highly trained locomotor patterns (in primate bipeds in ). In this regard, it may be very instructive to compare these data to the cortical activity patterns in adult rats spinalized as neonatal rats, where locomotion may require more cortical ‘skills’ and engagement to manage the locomotor task [12,13,14].
In summary, our data suggest that the trunk / hindlimb cortex in rats may be most involved in locomotion at the task level, representing proximal and task level variables, while leaving the detailed execution of locomotion to the low spinal level, except when corrections or adaptations become needed in a new context. The data presented here suggest that the activity in this part of cortex, at least in rats, will be best utilized to control proximal limb and trunk motions.
We thank anonymous reviewers, Dr. Emery N. Brown and Iahn Cajigas for valuable suggestions.
This work was supported by funds from NIH NINDS NS44564, NS54894, the Neilsen Foundation and the PA Dept. of Health (Tobacco Settlement research award).
Weiguo Song, Department of Neurobiology & Anatomy, Drexel University, College of Medicine, PA 19129 USA (Email: firstname.lastname@example.org).
Arun Ramakrishnan, School of Biomedical Engineering, Drexel University, PA 19104 USA (Email: ude.lexerd@88ra).
Ubong. I. Udoekwere, School of Biomedical Engineering, Drexel University, PA 19104 USA (Email: ude.lexerd@22uiu).
Simon. F. Giszter, Department of Neurobiology & Anatomy, Drexel University, College of Medicine, PA 19129, and with School of Biomedical Engineering, Drexel University, PA 19104 USA (Email: ude.lexerd@33gs, phone: 215-991-8412, fax: 215-843-9082).