|
|
Supplementary MaterialGeneral Setup The stereo image was controlled by an SGI Octane workstation and viewed through polarized lenses and a 96 Hz light-polarizing shutter screen. Three-dimensional wrist position was sent to the workstation at 100 Hz from an Optotrak 3020 motion tracking system. Cortical activity was collected via a Plexon Data Acquisition system. Spike times were transferred to the workstation, and a new brain-controlled cursor position was calculated every ~30msec. Open-loop vs. Closed-loop Experiment The left arm was restrained while the right arm was free to move during both hand- and brain-controlled movement blocks (Fig. S1). Cursor radius was 1 cm. Target and center radii were 1.2 cm. A liquid reward was given when the cursor boundary crossed the target boundary for ~300 ms or more. Radial distance (center start position to center of target) was 4.33 cm under brain-control. Since hand-controlled movements were quick, radial distance was increased to 8.66 cm during the hand-controlled movement blocks to increase the duration of cortical data collection. For a fair comparison, offline predicted trajectory (open-loop) hit rates were calculated with targets at the on-line brain-controlled distance of 4.33 cm. Each day's open-loop trajectories calculated offline were scaled, so the median radial endpoint distance was also 4.33 cm. A modified population vector (1, 2) was used to translate cortical signals into trajectories. Modifications to the population vector included scaling each cell's contribution by a measure of its quality of tuning-the R2 value from its preferred direction regression equation raised to a power of 0.8 (in previous offline tests, this value maximally reduced endpoint error). Additional scaling equalized the expected magnitude of movements in the X, Y, and Z directions. The on-line brain-controlled velocities were further scaled down to ensure only reasonably direct trajectories could get to the target with in the 10-15 sec time limit. Hand and brain-controlled movements were performed in alternating blocks of movements to all eight targets. Under brain-control, missed targets were repeated up to 5 times per block for practice. Performance statistics were compiled by first averaging statistics from repeated targets together in each block to ensure statistics represent all targets equally, even if the number of movement samples differed between targets. Monkeys 'L' and 'M' were chronically implanted in the motor and premotor areas of the left hemisphere with arrays consisting of fixed stainless steel and/or tungsten microwires insulated with Teflon/polyimide. Monitoring cortical activity during manual arm manipulation showed both animals' had units related to proximal and distal arm movements. Monkey 'M' also had some units related to upper back activity. Many electrodes detected waveforms from multiple cells, some of which could not be individually isolated. The terms 'cell' or 'unit' used in this paper refer to both individually-isolated cells and inseparable multi-cell groups. Often channels suspected of having very low amplitude cell activity (within the noise floor) were included in each day's experiment. Since each unit's X, Y, and Z weighting coefficients were scaled based on tuning quality, these noisy units usually had weights too small to affect the population vector. 'Cells used' in Table 1 is the number of units whose vector sum of the absolute value of their scaled weights accounted for at least 95% of the vector sum of the absolute values of all scaled weights. At the beginning of each day's experiments, 8-10 minutes of center-out arm movements and cortical data were collected to measure the baseline movement-related behavior of the recorded cells. This was used to calculate the preferred directions used in the modified population vector. R2 values shown in Table 1 were calculated from the baseline data for each cell by regressing firing rates, averaged by target, with target direction. Since new preferred directions and scaling factors were calculated each day, differing estimation errors (as well as any actual changes in the tuning properties) slightly changed the form of the population vector from day to day. Most of monkeys 'L' and 'M's units had fairly stable waveforms and tuning properties between days. Larger variations between days in the modified population vector coefficients sometimes happened due to intermittent hardware problems making some recorded channels noisy or unusable. Between the slow and ballistic brain-controlled experiments, subject 'M' had additional arrays implanted in the right motor and premotor cortices which detected units related to both proximal and distal left arm movements. In the ballistic open- vs. closed-loop test, the right arm was still free to move while the left arm was restrained, just as it had been in the slow task. Many of the right hemisphere cells showed significant tuning during right arm movements with the left arm restrained. Therefore, units spanning both hemispheres were used in the modified population vector (44 ± 4). Both brain- and hand-controlled movement blocks where at a radial distance of 8.66 cm. The allowed movement time was reduced to 800 msec. Target radius for the hand-controlled blocks was still 1.2 cm, but the target radius for the brain-control blocks was increased to 3.5 cm to decrease the task difficulty. Offline brain-predicted hit rates (open-loop) were calculated assuming a target radius of 3.5 cm as well. Each day's offline trajectories were scaled to make the median radial end-point distance match the radial target distance of 8.66 cm.
Coadaptive Brain-Control Task
The coadaptive algorithm was designed to be used with immobile patients with no prior knowledge of the recorded cell's tuning properties. In this case, the algorithm would start with random initial parameters and iteratively refine them as the patient attempted to make a series of brain-controlled movements. In these animal experiments, we tried random numbers for starting values, but we also attempted to speed up the coadaptive process by starting with the actual tuning properties obtained from that day's baseline hand-controlled recordings. This, however, did not improve initial performance or speed up the adaptation process. Cell-tuning functions obtained during normal arm movements were not good predictors of intended movement once both arms were restrained. The movement task was a 3D center-out task with eight targets (8.66 cm distance from center). The cursor (1 cm radius) was held in the center until 200 ms after a target appeared. Then the cursor was moved under brain-control until the target was hit or the time limit ran out (see algorithm details below). Because initial performance was so poor, we started each day's experiment with large, easy-to-hit targets (4 cm radius). As coadaptation progressed, and the subject's performance improved, the target size was decreased to encourage the development of more directional accuracy (range limited to between 4 and 1.2 cm radius; radius adjusted based on the accuracy of the most recent three blocks of movements). At the end of each day's experiment, tuning functions and statistics were calculated directly from the cortical discharge activity collected during 15-17 consecutive minutes of the best brain-controlled movements made at the 3 sec/movement time limit (see algorithm details below). Mean firing rates per movement were calculated for each unit over a fixed time window in the initial part of each movement starting at ~100 ms after the target appeared. If the target was hit before the end of the brain-control calculation window, the time until target hit was used instead. This was done for all brain-controlled movements whether they successfully hit the target or not. Mean firing rates per movement were also calculated from each day's baseline hand-controlled movements starting at ~100 ms after the target appeared and ending when the target was hit. Brain-control mean rates per movement were calculated over two different window durations-420 and 990 msec. The 420 msec rate calculations were used in all hand- and brain-control comparisons where variability was an important factor in the comparison (i.e. Fig. 1A, Fig. 2A-D, Paper Note 12). The 420 msec window was chosen because it matched the mean window length used in the hand-control calculations. In order to reduce estimation errors, the 990 msec brain-control rates were used in the calculation of preferred directions (Fig. 1B-D). Results reported in Figures 1 and 2 were calculated as follows. In Figure 1A, rates were normalized by the mean (dotted lines) and boxes were scaled to ±30 spikes/sec. In figures 2 C & D, each cell's rates were averaged by target, and the results from all targets were normalized by subtracting the mean and dividing by one standard deviation. ( is the angle between the cell's preferred direction and target direction. Lines in Figures 1B and 2B were based on regressing the individual cell's data points with days of practice-not just the daily averages shown in the plots. At the time of the coadaptive task, monkey 'M's cells were very stable from day to day (based on stability of waveforms and hand-controlled preferred directions)(3). The new subject, monkey 'O', had 64 Teflon/polyimide coated stainless steel microwire electrodes chronically implanted in the motor and premotor arm areas of the left hemisphere. Cell activity modulated by proximal and distal arm movements was observed. Individual cell recordings were very unstable from day to day and disappeared all together before the end of the experiment. This may be due, in part, to changes in the cells themselves, but it was most likely due to changes in our ability to record the activity (e.g. due to electrode motion and/or encapsulation). This instability made it difficult to track or train individual cells over time. Therefore, analysis of day-to-day learning was done only with monkey 'M's data. Usually, not more than 20% of the cortical units recorded were well-isolated individual cells in monkey 'M'. About half of the remaining units were multi-cell groups with overlapping wave forms (such as the light and dark gray wave forms in Fig. 1A). The other half of the remaining units were low amplitude 'noise' or 'hash' units comprised of background cell activity within the noise floor. All types of units were input into the coadaptive algorithm to see what would become useful with training. Many of the 'noise' units showed some significant tuning in both tasks, but tended to have lower R2s than the well isolated cells. The coadaptation process reduced the cell's weights based on their level of tuning and, therefore, reduced their effect on the movement. '# Units used' in Table 2 is the number of units whose vector sum of the absolute value of their average X, Y, and Z weights accounted for at least 95% of the vector sum of the absolute values of all averaged X, Y, and Z weights (averaged across positive and negative weights and across all blocks during the segment with a 3 second movement time limit, see algorithm details below). The magnitude of rotation in preferred direction between the two tasks was not related to whether the unit was a single cell, multi-cell group, or a tuned noise unit. All types of tuned units exhibited a wide range of changes in preferred directions between the hand-and brain-control tasks.
During the last two weeks of the coadaptive experiment, surface EMGs were collected at 1000 Hz over four days during both hand-control and brain-control tasks in monkey 'M'. Root mean squared values were low-pass filtered at 30 Hz. For each muscle and task, a baseline median level (SIGo) and three standard deviations (NOISEo) were calculated from data collected prior to all movements (from 100 msec before to 100 msec after the target appeared). Median values where also calculated by target at each time step during the movement, and a peak movement value (SIGm) was calculated for each muscle as the maximum of a 50 msec moving average of these median values. Relative magnitude of muscle modulation was compared between the hand-control (HC) and brain-control (BC) tasks by taking the ratio
Supplemental Figure 1. Real-time 3D brain-control experimental setup. (A) A mirror in front of the monkey's face reflects a 3D stereo image of the moving cursor (yellow sphere) and static target (blue sphere) projected from above. On the other side of the mirror, the monkey moves one arm with a position sensor (red) taped to the wrist. The 3D position of the cursor is determined by either the position sensor ("hand-control") or by the movement predicted by the subject's cortical activity ("brain-control"). The cursor and target are shown with dotted outlines in the monkey's workspace, but do not physically exist. The monkey sees the cursor and target as being in the workspace but cannot see its own arm. The movement task is a 3D center-out task. The cursor is held at the central position until a target appears at one of eight radial locations which form the corners of an imaginary cube. The center of the cube is located distal to the monkey's right shoulder. (B) shows the approximate relative cube and target size during the hand-control task from the same observer view point as in (A).
Medium version | Full size version
The Constant Parameter Prediction Algorithm (CPPA) Task: In the CPPA task, the positive and negative X, Y, and Z weights used to calculate movement were held constant throughout the task (see algorithm details below). Weights used were an average of the parameters from 2-10 of the most successful consecutive movement blocks in that day's coadaptation task. The task was designed to test if subjects could generalize their movement skills acquired in the coadaptive task. Instead of just center-out movements, the subjects had to move the cursor from the center to the target, then stop and move the cursor back to the center under brain-control to get a reward. In addition, six novel target positions were added-straight up, down, left, right, forward, and backward. All target radii were 2.0 cm and the radial distance to all targets' centers was 8.66 cm. Movement time per trial was limited to 4 seconds each way (Vscale = 33mm/sec, see algorithm details below). Monkey 'M' did the CPPA tasks on the last 13 days of the experiment when he was consistently able to do the coadaptive task at the highest precision level (i.e. smallest target size) allowed (task skipped on day 34). Monkey 'O' did the CPPA task for only 5 days before his cells were no longer recordable. Both monkeys improved their performance, but monkey 'O's improvement was not significant with only 5 days' data. ALGORITHM DETAILS: Real-Time Movement Prediction Algorithm Used in the Coadaptive and CPPA Tasks
Calculating Movement in X at time t, Mx(t):
Each cell had two 'weights' Wxpi and Wxni. If NRi(t) was positive it was multiplied by Wxpi. If it was negative, it was multiplied by Wxni. These products were then summed across all cells.
 iWx(n_or_p)i * NRi(t) = X(t)
