Area 2 of primary somatosensory cortex encodes kinematics of the whole arm

Proprioception, the sense of body position, movement, and associated forces, remains poorly understood, despite its critical role in movement. Most studies of area 2, a proprioceptive area of somatosensory cortex, have simply compared neurons’ activities to the movement of the hand through space. Using motion tracking, we sought to elaborate this relationship by characterizing how area 2 activity relates to whole arm movements. We found that a whole-arm model, unlike classic models, successfully predicted how features of neural activity changed as monkeys reached to targets in two workspaces. However, when we then evaluated this whole-arm model across active and passive movements, we found that many neurons did not consistently represent the whole arm over both conditions. These results suggest that 1) neural activity in area 2 includes representation of the whole arm during reaching and 2) many of these neurons represented limb state differently during active and passive movements.


23
Moving around in an uncontrolled environment is a remarkably complex feat. In addition to the 24 necessary computations on the efferent side to generate movement, an important aspect of motor 25 control is the afferent information we receive from our limbs, essential both for movement 26 planning and for the feedback it provides during movement. Of the sensory modalities we 27 receive, proprioception, or the sense of body position, movement and forces, is arguably the most 28 critical for making coordinated movements Gordon et al. 1995;29 Sainburg et al. 1995;Sainburg et al. 1993; Sanes et al. 1984). However, despite its importance, 30 few studies have explicitly addressed how proprioception is represented in the brain during 31 natural movement. In comparison, touch, vision, and the motor areas of the brain have received 32 far more attention. 33 Our conscious experience of proprioception typically focuses on where our hands are going, 34 rather than the rotations of our joints. This also matches with psychophysical data showing that 35 humans are better at estimating the location of the hand than estimating joint angles (Fuentes and 36 Bastian 2010). Perhaps consequently, of the few studies that examine central nervous system 37 (CNS) activity underlying proximal limb proprioception, most assume that individual neurons 38 represent the movement of the limb's endpoint, rather than joint angles or muscle lengths (Bosco 39 et The proximal arm representation within area 2 of primary somatosensory cortex (S1) receives a 44 combination of muscle and cutaneous information (Hyvärinen and Poranen 1978;Padberg et al. 45 2018; Pons et al. 1985) and is thought to be important for proprioception (Jennings et al. 1983; 46 Kaas et al. 1979; London and Miller 2013). Interestingly, recent computational studies have 47 shown that while neural activity may appear to be tuned to the state of the limb's endpoint, 48 features of this tuning might be a direct consequence of the biomechanics of the limb 49 (Chowdhury et al. 2017; Lillicrap and Scott 2013). Consistent with those results, we have 50 recently observed, using artificial neural networks, that that muscle lengths were better predictors 51 of area 2 activity than were hand kinematics (Lucas et al. 2019). Thus, in this study, we set out to 52 understand what information is represented in area 2, with two experiments that altered the 53 relationship between hand and whole limb kinematics. Using both classic single neuron analysis 54 techniques like tuning curves and preferred movement directions (Georgopoulos et al. 1982;55 Prud'homme and Kalaska 1994; Sergio and Kalaska 2003), as well as more recently developed 56 neural population analyses (Churchland et al. 2012; Cunningham and Byron 2014), we 57 discovered several features of neural activity that could not be explained by the classic hand-58 based model. However, using biomechanical modeling, we show that models relating neural 59 activity to kinematics of the whole arm can explain these features. Our results indicate that if 60 there is a transformation towards representing reaching in terms of only the hand, it likely occurs 61 beyond area 2. 62

63
For the experiments detailed in this paper, we recorded neural signals from three monkeys (H, C, 64 and L) using Utah multi-electrode arrays (Blackrock Microsystems) implanted in the arm 65 representation within Brodmann's area 2 of S1 ( Figure 1A). We trained each of these monkeys to 66 grasp a two-link planar manipulandum and make reaching movements to targets presented on a 67 screen in front of them ( Figure 1B). During these behavioral sessions, we also tracked the 68 locations of markers on the monkey's arm using a custom motion tracking system and registered 69 these markers to a musculoskeletal model in OpenSim (SimTK) to estimate joint angles and 70 muscle lengths. Our experiments included two components: one comparing active and passive 71 movements and one comparing reaching movements in two different workspaces. 72  80 We reported previously that the direction of maximal response (the "preferred direction"; PD) of 81 single neurons in area 2 is typically similar during active reaching movements and passive 82 perturbations of the hand (London and Miller 2013). However, this section shows that despite the 83 apparent similarity in directional tuning of individual neurons, information from a population of 84 neurons can be used to clearly distinguish active and passive movements, which is unexpected 85 given our previous results. We set out to explore whether this separation could be explained by 86 modeling the neural activity in area 2 purely in terms of behavioral variables. 87 In this experiment, the monkey performed a center-out reaching task to four targets. On half of 88 these reaching trials, the monkey's hand was bumped by the manipulandum during the center-89 hold period in one of the four target directions ( Figure 2 (PCA), which extracts the  129  sequential, orthogonal dimensions of highest covariance from the neural population activity, we  130  uncovered an unexpected linear boundary between active and passive movements within the first  131 three PCs (Figure 4). To quantify this separation, we trained a linear discriminant analysis (LDA) 132 model to classify the movement type in each of the four sessions. On average, for data not used 133 to train the LDA model, 84% of movements were correctly classified as active or passive, 134 indicating that despite the similarity of hand-movement coding by single neurons in area 2, the 135 neural population distinguishes the different movement types. 136 To explore potential causes of this separation in neural state space, we attempted to replicate it 137 using synthetic state spaces constructed from simplified neuron activity models that assume 138 different representations in area 2 (see Methods for modeling details). As a baseline, we fit a 139 purely hand kinematics model, in which neural activity was driven by the position and velocity 140 of hand movement. Classification rate of this simple model was slightly above chance but well 141 under the actual data, with only 65% correct movement classification ( Figure 4). Thinking that 142 the difference between active and passive movements may have been related to the differing 143 forces on the hand in the two conditions, we also fit a hand kinematics+force model, like the one 144 described in (Prud'homme and Kalaska 1994), where neural activity was driven by both the 145 velocity of the hand and the forces applied to it. Reaching only 64% correct classification, the 146 addition of force made no difference, indicating that area 2 contains some additional information 147 that allows it to distinguish between active and passive movements. 148 Finally, we fit a model that represented the Cartesian positions and velocities of the elbow as 149 well as the hand. Unlike the others, this model resulted in an average movement classification 150 accuracy of 82% ( Figure 4). This result provides a possible explanation of how area 2 separates 151 the two types of movement: while hand movement is similar, the movement of the whole arm 152 differs between active and passive conditions. While there might be additional explanations for 153 this separability, like an efference copy signal from motor areas specifically during active 154 movements (Bell 1981;London and Miller 2013;Nelson 1987), this finding suggests that the 155 parsimonious hand/elbow model is sufficient to explain it, unlike the hand-based models. 156