i (Wxpi - Wxni) * Ek[|NRi(k)|]/2
In the coadaptive task, Ek[|NRi(k)|] was calculated from the normalized rates observed in the most recent complete block of brain-controlled movements to all targets. In the CPPA task, it was calculated by averaging |NRi(k)| over the most recent 3 minutes of data and was updated every 500 msec. The expected drift was then subtracted from X(t) to get the predicted movement in X at time t:
Movements in Y and Z were predicted in a similar manner using the appropriate weights (i.e. Wypi, Wyni, Wzpi, Wzni). Additional normalization was then used to set the overall velocity scale of the movements, and to ensure the X, Y, and Z movement components had a similar average magnitude.
My(t) = Vscale*my(t)/Ek[|my(k)|] Mz(t) = Vscale*mz(t)/Ek[|mz(k)|] Ek[|mx(k)|] was the expected value or average of the absolute value of mx(k) taken over all calculation times, k, in a large time window. Time windows used were the same as those used to calculate Ek[|NRi(k)|]. Ek[|my(k)|] and Ek[|mz(k)|] were similarly calculated. Weight Updating Procedure Used in the Coadaptive Task: The following is the adaptive procedure used to adjust the weights Wxpi, Wxni, Wypi, Wyni, Wzpi, and Wzni at regular intervals: Initial Values: With an immobile patient, initial weight values can simply be set to small random numbers. We set the initial weight to either random numbers or to the cell's preferred directions calculated from that day's baseline hand-controlled movements. Unlike the traditional preferred directions calculated for the analysis shown in the journal article, separate positive and negative preferred direction weights were calculated for use as initial weights. These were calculated by separately regressing positive and negative normalized rates with the target directions (Dx, Dy, Dz), then normalizing the coefficients to a unit vector:
Only Negative NRi(t) = wo + wxniDx + wyniDy + wzniDz
Wxpi = wxpi / (wxpi2 + wypi2 + wzpi2) 1/2
Wxni = wxni / (wxni2 + wyni2 + wzni2) 1/2 Using these preferred directions, however, was not beneficial as the weights converged just as fast when initial values where set to small random numbers. Initial Ek[|NRi(k)|] values were also derived from the normalized rates observed during the baseline hand-controlled movements. Then, mx(k), my(k), and mz(k) were calculated for each time step, k, of these initial movements using the cells' positive and negative preferred directions as weights. Ek[|mx(k)|], Ek[|my(k)|], and Ek[|mz(k)|] were then derived from mx(k), my(k), and mz(k). In the immobile patient, initial Ek[|NRi(k)|] values can be derived using rate data recorded while the patient is at rest or is thinking about arm movements. Ek[|mx(k)|], Ek[|my(k)|], and Ek[|mz(k)|] values can initially be set to one, since they will quickly change to more appropriate values once the adaptation process begins. The above steps determined all initial parameters needed to start the adaptive phase Adaptive phase for iteratively improving parameters: During adaptation, the targets were presented in 'blocks'. A block consisted of all eight targets presented in random order. In this task the cursor first appeared at a central hold position. Then, after a 50 msec pause, a target appeared. After another 200 msec pause (allowing the subject time to become aware of the target) the cursor movement was under brain-control for a fixed amount of time or until the target was hit-which ever came first. The brain-controlled movements, Mx(t), My(t), Mz(t), were calculated as described in the 'real-time movement prediction algorithm' section above. Missed targets were repeated up to 3 times per block to ensured changes in weights more strongly emphasizes correcting for troublesome target directions. The following equation was used after each block to determine the weights for the next block. Terms used in the equation were calculated from the data collected during the brain-controlled parts of the most recently-completed block, S. Equations for Wxpi are shown. Similar equations exist for Wypi, Wzpi, Wxni, Wxni, Wxni.
 Wxpi(S)) + (1-Ps)Wxpi(Sbest) )