175
Given the importance of whole-arm kinematics for explaining neural population activity revealed 176 by the active vs. passive experiment, we conducted a second experiment to further examine how 177 single neurons represent reaching movements that are constrained to two disjoint areas in front of 178 the monkey, henceforth called workspaces ( Figure 5A). We found that whole-arm models and 179 hand-based models made different predictions of neural activity across these two conditions, 180 allowing us to further characterize the extent to which area 2 represented whole-arm movements 181 rather than just the movements of the hand. 182 In this experiment, we tested four different models of area 2 encoding, which we titled 183 "extrinsic", "egocentric", "muscle", and "hand-elbow", schematized in Figure 5D. In the hand-elbow model, the activity was related to the Cartesian kinematics (position and 193 velocity) of both the hand and the elbow markers. As in the active vs. passive experiment, we 194 aimed to test how well these models predicted features of neural activity during reaching. 195 However, a challenge in comparing these models of neural representation is that for the typical, 196 center-out reaching task in a small workspace, kinematic signals in the various coordinate frames 197 are highly correlated. Because a high correlation means that a linear transform can accurately 198 convert one coordinate frame into another, all four models make very similar predictions of 199 neural activity. To deal with this problem, we trained the monkeys to reach to randomly 200 generated targets presented in two different workspaces: one close to the body and contralateral 201 to the reaching hand, and one distant and ipsilateral. This had two important effects. First, the 202 random locations of the targets lessened the stereotopy of the movements, allowing for the 203 collection of more varied movement data than from the center-out paradigm. Second, the average 204 posture in each workspace was significantly different, such that while signals in the different 205 model coordinate frames were still correlated within each workspace, this correlation (and the 206 mapping between coordinate frames) changed significantly between workspaces. This forced the 207 models to make different predictions of neural activity across the two workspaces. Figure 5A  208 schematizes the two-workspace random target behavior, and Figure 5B shows an example raster 209 plot of neurons while a monkey performed the task. We recorded three sessions with each of 210 Monkeys C and H and two sessions with Monkey L. 211

227
Note that both the actual tuning curve and PD change between workspaces. Neither the egocentric nor 228 extrinsic models predict these changes, but the muscle and hand/elbow models can.

229
The simplest evaluation of these models is to compare how well they predicted actual neural 230 firing rates ( Figure 5C). To assess this, we used repeated k-fold cross-validation of a goodness-231 of-fit metric (see Methods for more details). Here, normal goodness-of-fit metrics like R 2 or 232 variance-accounted-for (VAF) are ill-suited to the Poisson-like statistics of neural activity; 233 instead, we used the likelihood-based pseudo-R 2 (Cameron and Windmeijer 1997;234 McFadden 1977). Like VAF, psuedo-R 2 has a maximum value of 1, but can also be negative for 235 models that do worse than predicting the mean firing rate during cross-validation. In general, the 236 values corresponding to what constitutes a good fit are lower for pR 2 than for either R 2 or VAF. 237 Of our four models, the whole-arm models out-performed the hand models ( Figure 6). Generally, 238 the hand-elbow model was the most predictive of actual neural firing rate across the two 239 workspaces, winning the great majority of the pairwise comparisons with the other models 240 ( Figure

261
We further tested our models on how well features (e.g., tuning curves and PDs) computed from 262 the model-predicted firing rates matched those of the actual data. Figure 5D shows the 263 directional tuning curves for an example neuron, along with the tuning curves predicted by each 264 model. We calculated the correlation between the predicted and the actual tuning curves in 265 different workspaces as a measure of their similarity. With this measure, the hand-elbow model 266 resulted in the best reconstruction of tuning curves, once again winning most of the pairwise 267 comparisons with the other models (Figure 7, S2). In this case, for only 37 neurons did either of 268 the hand-based models win any pairwise comparison against a whole-arm model (using = 269 0.05 and a Bonferroni correction to account for six total pairwise comparisons for each neuron) 270 For the other 251 cells, neither hand-based model could out-predict either whole-arm model. 271 Therefore, this experiment further suggests that area 2 of S1 encodes whole-arm rather than just 272 hand kinematics. 273

283
Of the 288 recorded neurons, 260 were significantly tuned to movement direction in both 284 workspaces. This fraction of tuned neurons is much higher than the fraction we found in the 285 active vs. passive task, likely for two reasons. First, unlike with the active vs. passive task, we 286 used the entirety of each trial in the two-workspace experiment to calculate PDs. Second, 287 monkeys tended to move faster during this random target task than during the center-out task, 288 which likely led to higher modulation depths and less measurement error, as predicted by 289 (Stevenson et al. 2011). Thus, in addition to the goodness-of-fit and tuning curve correlation 290 analyses, we were also able to examine the properties of the neural PDs in the two workspaces. 291 An interesting feature of this task is that for many neurons, the PD of movement changed 292 significantly between workspaces, exemplified by the difference between the vertical bars in the 293 leftmost panel of Figure 5D. Figure 8A shows the actual PD shifts for these neurons plotted 294 against the PD shifts predicted by each model. The large changes in PD, shown on the vertical 295 axes of the scatter plots are a clue that the extrinsic model does not explain neural activity 296 correctly; if it did, the preferred direction changes should have been insignificant (in principle,  297 zero), as shown by the generally small extrinsic model-predicted changes (second column of 298 Figure 8A). Additionally, and perhaps counterintuitively, the actual changes included both 299 clockwise and counter-clockwise rotations, so it is also unlikely that they arose from a rotation of 300 the extrinsic coordinate frame about the shoulder or the egocentric model. However, we found 301 that the whole-arm models did predict both clockwise and counter-clockwise PD changes. Based 302 on the circular VAF (cVAF) of the PD change prediction, Figure 8B shows that the hand-elbow 303 model once again out-predicted the other models, with hand/elbow having the highest average 304 cVAF (0.73), followed in order by muscle (0.63), extrinsic (0.55), and egocentric (0.38). We 305 made pairwise comparisons between model cVAF for each session. In every session but one, the 306 hand-based models lost pairwise comparisons to at least one of the whole-arm models, and in the 307 outlier session, no pairwise comparisons showed significance (again using = 0.05 and a 308 Bonferroni correction to account for six comparisons for each session). 309  activity. Even so, as shown through the main goodness-of-fit ( Figure 6, S1), tuning curve 336 correlation (Figure 7, S2), and cVAF ( Figure 8) analyses, both of these whole-arm models are 337 better than the hand-based models, supporting the fact that area 2 encodes whole arm kinematics. 338

340
In this study, we explored how somatosensory area 2 represents reaching movements using two 341 separate experiments. In the first experiment, we found that although single neuron directional 342 tuning is largely preserved between active and passive movements, the two types of movements 343 can be differentiated in the area 2 neural state space. While the classic hand-centric models 344 cannot explain this separation, we found that a model of neural activity built from the kinematics 345 of both the hand and elbow did. We further explored area 2 encoding of whole-arm kinematics in 346 a second experiment, in which the monkey reached to targets in a pair of workspaces that 347 changed the relationship between hand-based and whole-arm models. Across these two 348 conditions, features of neural activity, including the dynamics of neural discharge, tuning curve 349 shape and preferred direction, were better explained by a model incorporating the kinematics of 350 the whole arm than by the classic hand-based models. Altogether, these results suggest that area 351 2 of S1 represents movements in terms of not just the hand, but the whole arm. 352

353
A significant difference between the hand and whole-arm models is their number of parameters, 354 which make the whole-arm models more complex and expressible. There are two concerns with 355 testing models of differing complexity, the first dealing with model training and evaluation, and 356 the second with interpretation of the results. 357 In training and evaluating our models, we had to make sure that the complex models did not 358 overfit the data, resulting in artificially high performance on the training dataset but low 359 generalizability to new data. However, because we found through cross-validation that the more 360 complex models generalized to test data better than the simpler models, they were not 361 overfitting. Consequently, the hand-based models are clearly impoverished compared to the 362 whole-arm models. 363 The second concern is in interpreting what it means when the more complex models performed 364 better. One interpretation is that this is an obvious result; if the added degrees of freedom have 365 anything at all to do with area 2 neural activity, then the more complex models should perform 366 better. However, the choice of the two less complex hand-based models was not arbitrary. Both 367 are classic models within the canon of the cortical representation of limb state (Bosco et al. 2000;368 Prud'homme and Kalaska 1994). Additionally, as shown in the active/passive study, if the classic 369 extrinsic model is to be believed, then the separation between neural representations of active 370 and passive movement must have an explanation not rooted in the behavioral kinematics. Thus, 371 the point of this study was to determine whether either of these now classic models adequately 372 encapsulates what area 2 represents. Because the whole-arm models outperformed the classic 373 models, we conclude that the classic models are incomplete, and in the worst cases, potentially 374 incorrect or misleading. 375

376
One surprising result from the two-workspace experiment was that the hand/elbow model, based 377 on the addition of one arbitrary point on the proximal limb to the classic hand model, performed 378 as well or better than the muscle-based model. As proprioceptive signals originate in the 379 muscles, arising from muscle spindles and Golgi tendon organs, we expected to find that the 380 muscle model would outperform the other models. However, there are several potential reasons 381 why this was not so. The most important ones can be divided into two categories loosely tied to 382 1) errors in estimating the musclulotendon lengths, through motion tracking and musculoskeletal 383 modeling, and 2) the fidelity of the muscle model to the actual signals sent by the proprioceptors. 384 In the first category, the main issue is that of error propagation. The extra stages of analysis 385 required to compute musculotendon lengths (registering markers to a musculoskeletal model, 386 performing inverse kinematics to find joint angles, and using modeled moment arms to estimate 387 musculotendon lengths) introduce errors not present when simply using the positions of markers 388 on the arm. As a control, we ran the hand/elbow model through two of these extra steps by 389 computing the hand and elbow positions from the joint angles of the scaled model, estimated 390 from inverse kinematics. The results of this analysis showed that the performance of the 391 hand/elbow model with added noise dropped to that of the muscle model, indicating that there 392 are, in fact, errors introduced in even this portion of the processing chain. 393 The other potential source of error in this processing chain stems from the modeled moment 394 arms, which might not accurately reflect those of the actual muscles. In developing their 395 musculoskeletal model, Chan  have to record directly from gamma motor neurons or make assumptions of how gamma drive 417 changes over the course of reaching. In developing models of neural activity, one must carefully 418 consider the tradeoff between increased model complexity and the extra error introduced by 419 propagating through the additional requisite measurement and analysis steps. Given our data 420 obtained by measuring the kinematics of the arm with motion tracking, it seems that the 421 coordinate frame with which to best explain area 2 neural activity is simply the one with the 422 most information about the arm kinematics and the fewest steps in processing. However, this 423 does not rule out the idea that area 2 more nearly represents a different whole-arm model that 424 may be less abstracted from physiology, like musculotendon length or muscle spindle activity. 425