Wxpi(S) is the adjustments to the weights suggested by the errors seen during block S. It is calculates as:
Wxpi(S). It is calculated as:
Axpi = Ao( 1 + CA1(N[EMxpi(S)] + N[ECxpi(S)]) )
The term CA1(N[EMxpi(S)] + N[ECxpi(S)]) allows for each cell's Axpi to be increased or decreased individually depending on that particular cell's error magnitude, EMxpi(S), and error consistency, ECxpi(S). The operator, N[], normalizes the enclosed term across all cells to between -1 and 1. We normalized by subtracting out the median, dividing by 2 standard deviations, and truncating values which exceed the bounds to -1 or 1. CA1 is a small constant. We tried values ranging from 0.05 to 0.4. When cells have large and consistent errors, (N[EMxpi(S)] + N[ECxpi(S)]) equals 2, and Axpi is increased over Ao. When cells have small and inconsistent errors (like you would expect when the weights are just right), (N[EMxpi(S)] + N[ECxpi(S)]) equals -2, and Axpi is decreased over Ao. However, when cells have large and inconsistent errors, their entire weight, Wxpi(S + 1), should be scaled down in the next block. This is done through the Bxpiterm.
Bo = 1 - En[Phit(n)] (n: S-0, S-1,...,S-N + 1) The scale factor, Bo, goes to 1 if the subject has done poorly over the last N blocks, and goes to 0 when the subject has been doing well. If the subject has been doing well, Bxpi then goes to 1 and has no effect on the weights. CB is a constant usually set between 0.0 and 0.5. We tried values between 0 and 0.45. If the subject has been doing poorly, Bxpi can increase or decrease individually depending on that cell's performance in the most recent block. Bxpi decreases below 1 for cells with large inconsistent errors. Bxpi increases above one for cells with small consistent errors. If errors are small and inconsistent, or errors are large and consistent, Bxpi remains at about 1. The error magnitude measure, EMxpi(S), and error consistency measure, ECxpi(S), are calculated as follows:
ECxpi(S) = | Ek[ Wxpi(k) NRi(k) -(Tx(k)-Cx(k)) ]|/ EMxpi(S)
The above procedure is repeated block after block with new weights determined after each block. Also, the parameters Ek[|NRi(k)|] for the next block are calculated using data collected during the current block, S. Parameters Ek[|mx(k)|], Ek[|my(k)|], and Ek[|mz(k)|] are also calculated by block, but the values used in the next block are the weighted combinations of the current block, S, (scaled by Ps) and the recent best block, Sbest, scaled by 1-Ps. To further encourage improvement in the weights, the size of the targets can be reduced as the subject gets better at making brain-controlled movements. We started with a 4cm radius target and reduced it by 1mm in the next block if the average percentage of targets hit in the last 3 blocks exceeded 70%. Conversely, we also enlarged the radius by 1 mm in the next block if the percentage was less than 70%. Target radius was not allowed to go above 4 cm or below 1.2 cm. The time allowed to make each movement was limited. This encouraged straighter trajectories because the average distance the cursor could move in this time was held approximately constant by the normalization terms, VScale/Ek[|mx(k)|], VScale/Ek[|my(k)|], and VScale/Ek[|mz(k)|]. Weights which produced straighter paths had an increased chance of hitting the target in the time allowed. Therefore, a larger portion of these weights were then conserved in subsequent blocks through the (1-Ps)Wxpi(Sbest) term. We started the task with a 1 second movement time limit and Vscale = 100 mm/sec. This allowed many blocks (and, therefore, weight adjustments) to be completed in a short amount of time. After 10-15 minutes, we tripled the movement time to 3 seconds and cut the Vscale down to one third the value, 33 mm/sec. This allowed for more on-line error corrections by the subject, and helped further refine the weights to fit the brain's modulation capabilities.
1. A. P. Georgopoulos, R. E. Kettner, A. B. Schwartz, J. Neurosci. 8, 2931 (1988).
|
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Science. ISSN 0036-8075 (print), 1095-9203 (online)
News | Science Journals | Careers | Blogs and Communities | Multimedia | Collections | Help | Site Map | RSS
Subscribe | Feedback | Privacy / Legal | About Us | Advertise With Us | Contact Us
© 2009 American Association for the Advancement of Science. All Rights Reserved.
AAAS is a partner of HINARI, AGORA, PatientInform, CrossRef, and COUNTER.
You have reached the bottom of the page. Back to top