426
Because of their differing dimensionality, the signals from hand-based models and those from 427 whole-arm models do not have a one-to-one relationship: there are many different arm 428 configurations that result in a given hand position. Thus, a comparison between the hand-based 429 and whole-arm models is mainly a question of information content (do area 2 neurons have 430 information about more than just the hand?). In contrast, signals of the two whole-arm models do 431 have a one-to-one (albeit nonlinear) relationship to each other. Knowledge of the hand and elbow 432 position should completely determine the estimated musculotendon lengths, indicating that the 433 two models have the same informational content. As such, a comparison between the muscle 434 model and the hand/elbow is purely one of coordinate frame. While the interpretation for a 435 comparison of information content is straightforward, interpreting the results of a comparison 436 between coordinate frames is not. One major issue is that these comparisons only make sense 437 when using linear models to relate neural activity to behavior. This study shows that even after proprioceptive signals reach area 2, neural activity can still be 479 predicted well by a convergence of muscle-like signals, even though the signals have been 480 processed by several sensory areas along the way. One potential explanation for this is that at 481 each stage of processing, neurons simply spatially integrate information from many neurons of 482 the previous stage, progressively creating more complex response properties. This idea of 483 hierarchical processing was first used to explain how features like edge detection and orientation 484 tuning might develop within the visual system from spatial integration of the simpler 485 photoreceptor responses (Felleman and Van Essen 1991;Hubel and Wiesel 1959;1962 towards providing proprioceptive feedback remains relatively unexplored. At least one study did 504 use electrical stimulation in S1 for feedback during movement, using the stimulation to specify 505 target direction with respect to the evolving hand position (Dadarlat et al. 2015). In that study, 506 monkeys used the ICMS to reach to targets, even in the absence of visual feedback. However, 507 target-location information is very different from the information normally encoded by S1, and 508 the monkeys required several months to learn to use it. To our knowledge, no study has yet 509 shown a way to use ICMS to provide more biomimetic proprioceptive feedback during reaching. 510 Previously, our lab attempted to address this gap by stimulating a small number of electrodes in 511 area 2 based on neural activity recorded from them during normal reaching movements. In that 512 experiment, the monkey reported the direction of a mechanical bump to his arm that occurred 513 simultaneously with the ICMS. The ICMS biased one monkey's reports of the mechanical bump 514 direction toward the PDs of the stimulated electrodes. Key to this finding was the fact that any 515 bias in reporting actually decreased the reward rate, suggesting that the ICMS was 516 indistinguishable from the perception of the bump itself (Tomlinson and Miller 2016). 517 Unfortunately, the result could not be replicated in other monkeys; while the ICMS often biased 518 their reports, the direction of the bias could not be explained by the PDs of the stimulated 519 electrodes. One potential reason may be that the stimulation paradigm in those experiments was 520 derived from the classic, hand-based model and the assumption that area 2 represents active and 521 passive movements similarly. As this paper has shown, both of these assumptions have important 522 caveats. It is possible that a stimulation paradigm based on a whole-arm model may be more 523 successful, due to its greater accuracy at predicting neural PDs (Figure 8). It is also possible that 524 the stimulus model would need to include information about forces in addition to kinematics. 525 Regardless of the exact model, prospects for stimulating S1 to create natural proprioceptive 526 sensations would likely improve given a more accurate generative model of S1 activity. 527 In addition to developing better models for S1 activity, it will be important to consider the 528 implications of the difference between sensation for perception versus action. These two broad 529 purposes for sensation are thought to involve distinct pathways in both vision and touch 530 ( submovements that those guided by even a noisy visual signal, suggesting that monkeys used the 537 ICMS as a learned sensory substitute, rather than as a biomimetic replacement for 538 proprioception. As such, that study was also likely a cognitive one, engaging the perceptual 539 stream rather than the action stream of proprioception (see (Deroy and Auvray 2012; Elli et al.  540 2014) for discussion of the limits of sensory substitution). As we better characterize how S1 541 represents movements, we hope to develop a stimulation paradigm in which we can engage both 542 streams, to enable users of a BCI both to perceive their limb, and to respond rapidly to 543 movement perturbations. 544

545
This study began with an observation: the classic, hand-based cortical model of proprioception 546 could not explain the separability of active and passive movements we observed in the neural 547 state space within area 2. We found, however, that this feature could be explained by extending 548 the classic model to include the kinematics of the whole arm. In a second experiment, we found 549 that predictions of area 2 neural activity from such whole-arm models generalized to different 550 behavioral conditions better than those of the classic model. This suggests that even though our 551 perception of our arm is typically centered on the hand, area 2 of S1 still appears to represent 552 movement of the whole arm. 553

554
We would like to thank Brian London for initial discussions of the active vs. passive result and  555  Tucker Tomlinson, Christopher VerSteeg, and Joseph Sombeck for their help with training and  556 caring for the research animals. Additionally, we would like to thank them, along with Matt 557 Perich, Juan Gallego, Sara Solla, and the entire Miller Limb Lab for discussions and feedback 558 that greatly improved this work. 559 Northwestern University under protocol #IS00000367. 566

567
We recorded data from a monkey while it used a manipulandum to reach for targets presented on 568 a screen within a 20 cm x 20 cm workspace. After each successful reaching trial, the monkey 569 received a pulse of juice or water as a reward. We recorded the position of the handle using 570 encoders on the manipulandum joints. We also recorded the interaction forces between the 571 monkey's hand and the handle using a six-axis load cell mounted underneath the handle. 572 For the active vs. passive experiment, we had the monkey perform a classic center-out (CO) 573 reaching task, as described in (London and Miller 2013). Briefly, the monkey held in a target at 574 the center of the full workspace for a random amount of time, after which one of four outer 575 targets was presented. The trial ended in success once the monkey reached to the outer target. On 576 50% of the trials (deemed "passive" trials), during the center hold period, we used motors on the 577 manipulandum to deliver a 2 N perturbation to the monkey's hand in one of the four target 578 directions. After the bump, the monkey returned to the center target, after which the trial 579 proceeded like an active trial. From only the successful passive and active trials, we analyzed the 580 first 120 ms after movement onset. Movement onset was determined by looking for the peak in 581 handle acceleration either after the motor pulse (in the passive condition) or after 200 ms post-go 582 cue (in the active condition) and sweeping backwards in time until the acceleration was less than 583 10% of the peak. 584 For the two-workspace experiment, we partitioned the full workspace into four 10cm x 10cm 585 quadrants. Of these four quadrants, we chose the far ipsilateral one and the near contralateral one 586 in which to compare neural representations of movement. Before each trial, we chose one of the 587 two workspaces randomly, within which the monkey reached to a short sequence of targets 588 randomly positioned in the workspace. For this experiment, we only analyzed the portion of data 589 from the end of the center-hold period to the end of the trial. 590

591
Before each reaching experiment, we painted 10 markers on the outside of the monkey's arm, 592 marking bony landmarks and a few points in between, a la Chan and Moran 2006. Using a 593 custom motion tracking system built from a Microsoft Kinect, we recorded the 3D locations of 594 these markers with respect to the camera, synced in time to the other behavioral recordings. We 595 then aligned the Kinect-measured marker locations to the lab frame by aligning location of the 596 Kinect hand marker to the location of the handle in the manipulandum coordinate frame. Code 597 for motion tracking can be found at https://github.com/limblab/KinectTracking.git. 598

599
We registered the Kinect marker locations to a monkey arm musculoskeletal model in OpenSim 600 (SimTK), based on a model published by (Chan and Moran 2006

607
We implanted 100-electrode arrays (Blackrock Microsystems) into the arm representation of area 608 2 of S1 in these monkeys. For more details on surgical techniques, see (Weber et al. 2011). In 609 surgery, we roughly mapped S1 by recording from the cortical surface while manipulating the 610 arm and hand to localize their representations. To record neural data for our experiments, we 611 used a Cerebus recording system (Blackrock). This recording system sampled signals from each 612 of the 96 electrodes at 30 kHz. To conserve data storage space, the system detected spikes online 613 using a threshold set at -5 x signal RMS, and only wrote to disk a time stamp and the 1.6 ms 614 snippet of signal surrounding the threshold crossing. After data collection, we used Plexon 615 Offline Sorter to manually sort these snippets into putative single units, using features like 616 waveform shape and inter-spike interval. In addition to these recording sessions, we also 617 occasionally performed sensory mapping sessions to identify the neural receptive fields by 618 manipulating the monkey's arm while listening to neural activity. In all monkeys, we found a 619 roughly equal mix of cutaneous and deep (muscle) receptive fields, suggesting that we were 620 recording primarily from area 2 (Hyvärinen and We used a simple bootstrapping procedure to calculate PDs for each neuron. On each bootstrap 628 iteration, we randomly drew timepoints from the reaching data, making sure that the distribution 629 of movement directions was uniform to mitigate the effects of any potential bias. Then, as in 630 (Georgopoulos et al. 1982), we fit a cosine tuning function to the neural activity with respect to 631 the movement direction, using equations 1a-b. 632 ( ) = 0 + 1 * sin( ( )) + 2 * cos( ( )) (1 ) 633 where 635 = 2( 1 , 2 ) and = ( 1 2 + 2 2 ) 636 Here, ( ) is the average firing rate of neuron for a given time point , and ( ) is the 637 corresponding movement direction, which for the active/passive task was the target or bump 638 direction, and for the two-workspace experiment was the average movement direction over a 639 time bin. We took the circular mean of and mean of over all bootstrap iterations to 640 determine the preferred direction and the modulation depth respectively, for each neuron. 641 As the PD analysis is meaningless for neurons that don't have a preferred direction of movement, 642 we only analyzed the PDs of neurons that were significantly tuned. We assessed tuning through a 643 separate bootstrapping procedure, described in (Dekleva et al. 2018). Briefly, we randomly 644 sampled the timepoints from reaching data, again ensuring a uniform distribution of movement 645 directions, but this time also randomly shuffled the corresponding neural activity. We calculated 646 the for this shuffled data on each bootstrap iteration, thereby creating a null distribution of 647 modulation depths. We considered a neuron to be tuned if the true was greater than the 95 th 648 percentile of the null distribution. 649 7.5.2 Models of neural activity 650 For the active/passive analyses, we averaged behavioral variables and neural firing rates over the 651 120 ms period following movement onset in each trial. For the two-workspace analyses, both 652 behavioral variables and neural firing rate were averaged over 50 ms bins. We modeled neural 653 activity with respect to the behavior using Poisson generalized linear models (outline in 654 (Truccolo et al. 2005 We used repeated 5-fold cross-validation to evaluate our models of neural activity, given that the 683 models had different numbers of parameters, . On each repeat, we randomly split trials into five 684 groups (folds) and trained the models on four of them. We used these trained models to predict 685 neural firing rates (̂) in the fifth fold. We then compared the predicted firing rates from each 686 model to the actual firing rates in that test fold, using analyses described in the following 687 sections. This process (including random splitting) was repeated 20 times, resulting in n=100 688 sample size for each analysis result. Thus, if a more expressive model with more parameters 689 performs better than a simpler model, it would suggest that the extra parameters do provide 690 relevant information about the neural activity not accounted for by the simpler models. 691

Statistical tests and confidence intervals 692
To perform statistical tests on the output of repeated 5-fold cross-validation, we used a corrected 693 resampled t-test, outlined in (Ernst 2017) and (Nadeau and Bengio 2003). Here, sample mean 694 and variance are calculated as in a normal t-test, but a correction factor needs to be applied to the 695 standard error, depending on the nature of the cross-validation. Equation 3a-c shows a general 696 case of this correction for R repeats of K-fold cross-validation of some analysis result . 697 We then compare the t-statistic here ( ) to a t-distribution with × − 1 degrees of 701 freedom. Note that the correction applied is an extra term (i.e., The neural response to movement, whether active or passive, can be represented as a single 706 datapoint in a neural population space defined by the activity of all neurons in the relevant 120 707 ms period following movement onset. As described above, we estimated the separability of the 708 active and passive movements in neural space using 20x repeated, 5-fold cross-validation. We 709 did this in three steps. First, for each training set, we characterized the population response to 710 each trial by finding the first three modes of neural activity using principal component analysis 711 (PCA). We then projected the neural activity onto these three principal components and trained a 712 linear discriminant analysis (LDA) model to find the axis of maximal separation between active 713 and passive trials. Finally, we sequentially projected each test fold's neural data into the PC 714 space and then onto the LDA axis. This resulted in a scalar value for each trial, with the sign 715 indicating whether LDA classified the trial as active or passive. We took the average 716 classification accuracy of the test fold's data as the percent separability for the fold, giving us 717 100 total samples from the 20x5-fold cross-validation. By averaging these samples, we estimated 718 the overall neural separability of active and passive movements in a given session. 719

Estimating model-predicted separability 720
We also trained encoding models to predict neural firing from behavior (see equation 2a for 721 procedure) using three different models of neural activity: the extrinsic kinematics, extrinsic 722 kinematics+force, and hand/elbow kinematics models. See "Models of neural activity" section 723 for more details on the specific models. After training each model over the four training folds, 724 we estimated firing rates in both the training and the test folds. Subsequent analysis mirrored that 725 of the actual neural data: we found the three leading PCs and LDA axis of highest separability in 726 the training folds and then sequentially projected the test-fold data through the PC space and 727 onto the LDA axis. This resulted in 100 samples with which to estimate the model-predicted 728 separability of active and passive movements. 729

Neural space dimensionality reduction 730
To visualize the population neural activity for figures, we used a combination of LDA and PCA. 731 For the horizontal axis, we used LDA to find an axis in the three PCs of neural population space 732 along which active and passive trials were most separated. For the vertical axis, we projected all 733 activity onto the hyperplane orthogonal to the LDA axis and used PCA again to find the 734 remaining axis of highest variance. 735 7.5.5 Two-workspace analyses 736 These analyses examined how well models of neural activity could predict neural activity as the 737 monkey reached to targets in different workspaces. As such, we analyzed firing rate goodness-738 of-fit, along with how well the models could replicate the tuning curves and preferred directions 739 (PDs) of neurons. 740 7.5.5.1 Goodness-of-fit 741 We evaluated goodness-of-fit of these models for each neuron by using a pseudo-R 2 ( 2 ) 742 metric. We used a formulation of pseudo-R 2 based on a comparison between the deviance of the 743 full model and the deviance of a "null" model, i.e., a model that only predicts the overall mean 744 firing rate (Cameron and Windmeijer 1997; This pR 2 metric ranges from −∞ to 1, with a value of 1 corresponding to a perfectly fit model 751 and a value of 0 corresponding to a model that only fits as well as the "null" model. In contrast  752 with the general intuition for regular R 2 , a pR 2 of ~0.2 is considered a "good" fit (McFadden 753 1977). 754 7.5.5.2 Tuning curves 755 We binned the trajectory into 16 bins, each 22.5 degrees wide, based on the mean direction 756 across 50 ms of hand motion. For each directional bin, we calculated the sample mean and 95% 757 confidence interval of the mean. In figures, we plotted this mean firing rate against the center-758 point of the bin. 759

Preferred direction shift 760
We calculated PDs for each neuron in each workspace and found the predicted change in PD 761 from the contralateral workspace to the ipsilateral workspace, given each model. We compared 762 these changes to those observed for each neuron. The values of these PD shifts are shown in 763 Figure 8 for all neurons tuned to movements in both workspaces, averaged over all 100 test 764 folds. 765 We computed a variance-accounted-for (VAF) metric, here called the "circular VAF" (cVAF) 766 for each neuron ( ) in each fold as: 767 = ( , −̂, ) (5) 768 As the cVAF metric is essentially the inner product of unit vectors with direction , and 769 ̂, , it accounts for the circular domain of the PD shifts. Like regular VAF, the cVAF has a 770 maximum value of 1 when , and ̂, are the same, and decreases in proportion to the 771 squared difference between , and ̂, . We took the average cVAF over all neurons as 772 the cVAF for the fold. In total, given the 20 repeats of 5-fold cross-validation, this gave us 100-773 samples of the cVAF for each model in a given session. 774