Stability of motor representations after paralysis.

Neural plasticity allows us to learn skills and incorporate new experiences. What happens when our lived experiences fundamentally change, such as after a severe injury? To address this question, we analyzed intracortical population activity in the posterior parietal cortex (PPC) of a tetraplegic adult as she controlled a virtual hand through a brain-computer interface (BCI). By attempting to move her fingers, she could accurately drive the corresponding virtual fingers. Neural activity during finger movements exhibited robust representational structure similar to fMRI recordings of able-bodied individuals' motor cortex, which is known to reflect able-bodied usage patterns. The finger representational structure was consistent throughout multiple sessions, even though the structure contributed to BCI decoding errors. Within individual BCI movements, the representational structure was dynamic, first resembling muscle activation patterns and then resembling the anticipated sensory consequences. Our results reveal that motor representations in PPC reflect able-bodied motor usage patterns even after paralysis, and BCIs can re-engage these stable representations to restore lost motor functions.

Introduction

A central question in neuroscience is how experience affects the nervous system. Studies of this phenomenon, plasticity, were pioneered by Hubel and Wiesel, who found that temporary visual occlusion in kittens can induce lifelong reorganization of the visual cortex (Hubel and Wiesel, 1970). Their results demonstrated that the developing brain, rather than being genetically preprogrammed, is surprisingly malleable to external inputs.

Subsequent studies showed that other brain regions are also plastic during early development, but it is unclear how plastic the nervous system remains into adulthood. Visual occlusion in adult cats does not reorganize the visual cortex, and lesion studies of the adult visual cortex have arrived at competing conclusions of reorganization and stability (Smirnakis et al., 2005; Gilbert and Wiesel, 1992; Keck et al., 2008; Baseler et al., 2011). A similar discussion continues regarding the primary somatosensory cortex (S1). Classical studies posited that amputation and spinal cord injury modify the topography of body parts in S1, with intact body parts taking over cortical areas originally dedicated to the amputated part (Merzenich et al., 1984; Qi et al., 2000; Pons et al., 1991; Jain et al., 2008). However, recent human neuroimaging studies and sensory BCI studies have challenged the extent of this remapping, arguing that sensory topographies largely persist even after complete sensory loss (Makin and Bensmaia, 2017; Kikkert et al., 2021; Flesher et al., 2016; Armenta Salas et al., 2018). Thus, the level of plasticity in the adult nervous system is still an ongoing investigation.

Understanding plasticity is necessary to develop brain–computer interfaces (BCIs) that can restore sensorimotor function to paralyzed individuals (Orsborn et al., 2014). First, paralysis disrupts movement and blocks somatosensory inputs to motor areas, which could cause neural reorganization (Pons et al., 1991; Jain et al., 2008; Kambi et al., 2014). Second, BCIs bypass supporting cortical, subcortical, and spinal circuits, fundamentally altering how the cortex affects movement. Do these changes require paralyzed BCI users to learn fundamentally new motor skills (Sadtler et al., 2014), or do paralyzed participants use a preserved, pre-injury motor repertoire (Hwang et al., 2013)? Several paralyzed participants have been able to control BCI cursors by attempting arm or hand movements (Hochberg et al., 2006; Hochberg et al., 2012; Collinger et al., 2013; Gilja et al., 2015; Bouton et al., 2016; Ajiboye et al., 2017; Brandman et al., 2018), hinting that motor representations could remain stable after paralysis. However, the nervous system’s capacity for reorganization (Pons et al., 1991; Jain et al., 2008; Kikkert et al., 2021; Kambi et al., 2014) still leaves many BCI studies speculating whether their findings in tetraplegic individuals also generalize to able-bodied individuals (Flesher et al., 2016; Armenta Salas et al., 2018; Stavisky et al., 2019; Willett et al., 2020; Fifer et al., 2022). A direct comparison, between BCI control and able-bodied neural control of movement, would help address questions about generalization and plasticity.

Temporal dynamics provide another lens to investigate neural organization and its changes after paralysis. Temporal signatures can improve BCI classification (Willett et al., 2021) or provide a baseline for motor adaptation studies (Stavisky et al., 2017; Vyas et al., 2018). Notably, motor cortex (MC) activity exhibits quasi-oscillatory dynamics during arm reaching (Churchland et al., 2012). More generally, the temporal structure can depend on the movement type (Suresh et al., 2020) and the recorded brain region (Schaffelhofer and Scherberger, 2016). In this study, we recorded from the posterior parietal cortex (PPC), which is thought to compute an internal forward model for sensorimotor control (Mulliken et al., 2008; Wolpert et al., 1998; Desmurget and Grafton, 2000; Li et al., 2022). A forward-model overcomes inherent sensory delays to enable fast control by predicting the upcoming states. If PPC activity resembles a forward model after paralysis, this would suggest that even the temporal details of movement are preserved after injury.

Here, we investigate the neural representational structure of BCI finger movements in a tetraplegic participant. In able-bodied individuals, the cortical organization of finger representations follows the natural statistics of movements (Ejaz et al., 2015; Lillicrap and Scott, 2013). In a BCI task, the experimenter can instruct movement patterns unrelated to biomechanics or before-injury motifs. In this study, we tested whether the neural representational structure of BCI finger movements by a tetraplegic individual matches that of able-bodied individuals performing similar, overt movements, or whether the structure follows the task’s optimal representational structure (Bonnasse-Gahot and Nadal, 2008). If the BCI finger organization matches that of able-bodied movement, participants would likely be able to activate pre-injury motor representations, indicating that motor representations were preserved after paralysis.

We report that the neural representational structure of BCI finger movements in a tetraplegic individual matches that of able-bodied individuals. This match was stable across sessions, even though the measured representational structure contributed to errors in the BCI task. Furthermore, the neural representational dynamics matched the temporal profile expected of a forward model, first resembling muscle activation patterns and then resembling expected sensory outcomes. Our results suggest that adult motor representations in PPC remain even after years without use.

Results

Intracortical recordings during finger flexion

We recorded single and multineuron activity (95.8 ± standard deviation [SD] 6.7 neurons per session over 10 sessions) from participant NS while she attempted to move individual fingers of the right hand. We recorded from a microelectrode array implanted in the left (contralateral) PPC at the junction of the postcentral and intraparietal sulci (PC-IP, Figure 1—figure supplement 1). This region is thought to specialize in the planning and monitoring of grasping movements (Andersen et al., 2019; Orban and Caruana, 2014; Gallivan and Culham, 2015; Klaes et al., 2015).

Figure 1—figure supplement 1.

Multielectrode array implant location.

Figure 1—figure supplement 1 and legend text have been reproduced from Figure S1 of Aflalo et al., 2020. The original image and legend text are published under the terms of the Creative Commons Attribution-NonCommercial 4.0 International license (CC BY-NC 4.0) https://creativecommons.org/licenses/by-nc/4.0/. We used fuctional magnetic resonance imaging (fMRI) to identify cortical regions involved in imagined reaching and grasping actions. The participant performed two complementary tasks to ensure activation was robust across paradigms. (a) Event-related task design. Following an intertrial interval, the subject was cued with a specific imagined movement (precision grasp, power grasp, or reach without hand shaping). Following the cue, a cylindrical object was displayed. If the object was intact, the subject imagined performing the cued movement. If the object was broken, the subject withheld movement. (b) Block task design. Eight blocks were presented for 30 s per run. During the first 15 s of each block, common objects were presented every 3 s in varying spatial locations. Before each run, the subject was instructed to either imagine pointing at, imagine reaching and grasping, or look naturally at the object. During the last 15 s of each block, scrambled images were presented and the subject was instructed to guess the identity of the object. (c) Statistical parametric map showing voxels with significant activity for grasping (‘Go’ vs. ‘No-Go’) (p < 0.01, FDR-corrected), based on task (a). Array location and cortical landmarks are depicted in the legend. (d) Statistical parametric map showing voxels with significant activation (p < 0.01, FDR-corrected) for grasping versus looking, based on task (b).

Each recording session started with an initial calibration task (Figure 1—figure supplement 2, Methods). On each trial, we displayed a text cue (e.g., ‘T’ for thumb) on a computer screen, and the participant immediately attempted to flex the corresponding finger, as though pressing a key on a keyboard. Because participant NS previously suffered a C3–C4 spinal cord injury resulting in tetraplegia (AIS-A), her movement attempts did not generate overt motion. Instead, participant NS attempted to move her fingers as though she was not paralyzed.

Figure 1—figure supplement 2.

Calibration task.

Task structure, single trial. Each trial consisted of an intertrial interval (ITI) and a reaction-time Go phase. During the Go phase, green text specified which finger to flex. All letters were overlaid in gray to minimize visual differences between ITI and Go phases. T = thumb, I = index, M = middle, R = ring, P = pinky, X = no movement.

These attempted movements resulted in distinct neural activity patterns across the electrode array. To enable BCI control, we trained a linear classifier (Methods) to identify finger movements from neural firing rates. The participant subsequently performed several rounds of a similar finger flexion task, except that (1) the trained classifier now provided text feedback of its predicted finger and (2) the task randomized the visual cue location (Figure 1a and Methods). We repeated this online-control finger flexion task over multiple sessions (408 ± SD 40.8 trials/session over 10 sessions) and used this data for our offline analyses. Participant NS also performed a control task, identical in structure except that she attended to cues without performing the corresponding movements.

Figure 1.

Robust brain–computer interface (BCI) control of individual fingers.

(a) Main finger flexion task. When a letter was cued by the red crosshair, the participant looked at the cue and immediately attempted to flex the corresponding finger of her right (contralateral) hand. We included a null condition ‘X’, during which the participant looked at the target but did not move her fingers. Visual feedback indicated the decoded finger 1.5 s after cue presentation. To randomize the gaze location, cues were located on a grid (three rows, four columns) in a pseudorandom order. The red crosshair was jittered to minimize visual occlusion. (b) Confusion matrix showing robust BCI finger control (86% overall accuracy, 4016 trials aggregated over 10 sessions). Each entry (i, j) in the matrix corresponds to the ratio of movement i trials that were classified as movement j. (c–f) Mean firing rates for four example neurons, color-coded by attempted finger movement. Shaded areas indicate 95% confidence intervals (across trials of one session). Gaussian smoothing kernel (50 ms standard deviation [SD]).

Figure 1—figure supplement 1 and legend text have been reproduced from Figure S1 of Aflalo et al., 2020. The original image and legend text are published under the terms of the Creative Commons Attribution-NonCommercial 4.0 International license (CC BY-NC 4.0) https://creativecommons.org/licenses/by-nc/4.0/. We used fuctional magnetic resonance imaging (fMRI) to identify cortical regions involved in imagined reaching and grasping actions. The participant performed two complementary tasks to ensure activation was robust across paradigms. (a) Event-related task design. Following an intertrial interval, the subject was cued with a specific imagined movement (precision grasp, power grasp, or reach without hand shaping). Following the cue, a cylindrical object was displayed. If the object was intact, the subject imagined performing the cued movement. If the object was broken, the subject withheld movement. (b) Block task design. Eight blocks were presented for 30 s per run. During the first 15 s of each block, common objects were presented every 3 s in varying spatial locations. Before each run, the subject was instructed to either imagine pointing at, imagine reaching and grasping, or look naturally at the object. During the last 15 s of each block, scrambled images were presented and the subject was instructed to guess the identity of the object. (c) Statistical parametric map showing voxels with significant activity for grasping (‘Go’ vs. ‘No-Go’) (p < 0.01, FDR-corrected), based on task (a). Array location and cortical landmarks are depicted in the legend. (d) Statistical parametric map showing voxels with significant activation (p < 0.01, FDR-corrected) for grasping versus looking, based on task (b).

Task structure, single trial. Each trial consisted of an intertrial interval (ITI) and a reaction-time Go phase. During the Go phase, green text specified which finger to flex. All letters were overlaid in gray to minimize visual differences between ITI and Go phases. T = thumb, I = index, M = middle, R = ring, P = pinky, X = no movement.

Finger classification accuracy in the main task across 10 sessions. Session 1 accuracy excludes No-Go (X) trials (see Methods).

(a) Confusion matrix of offline finger classification, cross-validated within single sessions. 4080 trials of the main task aggregated over 10 sessions. (b) Confusion matrix of offline finger classification, cross-validated within single sessions. Five hundred and thirty trials of the calibration task aggregated over nine sessions. T = thumb, I = index, M = middle, R = ring, P = pinky, X = no movement. Each entry (i, j) in the matrix corresponds to the ratio of movement i trials that were classified as movement j.

All five fingers of the right (contralateral) hand were encoded within the population during movement execution. (a) Percentage of the population tuned significantly (p < 0.05, FDR-corrected) to flexion of each finger. Positive percentages indicate neurons that increased firing rate during finger movement and negative percentages indicates neurons that decreased firing rate. Error bars indicate a 95% bootstrap confidence interval. (b) Percentage of the population tuned best to each finger. (c) Cumulative distribution function of the population’s tuning significance p values. (d) Histogram of d′ (discriminability index) values across neurons. (a–d) Neurons were pooled across sessions.

(a) Linear regression could not decode target location [x, y] coordinates from the neural activity during the attempted-movement period. Violin plot shows that cross-validated regression r2 values are close to 0 across sessions, with each circle marking a single session. (b) A linear classifier (diagonal linear discriminant analysis [LDA]) could not classify the gaze location from neural activity during the attempted-movement period. Confusion matrix depicts cross-validated classifications of cue location. (c) Cross-validated classification accuracy for main and control tasks: a linear classifier (diagonal LDA) could not classify finger movements from neural activity during passive observation (orange) of the finger press task. Sliding bin width: 200 ms. The shaded region indicates ± standard error of the mean (SEM) (six sessions passive viewing, 10 sessions attempted flexion).

Accurately decoding fingers from PPC single-neuron activity

High classification accuracy during online control (86% ± SD 4% over 10 sessions; chance = 17%) (Figure 1b and Figure 1—figure supplement 3) and offline cross-validated classification (92% ± SD 2%; Figure 1—figure supplement 4a) demonstrated that the finger representations were reliable and linearly separable. During the calibration task, cross-validated classification was similarly robust (accuracy = 96% ± SD 3%; Figure 1—figure supplement 4b). These finger representations were robust across contexts and could be used in a range of environments, including to move the hand of a virtual reality avatar (Video 1).

Figure 1—figure supplement 3.

Brain–computer interface (BCI) classification accuracy across sessions.

Finger classification accuracy in the main task across 10 sessions. Session 1 accuracy excludes No-Go (X) trials (see Methods).

Figure 1—figure supplement 4.

Robust cross-validated finger classification during main and calibration tasks.

(a) Confusion matrix of offline finger classification, cross-validated within single sessions. 4080 trials of the main task aggregated over 10 sessions. (b) Confusion matrix of offline finger classification, cross-validated within single sessions. Five hundred and thirty trials of the calibration task aggregated over nine sessions. T = thumb, I = index, M = middle, R = ring, P = pinky, X = no movement. Each entry (i, j) in the matrix corresponds to the ratio of movement i trials that were classified as movement j.

Video 1.

Example brain–computer interface (BCI) control of a virtual reality hand.

Using a BCI, participant NS controls the individual fingers of a virtual reality hand. She views a virtual hand, table, and cues through an Oculus headset. Similar to the main finger movement task, she acquires green jewels by pressing the corresponding finger and avoids red amethysts by resting. Green jewels disappear when the correct finger is classified (or at the start of the next trial, if incorrectly classified). The screen copies the view that participant NS sees through the Oculus headset.

At the single-neuron level, most (89%) neurons were significantly tuned to individual finger-press movements (significance threshold: p < 0.05, FDR-corrected) (Figure 1—figure supplement 5). Figure 1c–f show the firing rates of example neurons, which were tuned to one or more fingers and change tuning profiles over the course of each movement.

Figure 1—figure supplement 5.

Single-neuron encoding of individual fingers.

All five fingers of the right (contralateral) hand were encoded within the population during movement execution. (a) Percentage of the population tuned significantly (p < 0.05, FDR-corrected) to flexion of each finger. Positive percentages indicate neurons that increased firing rate during finger movement and negative percentages indicates neurons that decreased firing rate. Error bars indicate a 95% bootstrap confidence interval. (b) Percentage of the population tuned best to each finger. (c) Cumulative distribution function of the population’s tuning significance p values. (d) Histogram of d′ (discriminability index) values across neurons. (a–d) Neurons were pooled across sessions.

To confirm that the observed neural responses could not be explained by visual confounds, we verified that we could not discriminate between fingers during the control task (Figure 1—figure supplement 6). Furthermore, we could not decode the gaze location during the finger classification time window in the standard online-control task (Figure 1—figure supplement 6). Thus, reliable finger representations emerged from the participant’s movement attempts.

Figure 1—figure supplement 6.

Gaze location did not affect finger decoding during the attempted-movement period.

(a) Linear regression could not decode target location [x, y] coordinates from the neural activity during the attempted-movement period. Violin plot shows that cross-validated regression r2 values are close to 0 across sessions, with each circle marking a single session. (b) A linear classifier (diagonal linear discriminant analysis [LDA]) could not classify the gaze location from neural activity during the attempted-movement period. Confusion matrix depicts cross-validated classifications of cue location. (c) Cross-validated classification accuracy for main and control tasks: a linear classifier (diagonal LDA) could not classify finger movements from neural activity during passive observation (orange) of the finger press task. Sliding bin width: 200 ms. The shaded region indicates ± standard error of the mean (SEM) (six sessions passive viewing, 10 sessions attempted flexion).

Finger representational structure matches the structure of able-bodied individuals

Having discovered that PC-IP neurons modulate selectively for finger movements, we next investigated how these neural representations were functionally organized and how this structure related to pre-injury movements. Here, we turned to the framework of representational similarity analysis (RSA) (Kriegeskorte et al., 2008a; Diedrichsen and Kriegeskorte, 2017). RSA quantifies neural representational structure by the pairwise distances between each finger’s neural activity patterns (Figure 2a). These pairwise distances form the representational dissimilarity matrix (RDM), a summary of the representational structure. Importantly, these distances are independent of the original feature types (e.g., electrode or voxel measurements), allowing us to compare finger organizations across subjects and across recording modalities (Kriegeskorte et al., 2008b).

Figure 2.

Representational structure during brain–computer interface (BCI) finger control matches the structure of able-bodied individuals.

(a) To calculate the representational dissimilarity matrix (RDM), a vector of firing rates was constructed for each trial. Repetitions were collected for each condition. Then, pairwise distances were estimated between conditions using a cross-validated dissimilarity metric. This process was repeated to generate an RDM for each session. We drop the No-Go condition (X) here to match previous finger studies (Ejaz et al., 2015; Kikkert et al., 2016). (b) Representational structure hypothesized by the preserved-representation hypothesis: average RDM for 36 able-bodied individuals performing a finger-press task. RDMs were measured at the junction of the postcentral and intraparietal sulci (PC-IP) using fMRI (Ejaz et al., 2015; Kieliba et al., 2021). Max-scaled to [0, 1]. (c) Representational structure hypothesized by the despecialization and task-optimal hypotheses: pairwise-equidistant RDM. Max-scaled to [0, 1]. (d) Finger representational structure measured in tetraplegic participant NS: cross-validated Mahalanobis distances (Methods) between neural activity patterns, averaged across 10 recording sessions. Max-scaled to [0, 1]. (e) Intuitive visualization of the distances in (d) using multidimensional scaling (MDS). Ellipses show mean ± standard deviation (SD) (10 sessions) after Generalized Procrustes alignment (without scaling) across sessions. (f) Measured RDMs (d) match the able-bodied PC-IP fMRI RDM (b) better than they match the task-optimal, unstructured model (c), as measured by the whitened unbiased cosine similarity (Diedrichsen et al., 2021) (WUC) (Methods). Mean differences were significant (able-bodied vs. unstructured, p = 5.7 × 10–5; two-tailed t-test, 1000 bootstrap samples over 10 sessions). Violin plot: solid horizontal lines indicate the median WUC over bootstrap samples, and dotted lines indicate the first and third quartiles. Noise ceiling: gray region estimates the best possible model fit (Methods). Asterisks denote a significant difference at ***p < 0.001. For convenience, a similar figure using a correlation-based similarity metric is shown in Figure 2—figure supplement 3.

fMRI representational dissimilarity matrices (RDMs) for three individual subjects and the group mean (N = 29). Intuitive visualization of distances using multidimensional scaling (MDS) and Generalized Procrustes alignment (without scaling); ellipses show mean ± standard deviation (SD) across subjects. Regions of interest (ROIs): motor cortex (MC, top row) and the junction of the postcentral and intraparietal sulci (PC-IP, bottom row).

RDMs across all sessions, using the cross-validated Mahalanobis distance. ‘Average’ RDM matches Figure 2d.

(a) Representational dissimilarity matrices (RDMs) calculated with an alternative dissimilarity metric: cross-validated Poisson KL-divergence (Schütt et al., 2019). Units: nats/neuron. Related to Figure 2—figure supplement 2a and Figure 2d. (b) Fit between measured RDMs and motor-intact BOLD data using alternative metrics. Distance metric: cross-validated Poisson KL-divergence. Similarity metric: whitened RDM Pearson correlation (Diedrichsen et al., 2021). Asterisks denote significant differences at ***p < 0.001.Similar to Figure 2f. (c) Representational dynamics calculated with an alternative dissimilarity metric: cross-validated Poisson KL-divergence. Similar to Figure 4e.

Gardner—Altman estimation plot (Ho et al., 2019) of the WUC similarity between same-region of interest (ROI) pairs of RDMs (N = 630 pairs between 36 subjects). Each circle on the swarm plot (left) marks the similarity for a pair of subjects. Horizontal black lines mark the mean pairwise similarity within each ROI. The curve (right) indicates the resampled (N = 5000) distribution of the effect size between ROIs, as measured by Cohen’s d. Cohen’s d of PC-IP minus MC: −2.1 (95% CI: [−2.22, −1.99]).

(a) Measured RDMs match the able-bodied MC fMRI RDM better than they match the able-bodied PC-IP fMRI RDM (p = 1.9 × 10–6; two-tailed t-test, 1000 bootstrap samples over 10 sessions), as measured by the whitened unbiased cosine similarity (Diedrichsen et al., 2021) (WUC) (Methods). Violin plot: solid horizontal lines indicate the median WUC over bootstrap samples, and dotted lines indicate the first and third quartiles. Noise ceiling: gray region estimates the best possible model fit (Methods). Asterisks denote significant differences at ***p < 0.001. Similar to Figure 2f. (b) Paired Gardner–Altman estimation plot (Ho et al., 2019) of the similarity (WUC) between participant NS (average RDM across sessions) and individual MC and PC-IP RDMs from able-bodied fMRI participants. The slopegraph’s connected points (left) show each fMRI participant’s (N = 36) MC and PC-IP similarities with participant NS’s mean finger RDM. Mean difference between PC-IP and MC similarities (right) presented as Cohen’s d (N = 5000 bootstrap samples).

Figure 2—figure supplement 3.

Representational structure during brain–computer interface (BCI) finger control matches the structure of able-bodied individuals when using alternative analysis parameters.

(a) Representational dissimilarity matrices (RDMs) calculated with an alternative dissimilarity metric: cross-validated Poisson KL-divergence (Schütt et al., 2019). Units: nats/neuron. Related to Figure 2—figure supplement 2a and Figure 2d. (b) Fit between measured RDMs and motor-intact BOLD data using alternative metrics. Distance metric: cross-validated Poisson KL-divergence. Similarity metric: whitened RDM Pearson correlation (Diedrichsen et al., 2021). Asterisks denote significant differences at ***p < 0.001.Similar to Figure 2f. (c) Representational dynamics calculated with an alternative dissimilarity metric: cross-validated Poisson KL-divergence. Similar to Figure 4e.

Figure 2—figure supplement 2.

Individual representational dissimilarity matrices (RDMs) for each session.

RDMs across all sessions, using the cross-validated Mahalanobis distance. ‘Average’ RDM matches Figure 2d.

Figure 4.

Representational dynamics analysis (RDA) dissociates neural processes over time.

(a) RDA performs representational similarity analysis (RSA) in a sliding window across time. Here, we model the measured representational structure as a non-negative linear combination of component model representational dissimilarity matrices (RDMs). (b–d) Hypothesized explanatory component RDMs: usage, muscle, and somatotopy (Ejaz et al., 2015). Max-scaled to [0, 1]. (e) RDA of the measured RDM over time shows an early fit to the muscle model and a late fit to the somatotopy model. Confidence intervals indicate ± standard error of the mean (SEM) bootstrapped across 10 sessions. Gray shaded region indicates the approximate onset time of the saccade to cue (interquartile range across trials). Difference in model start time (170 ms, Methods) was significant (p = 0.002, two-sided Wilcoxon signed-rank test). RDM snapshots (bottom, each max-scaled to [0, 1]) intuitively visualize the change in representational structure over time from muscle-like to somatotopic.

Violin plot of WUC similarity between the measured RDM (N = 1000 bootstrap samples over 10 sessions) and the corresponding model combination. Violin plot: solid horizontal lines indicate the mean WUC over bootstrap samples, and dotted lines indicate the first and third quartiles. Horizontal lines (above) indicate significance groups, where the circle-indicated model is significant over the vertical-tick-indicated models (two-tailed t-test, q < 0.01, FDR-corrected for 28 model-pair comparisons). For example, the muscle + somatotopy combined model is significant over the individual muscle, hand usage, somatotopy, combined muscle + hand usage, and pairwise-equidistant/unstructured (null) models.

When linear modeling within single sessions, the muscle model (blue) consistently preceded the somatotopy model (orange). Time difference: 170 ± 66 ms (SD across sessions) (p = 0.002, two-sided Wilcoxon signed-rank test). Line styles indicate session. Related to Figure 4e.

(a) Representational dynamics analysis shows an early fit to the hand-usage model and a late fit to the somatotopy model. Confidence intervals indicate ± standard error of the mean (SEM) across sessions. Related to Figure 4e. (b) Representational dynamics analysis shows a consistent delay between models during the calibration task. Note: The absolute timing differs from the main task because the calibration task does not require an initial saccade to read the cue. Related to Figure 4e.

(a) Histogram of L-ratio, a spike-sorting cluster metric. Threshold for well-isolated units: 33% quantile (Lratio < 10−1.1). (b) Representational dissimilarity matrices (RDMs) calculated only using well-isolated units, using the cross-validated Mahalanobis distance. Similar to Figure 2d and Figure 2—figure supplement 2a. (c) Whitened unbiased similarity (WUC) between measured (b) RDMs (using only well-isolated units) and model predictions (Figure 2b, c), showing that the measured RDMs match the able-bodied fMRI RDM significantly better than they match the unstructured model (p = 3.1 × 10–10, two-tailed t-test) and the SPLa fMRI RDM (p = 1.7 × 10–8). Error bars indicate ± standard error of the mean (SEM). Noise ceiling: gray region estimates the best possible model fit (Methods). Gray downward-semicircle indicates that the noise ceiling is significantly higher (p < 0.001) than the fit of the SPLa fMRI RDM and the unstructured model. Asterisks denote significant differences at ***p < 0.001. Similar to Figure 2f. (d) Representational dynamics analysis, repeated using only well-isolated units, shows an early fit to the muscle model and a late fit to the somatotopy model. Confidence intervals indicate ± SEM across sessions. Similar to Figure 4e.

We used RSA to test three hypotheses: (1) the BCI finger representational structure could match that of able-bodied individuals (Ejaz et al., 2015; Kieliba et al., 2021; Figure 2b and Figure 2—figure supplement 1), which would imply that motor representations did not reorganize after paralysis. This hypothesis would be consistent with recent functional magnetic resonance imaging (fMRI) studies of amputees, which showed that sensorimotor cortex representations of phantom limb finger movements match the same organization found in able-bodied individuals (Kikkert et al., 2016; Wesselink et al., 2019). We note that our able-bodied model was recorded from human PC-IP using fMRI, which measures fundamentally different features (millimeter-scale blood oxygenation) than microelectrode arrays (sparse sampling of single neurons). Another possibility is that (2) the participant’s pre-injury motor representations had despecialized after paralysis, such that finger activity patterns are unstructured and pairwise independent (Figure 2c). However, this hypothesis would be inconsistent with results from fMRI studies of amputees’ sensorimotor cortex (Kikkert et al., 2016; Wesselink et al., 2019). Lastly, (3) the finger movement representational structure might optimize for the statistics of the task (Lillicrap and Scott, 2013; Clancy et al., 2014). Our BCI task, as well as previous experiments with participant NS, involved no correlation between individual fingers, so the optimal structure would represent each finger independently to minimize confusion between fingers. In other words, the task-statistics hypothesis (3) would predict that, with BCI usage, the representational structure would converge toward the task-optimal, unstructured representational structure (Figure 2c).

Figure 2—figure supplement 1.

fMRI representational structure for finger movements, from Kieliba et al., 2021.

fMRI representational dissimilarity matrices (RDMs) for three individual subjects and the group mean (N = 29). Intuitive visualization of distances using multidimensional scaling (MDS) and Generalized Procrustes alignment (without scaling); ellipses show mean ± standard deviation (SD) across subjects. Regions of interest (ROIs): motor cortex (MC, top row) and the junction of the postcentral and intraparietal sulci (PC-IP, bottom row).

Does the finger representational structure in a tetraplegic individual match that of able-bodied individuals? We quantified the finger representational structure by measuring the cross-validated Mahalanobis distance (Methods) between each finger pair, using the firing rates from the same time window used for BCI control. The resulting RDMs are shown in Figure 2d (average across sessions) and Figure 2—figure supplement 2 (all sessions). For visual intuition, we also projected the representational structure to two dimensions in Figure 2e, which shows that the thumb is distinct, while the middle, ring, and pinky fingers are close in neural space. We then compared the measured RDMs against the able-bodied fMRI and unstructured models, using the whitened unbiased RDM cosine similarity (WUC) (Diedrichsen et al., 2021). The measured representational structure matched the able-bodied representational structure significantly over the unstructured model (p = 5.7 × 10–5, two-tailed t-test) (Figure 2f), ruling out the despecialization hypothesis (2). Our findings were robust to different choices of distance and model-similarity metrics (Figure 2—figure supplement 3).

We note that we constructed the able-bodied fMRI model from the mean of PC-IP fMRI RDMs across multiple able-bodied participants (N = 29). When compared among the RDM distribution of individual able-bodied participants, participant NS’s average PC-IP RDM was statistically typical (permutation shuffle test, p = 0.55), in part because PC-IP fMRI RDMs were relatively variable across able-bodied participants (Figure 2—figure supplement 4).

Figure 2—figure supplement 4.

fMRI finger representational dissimilarity matrices (RDMs) are more consistent across able-bodied participants in motor cortex (MC) than in the junction of postcentral and intraparietal sulci (PC-IP).

Gardner—Altman estimation plot (Ho et al., 2019) of the WUC similarity between same-region of interest (ROI) pairs of RDMs (N = 630 pairs between 36 subjects). Each circle on the swarm plot (left) marks the similarity for a pair of subjects. Horizontal black lines mark the mean pairwise similarity within each ROI. The curve (right) indicates the resampled (N = 5000) distribution of the effect size between ROIs, as measured by Cohen’s d. Cohen’s d of PC-IP minus MC: −2.1 (95% CI: [−2.22, −1.99]).

We also compared the PC-IP BCI RDM with able-bodied fMRI MC RDMs, which have been previously shown to match the patterns of natural hand use (Ejaz et al., 2015). The able-bodied MC and PC-IP fMRI finger organizations are similar in that they represent the thumb distinctly from the other fingers, but PC-IP fMRI signals represent each of the non-thumb fingers similarly while MC distinguishes between all five fingers (Figure 2—figure supplement 1). Interestingly, PC-IP BCI finger representations matched the able-bodied fMRI finger representational structure in the MC (Figure 2—figure supplement 1) even better than that of able-bodied PC-IP (Figure 2—figure supplement 5). The WUC similarity with the MC RDM was close to the noise ceiling (Methods), indicating that the MC RDM matches participant NS’s data better than almost any other model could (see Discussion).

Figure 2—figure supplement 5.

Finger representational structure of the tetraplegic individual, measured at the junction of the postcentral and intraparietal sulci (PC-IP), matches fMRI representational dissimilarity matrices (RDMs) from motor cortex (MC) even better than fMRI RDMs from PC-IP.

(a) Measured RDMs match the able-bodied MC fMRI RDM better than they match the able-bodied PC-IP fMRI RDM (p = 1.9 × 10–6; two-tailed t-test, 1000 bootstrap samples over 10 sessions), as measured by the whitened unbiased cosine similarity (Diedrichsen et al., 2021) (WUC) (Methods). Violin plot: solid horizontal lines indicate the median WUC over bootstrap samples, and dotted lines indicate the first and third quartiles. Noise ceiling: gray region estimates the best possible model fit (Methods). Asterisks denote significant differences at ***p < 0.001. Similar to Figure 2f. (b) Paired Gardner–Altman estimation plot (Ho et al., 2019) of the similarity (WUC) between participant NS (average RDM across sessions) and individual MC and PC-IP RDMs from able-bodied fMRI participants. The slopegraph’s connected points (left) show each fMRI participant’s (N = 36) MC and PC-IP similarities with participant NS’s mean finger RDM. Mean difference between PC-IP and MC similarities (right) presented as Cohen’s d (N = 5000 bootstrap samples).

Representational structure did not trend toward task optimum

Next, we investigated whether the BCI finger representational structure matched that of able-bodied individuals consistently or whether the representational structure changed over time to improve BCI performance. The task-optimal structure hypothesis (3) predicted that the BCI RDMs would trend to optimize for the task statistics (unstructured model, Figure 2c) as the participant gained experience with the BCI task. However, we did not find conclusive evidence for a trend from the able-bodied model toward the unstructured model (linear-model session × model interaction: t(6) = 0.50, one-tailed t-test p = 0.32, Bayes factor [BF] = 0.66) (Figure 3a). Indeed, participant NS’s finger RDMs were largely consistent across different recording sessions (average pairwise correlation, excluding the diagonal: r = 0.90 ± SD 0.04, min 0.83, max 0.99).

Figure 3.

Hand representation changed minimally after weeks of brain–computer interface (BCI) control.

(a) Slope comparison shows that the model fit did not trend toward the unstructured model over sessions (p = 0.32). (b) The distance between high-error finger pairs (middle-ring and ring-pinky) did not increase across sessions or runs (within sessions), as shown by partial regression plots. Distance metric: cross-validated Mahalanobis, averaged across runs (for the session plot) or averaged across sessions (for the run plot). The black line indicates linear regression. The gray shaded region indicates a 95% confidence interval. Each run consisted of 8 presses per finger. (c) Minimal change in representational structure between early and late sessions or between early and late runs. Mean representational dissimilarity matrix (RDM), when grouped by sessions (top row) or individual runs (bottom row). Grouped into early half (left column) or late half (center column). Multidimensional scaling (MDS) visualization (right column) of early (opaque) and late (translucent) representational structures after Generalized Procrustes alignment (without scaling, to allow distance comparisons).

Brain–computer interface (BCI) classification errors could have encouraged inter-finger distances to increase to improve separability, but this did not occur. Inter-finger distances instead decreased slightly (across sessions: t(8) = −4.0, two tailed t-test p = 0.004; across runs within sessions: t(82) = −2.4, two-tailed t-test p = 0.019), although the effect size was very small (across sessions: Cohen’s = 0.008; across runs within sessions: =0.005). Markers indicate average pairwise distance for each finger pair and session (top) or run-within-session (bottom).

We considered whether learning, across sessions or within sessions, could have caused smaller-scale changes in the representational structure. The observed representational structure, where middle-ring and ring-pinky pairs had relatively small distances, was detrimental to classification performance. The majority (70%) of the online classification errors were middle-ring or ring-pinky confusions (Figure 1b). Due to these systematic errors, one might reasonably predict that plasticity mechanisms would improve control by increasing the inter-finger distances between the confused finger pairs. Contrary to this prediction, the middle-ring and ring-pinky distances did not increase over the course of the experiment (across sessions: t(8) = −4.5, one-tailed t-test p > 0.99, BF = 0.03; across runs within sessions: t(82) = −0.45, one-tailed t-test p = 0.67, BF = 0.12) (Figure 3b). When analyzing all finger pairs together, the inter-finger distances also did not increase (across sessions: t(8) = −4.0, one-tailed t-test p = 0.98, BF = 0.01; across runs within sessions: t(74) = −2.4, one-tailed t-test p = 0.99, BF = 0.02), as visualized by the similarity between the average early-half RDM and the average late-half RDM (Figure 3c). These analyses demonstrate that the representational structure did not trend toward the task optimum (Figure 2c) with experience, ruling out the task-statistics hypothesis (3).

Finger representational structure is motor-like and then somatotopic

PPC is hypothesized to overcome inherent sensory delays by computing an internal forward model for rapid sensorimotor control (Mulliken et al., 2008; Wolpert et al., 1998; Desmurget and Grafton, 2000.) The forward model integrates an efference copy of motor signals and delayed sensory feedback to dynamically predict the state of the body. The hypothesized forward-model role would predict that the representational structure changes over the time course of each movement, with an early motor-command-like component during movement initiation. To investigate this temporal evolution, we modeled the representational structure of finger movements at each timepoint as a non-negative linear combination (Kietzmann et al., 2019) of potentially predictive models (Figure 4a).

We considered three models (Ejaz et al., 2015) that could account for representational structure: hand usage, muscle activation, and somatotopy. The hand-usage model (Figure 4b) predicts that the neural representational structure should follow the correlation pattern of finger kinematics during natural hand use. The muscle activation model (Figure 4c) predicts that the representational structure should follow the coactivation patterns of muscle activity during individual finger movements. The somatotopy model (Figure 4d) predicts that the representational structure should follow the spatial layout of the body, with neighboring fingers represented similarly to each other (Ejaz et al., 2015; Schellekens et al., 2018).Schellekens et al., 2018 While somatotopy usually refers to physical spaces that resemble the body, here we use the term broadly to describe encoding spaces that resemble the body.

Because the hand-usage model is nearly multicollinear with the muscle and somatotopy models (variance inflation factor: VIFusage,OLS = VIFusage,NNLS = 20.9, Methods), we first reduced the number of component models. Through a model selection procedure (Methods), we found that the hand usage + somatotopy and muscle + somatotopy model combinations matched the data best (Figure 4—figure supplement 1), with the muscle + somatotopy model matching the data marginally better. Thus, in the main text, we present our temporal analysis using the muscle and somatotopy component models.

Figure 4—figure supplement 1.

Fit between measured representational dissimilarity matrix (RDM) and linear combinations of models.

Violin plot of WUC similarity between the measured RDM (N = 1000 bootstrap samples over 10 sessions) and the corresponding model combination. Violin plot: solid horizontal lines indicate the mean WUC over bootstrap samples, and dotted lines indicate the first and third quartiles. Horizontal lines (above) indicate significance groups, where the circle-indicated model is significant over the vertical-tick-indicated models (two-tailed t-test, q < 0.01, FDR-corrected for 28 model-pair comparisons). For example, the muscle + somatotopy combined model is significant over the individual muscle, hand usage, somatotopy, combined muscle + hand usage, and pairwise-equidistant/unstructured (null) models.

Figure 4e shows the decomposition of the representational structure into the muscle and somatotopy component models. The results show a dynamic structure, with the muscle model emerging 170 ms earlier than the somatotopy model (p = 0.002, two-sided Wilcoxon signed-rank test). This timing difference was consistent across individual sessions (Figure 4—figure supplement 2) and task contexts, such as the calibration task (Figure 4—figure supplement 3). Indeed, the transition from the muscle model (Figure 4c) to the somatotopy model (Figure 4d) is visually apparent when comparing the average RDMs at the early (muscle-model-like) late (somatotopic) phases of movement (Figure 4e).

Figure 4—figure supplement 2.

Temporal delays between component models are consistent across single sessions.

When linear modeling within single sessions, the muscle model (blue) consistently preceded the somatotopy model (orange). Time difference: 170 ± 66 ms (SD across sessions) (p = 0.002, two-sided Wilcoxon signed-rank test). Line styles indicate session. Related to Figure 4e.

Figure 4—figure supplement 3.

Representational dynamics are robust across tasks and model combination choices.

(a) Representational dynamics analysis shows an early fit to the hand-usage model and a late fit to the somatotopy model. Confidence intervals indicate ± standard error of the mean (SEM) across sessions. Related to Figure 4e. (b) Representational dynamics analysis shows a consistent delay between models during the calibration task. Note: The absolute timing differs from the main task because the calibration task does not require an initial saccade to read the cue. Related to Figure 4e.

These temporal dynamics were robust to our model selection procedure, demonstrating a similar timing difference for the hand usage + somatotopy model combination (Figure 4—figure supplement 3).

Discussion

Neural prosthetic control of individual fingers using recordings from PC-IP

We found that participant NS could robustly control the movement of individual fingers using a neural prosthetic in a variety of contexts (Figure 1, Figure 1—figure supplement 4, Video 1), even after years of paralysis. Her BCI control accuracy exceeded the previous best of other five-finger, online BCI control studies (Hotson et al., 2016; Jorge et al., 2020). These results establish PC-IP as a candidate implant region for dexterous neural prostheses.

Connecting BCI studies to basic neuroscience

Although previous studies have shown that the anterior intraparietal (AIP) area of PPC is involved in whole-hand grasping (Schaffelhofer and Scherberger, 2016; Klaes et al., 2015; Murata et al., 2000), our work is the first to discover individual finger representations in PPC (Figure 1—figure supplement 5). Likewise, many other BCI studies with tetraplegic participants have contributed to basic neuroscience, deepening our understanding of the human cortex (Stavisky et al., 2019; Willett et al., 2020; Rutishauser et al., 2018; Zhang et al., 2017; Aflalo et al., 2020; Chivukula et al., 2021). A frequent (Flesher et al., 2016; Armenta Salas et al., 2018; Stavisky et al., 2019; Willett et al., 2020; Fifer et al., 2022; Chivukula et al., 2021; Andersen and Aflalo, 2022) discussion question is: how well do these findings generalize to the brains of able-bodied individuals? Specifically, do the observed phenomena result from partial reorganization (Kambi et al., 2014; Nardone et al., 2013) after spinal cord injury, or do they reflect intact motor circuits, preserved from before injury (Makin and Bensmaia, 2017)?

Early human BCI studies (Hochberg et al., 2006; Collinger et al., 2013) recorded from the MC and found that single-neuron directional tuning is qualitatively similar to that of able-bodied non-human primates (NHPs) (Hochberg et al., 2006; Georgopoulos et al., 1982). Many subsequent human BCI studies have also successfully replicated results from other classical NHP neurophysiology studies (Hochberg et al., 2012; Collinger et al., 2013; Gilja et al., 2015; Bouton et al., 2016; Ajiboye et al., 2017; Brandman et al., 2018; Aflalo et al., 2015), leading to the general heuristic that the sensorimotor cortex retains its major properties after spinal cord injury (Andersen and Aflalo, 2022). This heuristic further suggests that BCI studies of tetraplegic individuals should generalize to able-bodied individuals. However, this generalization hypothesis had lacked direct, quantitative comparisons between tetraplegic and able-bodied individuals. Thus, as human BCI studies expanded beyond replicating results and beganto challenge conventional wisdom, neuroscientists questioned whether cortical reorganization could influence these novel phenomena (see Discussions of Flesher et al., 2016; Armenta Salas et al., 2018; Stavisky et al., 2019; Willett et al., 2020; Fifer et al., 2022; Chivukula et al., 2021; Andersen and Aflalo, 2022). As an example of a novel discovery, a recent BCI study found that the hand knob of tetraplegic individuals is directionally tuned to movements of the entire body (Willett et al., 2020), challenging the traditional notion that primary somatosensory and motor subregions respond selectively to individual body parts (Penfield and Boldrey, 1937). Given the brain’s capacity for reorganization (Jain et al., 2008; Kambi et al., 2014), could these BCI results be specific to cortical remapping? Detailed comparisons with able-bodied individuals, as shown here, help shed light on this question.

Matching finger organization between tetraplegic and able-bodied participants

We asked whether participant NS’s BCI finger representations resembled that of able-bodied individuals or whether her finger representations had reorganized after paralysis. Single-neuron recordings of PC-IP during individuated finger movements for able bodied humans are not available for comparison. However, many fMRI studies have characterized finger representations (Kikkert et al., 2021; Ejaz et al., 2015; Kikkert et al., 2016; Yousry et al., 1997), and RSA has previously shown RDM correspondence between fMRI and single-neuron recordings of another cortical region (inferior temporal cortex) (Kriegeskorte et al., 2008b). This match was surprising because single-neuron and fMRI recordings differ fundamentally; single-neuron recordings sparsely sample 102 neurons in a small region, while fMRI samples 104–106 neurons/voxel (Kriegeskorte and Diedrichsen, 2016; Guest and Love, 2017). The correspondence suggested that RSA might identify modality-invariant neural organizations (Kriegeskorte et al., 2008b), so here we used fMRI recordings of human PC-IP as an able-bodied model.

We found that participant NS exhibited a consistent finger representational structure across sessions, and this representational structure matched the able-bodied fMRI model better than the task-optimal, unstructured model (Figure 2). When compared with individual able-bodied participants, participant NS’s finger organization was also quite typical, in part due to the relative variability in PC-IP fMRI representational structure across able-bodied participants.

The MC fMRI finger representation is well studied and has been shown to reflect the patterns of natural hand use (Ejaz et al., 2015; Kikkert et al., 2016; Wesselink et al., 2019), so we also considered a model constructed from MC fMRI recordings. Compared to the PC-IP fMRI finger representation, MC represents the non-thumb fingers more distinctly from each other (Figure 2—figure supplement 1). Interestingly, participant NS’s finger RDMs more strongly matched the able-bodied MC fMRI model, reaching similarities close to the theoretical maximum (Figure 2—figure supplement 3 and Figure 2—figure supplement 5). This result does obscure a straightforward interpretation of the RSA results – why does our recording area match MC better than the corresponding implant location? Several factors might contribute, including differing neurovascular sensitivity to the early and late phases of the neural response (Figure 4e), heterogeneous neural organizations across the single-neuron and voxel spatial scales (Kriegeskorte and Diedrichsen, 2016; Guest and Love, 2017; Arbuckle et al., 2020), or mismatches in functional anatomy between participant NS and standard atlases (Eickhoff et al., 2018). Furthermore, fMRI BOLD contrast is thought to reflect cortical inputs and intracortical processing (Logothetis et al., 2001). Thus, the match between PC-IP spiking output and MC fMRI signals could also suggest that PC-IP sends signals to MC, thereby driving the observed MC fMRI structure.

Even so, it is striking that participant NS’s finger representation matches the neural and hand use patterns (Figure 4b and Figure 4—figure supplement 1) of able-bodied individuals. Despite the lack of overt movement or biomechanical constraints (Lang and Schieber, 2004), the measured finger representation still reflected these usage-related patterns. This result matches recent sensorimotor cortex studies of tetraplegic individuals, where MC decoding errors (Jorge et al., 2020) and S1 finger somatotopy (Kikkert et al., 2021) appeared to reflect able-bodied usage patterns. Taken together with our dynamics analyses (see Discussion), the evidence supports the interpretation that motor representations are preserved after paralysis. Comparisons with single-neuron recordings from able-bodied participants would validate this interpretation. although such recordings may be difficult to acquire.

Able-bodied-like finger representation is not explained by learning

Hand use patterns shape neural finger organization (Ejaz et al., 2015; Kikkert et al., 2016; Wesselink et al., 2019), so we also considered the possibility that participant NS’s able-bodied-like representational structure emerged from BCI usage patterns after paralysis. Contrary to this hypothesis, her BCI finger representational structure changed minimally over weeks (Figure 3). Furthermore, even though participant NS’s representational structure contributed to BCI errors (Figure 1b) and she was anecdotally cognizant of which fingers would get confused, she did not increase the neural distance between fingers with experience. This relative stability suggests that the measured representational structure has been stable after paralysis, rather than emergent from BCI learning.

The stability of finger representations here suggests that BCIs can benefit from the pre-existing, natural repertoire (Hwang et al., 2013), although learning can play an important role under different experimental constraints. In our study, the participant received only a delayed, discrete feedback signal after classification (Figure 1a). Because we were interested in understanding participant NS’s natural finger representation, we did not artificially perturb the BCI mapping. When given continuous feedback, however, participants in previous BCI studies could learn to adapt to within-manifold perturbations to the BCI mapping (Sadtler et al., 2014; Vyas et al., 2018; Ganguly and Carmena, 2009; Sakellaridi et al., 2019). BCI users could even slowly learn to generate off-manifold neural activity patterns when the BCI decoder perturbations were incremental (Oby et al., 2019). Notably, learning was inconsistent when perturbations were sudden, indicating that learning is sensitive to specific training procedures.

To further understand how much finger representations can be actively modified, future studies could benefit from perturbations (Kieliba et al., 2021; Oby et al., 2019), continuous neurofeedback (Vyas et al., 2018; Ganguly and Carmena, 2009; Oby et al., 2019), and additional participants. Additionally, given that finger representations were dynamic (Figure 4e), learning could occur separately in the early and late dynamic phases. Time-variant BCI decoding algorithms, such as recurrent neural networks (Willett et al., 2021; Sussillo et al., 2012), could also help facilitate learning specific to different time windows of finger movement.

Representational dynamics are consistent with PPC as a forward model

In able-bodied individuals, PPC is thought to maintain a forward estimate of movement state (Mulliken et al., 2008; Wolpert et al., 1998; Desmurget and Grafton, 2000; Aflalo et al., 2015; McNamee and Wolpert, 2019). As such, PPC receives delayed multimodal sensory feedback and is hypothesized to receive efference copies of motor command signals (Mulliken et al., 2008; Andersen et al., 1997). This hypothesized role predicts that PPC houses multiple functional representations, each engaged at different timepoints of motor production.

To dissociate these neural processes, we performed a time-resolved version of RSA (Figure 4). We considered three component models: muscle, usage, and somatotopy (Ejaz et al., 2015). Our temporal analysis showed a consistent ordering: early emergence of the muscle model followed by the somatotopy model.

This ordering was consistent when exchanging the muscle and hand-usage component models (Figure 4 and Figure 4—figure supplement 3), as hand-usage and muscle activation patterns are strongly correlated for individual finger movements (Overduin et al., 2012). Therefore, we group these two models under the single concept of motor production. In the future, more complex multi-finger movements (Ejaz et al., 2015) would help distinguish between muscle and hand-usage models.

The somatotopy model predicts that neighboring fingers will have similar cortical activity patterns (Ejaz et al., 2015), analogous to overlapping Gaussian receptive fields (Schellekens et al., 2018). Gaussian receptive fields have been useful tools for understanding finger topographies within the sensorimotor cortex (Schellekens et al., 2018; Schellekens et al., 2021). In another study with participant NS, we found that the same PC-IP population encodes touch (Chivukula et al., 2021) with Gaussian-like receptive fields. Based on these results, the somatotopy model can be thought of as a sensory-consequence model. However, because participant NS has no sensation below her shoulders, we interpret the somatotopy model as the preserved prediction of the sensory consequences of a finger movement. These sensory outcome signals could be the consequence of internal computations within the PPC or could come from other structures important for body-state estimation, such as the cerebellum (McNamee and Wolpert, 2019).

The 170 ms timing difference we found roughly matches the delay between feedforward muscle activation and somatosensory afferents (Scott, 2016; Sollmann et al., 2017) in able-bodied individuals. Given PPC’s hypothesized role as a forward model, PPC likely integrates motor planning and production signals to predict sensory outcomes at such a timing (Mulliken et al., 2008; Wolpert et al., 1998; Desmurget and Grafton, 2000; McNamee and Wolpert, 2019). Alternatively, because participant NS cannot move overtly, the sensory-consequence model could also reflect the error between the internal model’s expected sensory outcomes and the actual (lack of) sensory feedback (Adams et al., 2013). In either scenario, the match in timing between BCI control and able-bodied individuals provides further evidence that the recorded motor circuits have preserved their functional role.

Stability of sensorimotor representations after paralysis

A persistent question in neuroscience has been how experience shapes the brain, and to what extent existing neural circuits can be modified. Early studies by Merzenich et al. showed that the primary somatosensory cortex reorganized after amputation, with intact body parts invading the deprived cortex (Merzenich et al., 1984; Qi et al., 2000; Pons et al., 1991). However, the authors also recognized that the amputated body part might persist in latent somatosensory maps. Since then, preserved, latent somatosensory representations have been demonstrated in studies of amputation (Makin and Bensmaia, 2017; Kikkert et al., 2016; Wesselink et al., 2019; Bruurmijn et al., 2017) and even paralysis (Kikkert et al., 2021; Flesher et al., 2016; Armenta Salas et al., 2018; Fifer et al., 2022). Overall, deafferentation appears to expand the remaining regions slightly, even while the pre-injury structure persists in the deafferented cortex (Makin and Bensmaia, 2017). Fewer studies have investigated sensorimotor plasticity beyond the primary somatosensory cortex and MC, but our results in PC-IP indicate that association areas can also remain stable after paralysis.

The topic of cortical reorganization has long been significant to the development of BCIs, particularly when deciding where to implant recording electrodes. If, as previously thought, sensory deprivation drives cortical reorganization and any group of neurons can learn to control a prosthetic (Fetz, 1969; Moritz and Fetz, 2011), the specific implant location would not affect BCI performance. However, our results and others (Smirnakis et al., 2005; Makin and Bensmaia, 2017; Kikkert et al., 2021; Hwang et al., 2013; Kikkert et al., 2016; Wesselink et al., 2019; Bruurmijn et al., 2017) suggest that the pre-injury properties of brain regions do affect BCI performance. Even though experience shapes neural organization (Merzenich et al., 1984; Ejaz et al., 2015; Wesselink et al., 2019), representations may be remarkably persistent once formed (Kikkert et al., 2021; Wesselink et al., 2019). Thus, even though BCIs bypass limbs and their biomechanical constraints (Lang and Schieber, 2004), BCIs may still benefit from tapping into the preserved, natural (Hwang et al., 2013) movement repertoire of motor areas.

As BCIs enable more complex motor skills, such as handwriting (Willett et al., 2021), future studies could investigate whether these complex skills also retain their pre-injury representational structure. For example, does a tetraplegic participant’s BCI handwriting look like their physical, pre-injury handwriting? These results will have important implications for the design of future neural prosthetics.

Materials and methods

Data collection

Study participant

The study participant NS has an AIS-A spinal cord injury at cervical level C3–C4 that she sustained approximately 10 years before this study. Participant NS cannot move or feel her hands. As part of a BCI clinical study (ClinicalTrials.gov identifier: NCT01958086), participant NS was implanted with two 96-channel Neuroport Utah electrode arrays (Blackrock Microsystems model numbers 4382 and 4383). She consented to the surgical procedure as well as to the subsequent clinical studies after understanding their nature, objectives, and potential risks. All procedures were approved by the California Institute of Technology, Casa Colina Hospital and Centers for Healthcare, and the University of California, Los Angeles Institutional Review Boards.

Multielectrode array implant location

The recording array was implanted over the hand/limb region of the left PPC at the junction of the intraparietal sulcus with the postcentral sulcus (Figure 1—figure supplement 1; Talairach coordinates [−36 lateral, 48 posterior, 53 superior]). We previously (Klaes et al., 2015; Zhang et al., 2017; Aflalo et al., 2015) referred to this brain area as the AIP area, a region functionally defined in NHPs. Here, we describe the implanted area anatomically, denoting it the PC-IP area. More details regarding the methodology for functional localization and implantation can be found in Aflalo et al., 2015.

Neural data preprocessing

Using the NeuroPort system (Blackrock Microsystems), neural signals were recorded from the electrode array, amplified, analog bandpass-filtered (0.3 Hz to 7.5 kHz), and digitized (30 kHz, 250 nV resolution). A digital high-pass filter (250 Hz) was then applied to each electrode.

Threshold crossings were detected at a threshold of −3.5× RMS (root-mean-square of an electrode’s voltage time series). Threshold crossings were used as features for in-session BCI control. For all other analyses, we used k-medoids clustering on each electrode to spike-sort the threshold crossing waveforms. The first principal components were used as input features to k-medoids, where was selected for each electrode to account for 95% of waveform variance. The gap criteria (Tibshirani et al., 2001) were used to determine the number of waveform clusters for each electrode.

Experimental setup

Recording sessions

Experiments were conducted in 2–3 hr recording sessions at Casa Colina Hospital and Centers for Healthcare. All tasks were performed with participant NS seated in her motorized wheelchair with her hands resting prone on the armrests. Participant NS viewed text cues on a 27-inch LCD monitor that occupied approximately 40 degrees of visual angle. Cues were presented using the psychophysics toolbox (Brainard, 1997) for MATLAB (Mathworks).

The data were collected on 9 days over 6 weeks. Almost all experiment days were treated as individual sessions (i.e., the day’s recordings were spike-sorted together). The second experiment day (2018-09-17) was an exception, with data being recorded in a morning period and an afternoon period with a sizable rest in between. To reduce the effects of recording drift, we treated the two periods as separate sessions (i.e., spike-sorted each separately) for a total of 10 sessions. Each session can thus be considered a different resampling of a larger underlying neural population, with both unique and shared neurons each session. We did not rerun the calibration task for the afternoon session of the second experiment day (2018-09-17), resulting in nine sessions of the calibration task for Figure 1—figure supplement 4b.

Each session consisted of a series of 2–3 min, uninterrupted ‘runs’ of the task. The participant rested for a few minutes between runs as needed.

Calibration task

At the beginning of each recording day, the participant performed a reaction-time finger flexion task (Figure 1—figure supplement 2; denoted ‘calibration task’ in the Results) to train a finger classifier for subsequent runs of the primary task. On each trial, a letter appeared on the screen (e.g., ‘T’ for thumb). The participant was instructed to immediately flex the corresponding finger on the right hand (contralateral to the implant), as though pressing a key on a keyboard. The condition order was block-randomized, such that each condition appeared once before repetition.

Finger flexion grid task

In the primary task, movement cues were arranged in a 3 × 4 grid of letters on the screen (Figure 1a). Each screen consisted of two repetitions each of T (thumb), I (index), M (middle), R (ring), P (pinky/little), and X (No-Go) arranged randomly on the grid. Each trial lasted 3 s. At the beginning of each trial, a new cue was randomly selected with a crosshairs indicator, which jittered randomly to prevent letter occlusion. Each cue was selected once (for a total of 12 trials) before the screen was updated to a new arrangement. Each run consisted of three to four screens.

On each trial, the participant was instructed to immediately (1) saccade to the cued target, (2) fixate, and (3) attempt to press the corresponding finger. During both movement and No-Go trials, the participant was instructed to fixate on the target at least until the visual classification feedback was shown. The cue location randomization was used to investigate whether cue location would affect movement representations.

On each trial, 1.5 s after cue presentation, the classifier decoded the finger movement and presented its prediction via text feedback overlaid on the cue.

No-movement control task

The control task was similar to the primary task, except that the subject was instructed to saccade to each cued letter and fixate without attempting any finger movements. No classification feedback was shown.

Statistical analysis

Unit selection

Single-unit neurons were identified using the k-medoids clustering method, as described in the Neural Data Preprocessing section. Analyses in the main text used all identified units, regardless of sort quality. With spike-sorting, there is always the possibility that a single waveform cluster corresponds to activity from multiple neurons. To confirm that potential multiunit clustering did not bias our results, we repeated our analyses using only well-isolated units (Figure 4—figure supplement 4).

Figure 4—figure supplement 4.

Well-isolated single neurons of the tetraplegic participant match the finger representational structure of able-bodied individuals.

(a) Histogram of L-ratio, a spike-sorting cluster metric. Threshold for well-isolated units: 33% quantile (Lratio < 10−1.1). (b) Representational dissimilarity matrices (RDMs) calculated only using well-isolated units, using the cross-validated Mahalanobis distance. Similar to Figure 2d and Figure 2—figure supplement 2a. (c) Whitened unbiased similarity (WUC) between measured (b) RDMs (using only well-isolated units) and model predictions (Figure 2b, c), showing that the measured RDMs match the able-bodied fMRI RDM significantly better than they match the unstructured model (p = 3.1 × 10–10, two-tailed t-test) and the SPLa fMRI RDM (p = 1.7 × 10–8). Error bars indicate ± standard error of the mean (SEM). Noise ceiling: gray region estimates the best possible model fit (Methods). Gray downward-semicircle indicates that the noise ceiling is significantly higher (p < 0.001) than the fit of the SPLa fMRI RDM and the unstructured model. Asterisks denote significant differences at ***p < 0.001. Similar to Figure 2f. (d) Representational dynamics analysis, repeated using only well-isolated units, shows an early fit to the muscle model and a late fit to the somatotopy model. Confidence intervals indicate ± SEM across sessions. Similar to Figure 4e.

Well-isolated single units were identified using the L-ratio metric (Schmitzer-Torbert et al., 2005). The neurons corresponding to the lowest third of L-ratio values (across all sessions) were selected as ‘well-isolated’. This corresponded to a threshold of dividing well-isolated single units and potential multiunits (Figure 4—figure supplement 4).

Single-unit tuning to finger flexion

We calculated the firing rate for each neuron in the window [0.5, 1.5] s after cue presentation. To calculate significance for each neuron (Figure 1—figure supplement 5), we used a two-tailed t-test comparing each movement’s firing rate to the No-Go firing rate. A neuron was considered significantly tuned to a movement if p < 0.05 (after FDR correction). We also computed the mean firing rate change between each movement and the No-Go condition. If a neuron was significantly tuned to at least one finger, we denoted the neuron’s ‘best finger’ as the significant finger with the largest effect size (mean firing rate change). For each neuron and finger, we also calculated the discriminability index (d′, RMS SD) between the baseline (No-Go) firing rate and the firing rate during finger movement.

In Figure 1—figure supplement 5, neurons were pooled across all 10 sessions. Neurons with mean firing rates less than 0.1 Hz were excluded to minimize sensitivity to discrete spike counting.

Finger classification

To classify finger movements from firing rate vectors, we used linear discriminant analysis (LDA) with diagonal covariance matrices (Dudoit et al., 2002) (a form of regularization); diagonal LDA is also equivalent to Gaussian Naive Bayes (GNB) when GNB assumes that all classes share a covariance matrix.

We used data from the calibration task to fit the BCI classifier. Input features (firing rate vectors) were calculated by counting the number of threshold crossings on each electrode during a 1-s time window within each trial’s movement execution phase. The exact time window was a hyperparameter for each session and was chosen to maximize the cross-validated accuracy on the calibration dataset. To prevent low-firing rate discretization effects, we excluded electrodes with mean firing rates less than 1 Hz. This classifier was then used in subsequent online BCI control for the main task (finger flexion grid).

During online control of the finger flexion grid task, input features were similarly constructed by counting the threshold crossings from each electrode in a 1-s time window. This time window was fixed to [0.5, 1.5] s after cue presentation. The window start time was chosen based on the estimated saccade latency in the first experimental session. The saccade latency was estimated by taking the median latency for the subject to look >80% of the distance between targets. The analysis window was a priori determined to be 1 s; this choice was supported post hoc by a sliding window analysis (not shown), which confirmed that finger movements could be accurately classified up to 1.6 s after cue. The online classifier was occasionally retrained using data from this main task, usually every four run blocks.

Offline classification accuracy (Figure 1—figure supplement 4) was computed using leave-one-out cross-validation (within each session). We used features from the same time window as the online-control task. However, offline analyses used firing rates after spike-sorting, instead of raw threshold crossings.

In the Results, reported classification accuracies aggregate trials over all sessions (as opposed to averaging the accuracies across sessions with different numbers of trials). Reported SDs indicate variability across sessions, weighted by the number of trials in each session. To visualize confusion matrices, trials were pooled across sessions. Confusion matrix counts were normalized by row sum (true label) to display confusion percentages.

In the first session (2018-09-10), the No-Go condition (X) was not cued in the calibration task, so the classifier did not output No-Go predictions during that session. However, No-Go was cued in the main task; these 84 No-Go trials were thus excluded from the online-control accuracy metrics (Figure 1b and Figure 1—figure supplement 3), but they were included in the offline cross-validated confusion matrix (Figure 1—figure supplement 4).

Cross-validated neural distance

We quantified the dissimilarity between the neural activity patterns of each finger pair , using the cross-validated (squared) Mahalanobis distance (Schütt et al., 2019; Nili et al., 2014):

where and denote independent partitions of the trials, are the partition-specific noise covariance matrices, are the trial measurements of firing rate vectors for conditions , and normalizes for the number of neurons. The units of are .

The cross-validated Mahalanobis distance, also referred to as the ‘crossnobis’ distance (Schütt et al., 2019), measures the separability of multivariate patterns, analogous to LDA classification accuracy (Nili et al., 2014). To generate independent partitions and for each session, we stratified-split the trials into five non-overlapping subsets. We then calculated the crossnobis distance for each possible combination of subsets and averaged the results. Cross-validation ensures that the (squared) distance estimate is unbiased; when the underlying distributions are identical (Walther et al., 2016). The noise covariance was regularized (Ledoit and Wolf, 2003) to guarantee invertibility.

Similar results were also obtained when estimating neural distances with the cross-validated Poisson symmetrized KL-divergence (Schütt et al., 2019; Figure 2—figure supplement 3).

Representational models

We used RDMs to describe both the type and format of information encoded in a recorded population. To make these RDMs, we calculated the distances between each pair of finger movements and organized the 10 unique inter-finger distances into a -sized RDM (Figure 2d). Conveniently, the RDM abstracts away the underlying feature types, enabling direct comparison with RDMs across brain regions (Kietzmann et al., 2019), subjects, or recording modalities (Kriegeskorte et al., 2008b).

We also used RDMs to quantify hypotheses about how the brain might represent different actions. In Figure 2b, we generated an able-bodied model RDM using fMRI data from two independent studies, Kieliba et al., 2021 (N = 29, pre-intervention, right hand, 3T scans) and Ejaz et al., 2015 (N = 7, no intervention, right hand, 7T scans). The fMRI ROI was selected to match participant NS’s anatomical implant location (PC-IP). Specifically, a 4-mm geodesic distance around vertex 7123 was initially drawn in fs_LR_32k space, then resampled onto fsaverage. The RDM for each subject was then calculated using the cross-validated (squared) Mahalanobis distance between fMRI activity patterns. Based on a permutation shuffle test, RDMs were similar between the studies’ groups of participants, so we aggregated the RDMs into a single dataset here. The MC RDMs (Figure 2—figure supplement 1) used data from the same scans (Ejaz et al., 2015; Kieliba et al., 2021), with ROIs covering Brodmann area 4 near the hand knob of the precentral gyrus.

In Figure 4 and its supplemental figures, we decomposed the data RDMs into model RDMs borrowed from Ejaz et al., 2015. The hand-usage model was constructed using the velocity time series of each finger’s MCP joint during everyday tasks (Ingram et al., 2008). The muscle activity model was constructed using EMG activity during single- and multi-finger tasks. The somatotopic model is based on a cortical sheet analogy and assumes that finger activation patterns are linearly spaced Gaussian kernels across the cortical sheet. Further modeling details are available in the methods section of Ejaz et al., 2015.

Comparing representational structures

We used the rsatoolbox Python library (Schütt et al., 2019) to calculate data RDMs and compare them with model RDMs (RSA) (Kriegeskorte et al., 2008a).

To quantify model fit, we used the whitened unbiased RDM cosine similarity (WUC) metric (Diedrichsen et al., 2021), which (Diedrichsen et al., 2021) recommend for models that predict continuous real values. We used WUC instead of Pearson correlation for two reasons (Diedrichsen et al., 2021). First, cosine similarity metrics like WUC properly exploit the informative zero point; because we used an unbiased distance estimate, indicates that the distributions are identical. Second, Pearson correlation assumes that observations are independent, but the elements of each RDM covary (Diedrichsen et al., 2021) because the underlying dataset is shared. For example, the (thumb, middle)-pairwise dissimilarity uses the same thumb data as the (thumb, ring)-pairwise dissimilarity.

Like correlation similarities, a larger WUC indicates a better match, and the maximum WUC value is 1. However, cosine similarities like WUC are often larger than the corresponding correlation values or are even close to 1 (Diedrichsen et al., 2021). Thus, while correlation values can be compared against a null hypothesis of 0-correlation, WUC values should be interpreted by comparing against a baseline. The baseline is usually (Diedrichsen et al., 2021) chosen to be a null model where all conditions are pairwise-equidistant (and would thus correspond to 0-correlation). In this study, this happens to correspond to the unstructured model. For more details about interpreting the WUC metric, see Diedrichsen et al., 2021.

To demonstrate that our model comparisons were robust to the specific choice of RDM similarity metric, we also show model fits using whitened Pearson correlation in Figure 2—figure supplement 3. Whitened Pearson correlation is a common alternative to WUC (Diedrichsen et al., 2021).

Noise ceiling for model fits

Measurement noise and behavioral variability cause data RDMs to vary across repetitions, so even a perfect model RDM would not achieve a WUC similarity of 1. To estimate the noise ceiling (Nili et al., 2014) (the maximum similarity possible given the observed variability between data RDMs), we assume that the unknown, perfect model resembles the average RDM. Specifically, we calculated the average similarity of each individual-session RDM (Figure 2—figure supplement 2) with the mean RDM across all other sessions (i.e., excluding that session):

where is the WUC similarity function, is the number of RDMs, refers to a single RDM from an individual session, and is the ‘lower’ noise ceiling. This noise ceiling is analogous to leave-one-out-cross-validation. If a model achieves the noise ceiling, the model fits the data well (Nili et al., 2014).

Measuring changes in the representational structure

To assess the effect of BCI task experience on the inter-finger distances, we performed a linear regression analysis (Figure 3b and Figure 3—figure supplement 1). We first subdivided each session’s dataset into individual runs and calculated separate RDMs for each (session, run) index. We then used linear regression to predict each RDM’s (squared) inter-finger distances from the session index, run index, and finger pair:

Figure 3—figure supplement 1.

Inter-finger distances did not increase across sessions or within sessions.

Brain–computer interface (BCI) classification errors could have encouraged inter-finger distances to increase to improve separability, but this did not occur. Inter-finger distances instead decreased slightly (across sessions: t(8) = −4.0, two tailed t-test p = 0.004; across runs within sessions: t(82) = −2.4, two-tailed t-test p = 0.019), although the effect size was very small (across sessions: Cohen’s = 0.008; across runs within sessions: =0.005). Markers indicate average pairwise distance for each finger pair and session (top) or run-within-session (bottom).

where is the average inter-finger distance, is the coefficient for finger-pair , is the session index, and is the run index. would suggest that RDMs are dependent on experience across sessions. would suggest that RDMs depend on experience across runs within a session. For t-tests, we conservatively estimated the degrees-of-freedom as the number of RDMs, because the individual elements of each RDM covary and thus are not independent (Diedrichsen et al., 2021). The effect sizes for the session-index predictor and the run-index predictor were quantified using Cohen’s (Cohen, 1988), comparing against the finger-pair-only model as a baseline.

For t-tests without significant differences, we also calculated BFs to determine the likelihood of the null hypothesis, using the common threshold that BF <1/3 substantially supports the null hypothesis (Dienes, 2014). BFs were computed using the R package BayesFactor (Morey et al., 2015) with default priors. To calculate BFs for one-sided t-tests (for example, ), we sampled (N = 106) from the posterior of the corresponding two-sided t-test (), calculated the proportion of samples that satisfied the one-sided inequality, and divided by the prior odds (Morey and Wagenmakers, 2014) ().

Linear combinations of models

We modeled the finger RDM (in vector form) as a zero-intercept, non-negative linear combination (Jozwik et al., 2016) of potentially predictive model RDMs: usage, muscle, and somatomorphic (Figure 4).

First, we used the VIF to assess multicollinearity between the hypothesized models. For each model (e.g., usage), we calculated the standard, ordinary least squares (OLS)-based VIF (VIFusage,OLS), and we also calculated a modified VIF (VIFusage,NNLS) based on non-negative least squares (NNLS).

where is the from an OLS regression predicting RDM from all other RDMs . VIFOLS values above a threshold indicate that multicollinearity is a problem; VIF >5 or VIF >10 are common thresholds (James et al., 2013). Here, we constrained the linear combination coefficients to be non-negative, which can sometimes mitigate multicollinearity. Thus, we also calculated VIFNNLS, which follows the same equation above, except that we use NNLS to predict from .

Because multicollinearity was a problem here, we next determined the best subset of model RDMs to use. We used NNLS to predict the data RDM from the model RDMs. We estimated the model fits using leave-one-session-out cross-validation. To estimate model-fit uncertainty, we bootstrapped RDMs (sessions) over the cross-validation procedure (Schütt et al., 2019). We then used the ‘one-standard error’ rule (James et al., 2013) to select the best parsimonious model, choosing the simplest model within one standard error of the best model fit.

Representational dynamics analysis

To investigate how the finger movement representational structure unfolds over time, we used a time-resolved version of RSA (Kietzmann et al., 2019; Figure 4a). At each timepoint within a trial, we computed the instantaneous firing rates by binning the spikes in a 200-ms time window centered at that point. These firing rates were used to calculate cross-validated Mahalanobis distances between each pair of fingers and generate an RDM. Snapshots (Figure 4e) show single-timepoint RDMs averaged across sessions.

The temporal sequence of RDMs constitutes an RDM movie (size ) that visualizes the representational trajectory across the trial duration. RDM movies were computed separately for each recording session. At each timepoint, we linearly decomposed the data RDM into the component models using non-negative least squares. Because the component models were multicollinear, component models were limited to the subsets chosen in the previous model reduction step. Each component RDM was normalized by its vector length (ℓ2-norm) before decomposition to allow comparison between coefficient magnitudes. We used bootstrapped sampling of RDMs across sessions and decomposed the bootstrap-mean RDM to generate a confidence intervals on the coefficients.

We computed the start time of each model component as the time at which the corresponding mixture coefficient exceeded 0.2 (about 25% of the median peak coefficient across models and sessions).

Using data from an tetraplegic individual, the authors show that the neural representations for attempted single finger movements after multiple years after the injury is still organized in a way that is typical for healthy participants. They also show that the representational structure does not change during task training on a simple finger classification task – and that the representational structure, even without active motor outflow or sensory inflow, switches from a motor representation to a sensory representation during the trial. The analyses are convincing, and the results have important implications for the use and training of BCI devices in humans.

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Preserved motor representations after paralysis" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, including Jörn Diedrichsen as Reviewing Editor and Reviewer #1, and the evaluation has been overseen by Richard Ivry as the Senior Editor.

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

1) One important question to address clearly in the revision is the mismatch of the representational structure of the recorded PPC site with the SPLa ROI from the imaging studies (and the match with M1). Reviewer #1 raises the important questions that fMRI ROI is likely not equivalent to the recorded site and should be modified and the finding reinvestigated. If the discrepancy persists after the reanalysis, Reviewer #2, and #3 make suggestions for a clearer discussion of the findings. Overall, the claim of the stability of pre-injury organization probably needs to be tempered a bit.

2) The stability of representations across learning is an interesting finding, but Reviewer #2 and #3 both made the valid request that you should more clearly discuss possible reasons for why the basic manifold of neural activity was preserved, by comparing it to previous studies that have found "off-manifold" learning (delayed feedback, recalibration, training protocol – see reviews). In this context you should also discuss the limitations of having N=1. Although the reviewers agree that case studies in this scenario can be valuable, the limitations in generalization should be noted and discussed adequately.

3) All reviewers agreed also that the temporal change of the representational structure was a very interesting point. Please address the specific comments raised in the reviews for this analysis. As this is likely the most novel contribution of the paper, we believe it should be more explicitly mentioned in the abstract / introduction and more clearly discussed at the end of the paper.

Reviewer #1 (Recommendations for the authors):

I found the investigation very interesting and illuminating. Personally, I believe that the third result (i.e., the time-dependency of the representational structure) is one of the most interesting results of the paper, and should be mentioned in the abstract / introduction.

My only main comment concerns the first result – the comparison of the representational structure with previous fMRI results. Judging from Figure 1 – supp 1, the implantation site is located right on the boundary between the postcentral gyrus and posterior parietal cortex, seemingly ~5mm anterior to the peak of activation observed during grasping (AIP). Thus, based on the ROI definition used in both Ejaz et al., (2015) and Kieliba at al., (2021), it would be right on the boundary between SPLa and S1 – and cycto-architectonically, right on the boundary between BA2 and BA5. The SPLa ROI in Kieliba at al. extends substantially posterior to this. Thus, it is not clear whether the observed mismatch between neuronal recordings and fMRI reveals a real discrepancy (see Arbuckle et al., 2020), a function of the low reliability of the SPLa RDM (as discussed), or the mismatch in location. To improve this analysis, I would be happy to provide an RD from a more tightly matched ROI from the Ejaz et al., (2015) data, if needed.

Reviewer #2 (Recommendations for the authors):

Suggestions:

1) As described in the public review, the claim that motor representations are preserved after paralysis needs to be revised and described with more nuance to reflect the complex results that were presented in the manuscript. These revisions should be made to the title and throughout the discussion.

2) The third hypothesis (that pre-injury motor representations could have de-specialized) references Figure 2C, but that same figure panel was referenced for the second hypothesis (that the PC-IP representation would match that from SPLa in the imaging studies). It seems that the third hypothesis is not represented in Figure 2 so it should be added.

3) In figure 3a, it does appear that the single unit model became slightly less similar to the M1 model over time and in Figure 3c, the MDS plot suggests that the activity associated with each finger became slightly more variable (larger standard deviation) in the later sessions. While not a dramatic effect, could this be indicative of learning? The BCI performance should be shown for each session. Was there any trend in the overall success rate or error pattern over time?

4) The somatotopic model fit appeared to be stronger (in terms of model coefficient) than the muscle or usage model. Is this surprising given the interpretation that the somatotopic representation is primarily driven by sensory feedback, and this is a covert task. Might one expect the feedforward efferent signal to be stronger? Further, the discussion mentions that there is no topography in PPC, These seem to be in conflict. This may require a more nuanced discussion of the limitations of fMRI temporal resolution that may impact the ability to measure a consistent pattern. Further, the authors reference a previous study, Chivukula et al., as evidence of touch representations in PPC.

5) The discussion about limitations in the learning aspects of the study should be expanded as described in the weaknesses section of the public review. In addition, please comment on how the choice of the decoding window and features should impact the results given that there appears to be two temporally distinct representations measured in parietal cortex. Perhaps the participant converged on a strategy that was essentially a compromise between these two patterns, which both would have been included during the decoding window. Nevertheless, that no learning was observed in the current study is an important finding worthy of dissemination.

Reviewer #3 (Recommendations for the authors):

I'll avoid repeating concerns from the public review; most of those can be addressed by adding additional caveats or changing the introduction framing. As just one specific example, the manuscript would thus be substantially improved by appropriately narrowing the scope of claims. E.g., the last line of the Introduction would be only marginally longer, but would much better reflect the overall state of knowledge, if it were edited to say "Our results reveal that adult motor representations are preserved **in PPC** even after years without use".

Page 26 line 481: I really didn't understand this paragraph about estimating noise ceiling (and thus also I struggled to assess the claim that the model fit was "close to the theoretical maximum" on page 8 line 125). I understand at a high level what you're trying to do, but there are a lot of ways to estimate what even a perfect fit would look like given the single trial / spike measurement noise present. Can you explain the motivation/specific approach more (what do you mean by noise of the RDM?) and in more detail explain this calculation? Both more explanatory text, and adding the equations to Methods, would help. In general I found that the Methods could be more clear throughout and would recommend asking a colleague in a different systems neuroscience lab to read the Methods and explain what you did back to you (to reveal places where a reader would not be able to replicate the steps).

Page 26 line 491: Can you first briefly explain what the "directional hypothesis" is? Or else the reader would have had to read the Ejaz et al., methods to follow this section.

Can authors explain how hypothesis #1 (RSA structure matches natural statistics of use, as per ref 14) is distinct from hypothesis #2 (the BCI finger representation in PC-IP might instead match the representation of able-bodied individuals in the same brain area")? Couldn't #1 and #2 describe the same RSA structure (it seems that in fact they are highly correlated)? Was there past work showing they are different? Is there a reference to help unpack hypothesis #2? Right now hypothesis #1 vs #2 are really just illustrated by the two different RSA matrices in Figure 2b vs 2c, but that's asking the reader to squint and judge how different those matrices are. Some more hand-holding up front of where these hypotheses differ (and why) would help set the stage for this key question and result.

"In the first human BCI study of neural-population plasticity" (page 18, line 330) is a questionable claim (and unnecessarily provocative, in my view). For example, Wodlinger et al., (2015) examined BCI control of a robot arm over many sessions across ~250 days of use (Figure 5, and "showing that learning took place consistently over a long period of time") – this would seem to be an arm BCI learning study (asking, does neural activity change in a way that would lead to improved performance over time, by examining the decoder output which is, by definition, the task-relevant readout of the motor cortical population). That is not that different from how this study examines finger BCI learning. Another example is Eichenlaub, Jarosiewicz et al., (2020), which looked at learning-associated replay events in a BCI-controlled sequence task.

The page 18 line 331-332 sentence (about handwriting) also puzzled me: didn't that study explicitly find that yes, handwriting remains preserved through at least a decade of injury? Perhaps the present manuscript can be more specific about what future studies should ask about these complex motor skills

Essential revisions:

1) One important question to address clearly in the revision is the mismatch of the representational structure of the recorded PPC site with the SPLa ROI from the imaging studies (and the match with M1). Reviewer #1 raises the important questions that fMRI ROI is likely not equivalent to the recorded site and should be modified and the finding reinvestigated.

In the original submission, we found that PC-IP single neurons more strongly matched the MC fMRI RDMs than the SPLa fMRI RDMs, as measured by the whitened unbiased cosine (WUC) similarity (Diedrichsen et al., 2021). After narrowing the SPLa fMRI ROI to the implant region (PC-IP), a similar pattern of results still held. The PC-IP fMRI RDM matched the tetraplegic participant's data better than the task-optimal, unstructured model. Furthermore, participant NS's finger organization still matched the MC fMRI data more strongly than the PC-IP fMRI data from the same participants.

If the discrepancy persists after the reanalysis, Reviewer #2, and #3 make suggestions for a clearer discussion of the findings. Overall, the claim of the stability of pre-injury organization probably needs to be tempered a bit.

Given this discrepancy, we have tempered the preservation claims in the manuscript; we changed the title from "Preserved motor representations after paralysis" to "Stability of motor representations after paralysis." We comment on possible causes for the discrepancy as well:

"This result does obscure a straightforward interpretation of the RSA results – why does our recording area match MC better than the corresponding implant location? Several factors might contribute, including differing neurovascular sensitivity to the early and late phases of the neural response (Figure 4e), heterogeneous neural organizations across the single-neuron and voxel spatial scales (Arbuckle et al., 2020; Guest and Love, 2017; Kriegeskorte and Diedrichsen, 2016), or mismatches in functional anatomy between participant NS and standard atlases (Eickhoff et al., 2018)."

Still, we find it compelling that participant NS's finger representation strongly matches the MC fMRI RDM, which is known to reflect the natural hand use patterns of able-bodied individuals (Ejaz et al., 2015). Participant NS has not moved her own hands in 10 years, and BCIs remove any biomechanical constraints that might encourage correlated representations.

Revised text in the Discussion:

"Despite the lack of overt movement or biomechanical constraints, the measured finger representation still reflected these usage-related patterns."

Based on our RSA and learning analyses (Figures 2 and 3), such a structure is unlikely to have emerged after paralysis. The apparent preservation of the forward model (Figure 4e), as well as other recent studies of tetraplegic individuals, also corroborate an interpretation that motor representations are preserved. However, we acknowledge that single-neuron recordings in participants would provide the most direct comparison baseline, although these recordings are difficult to acquire.

Revised text in the Discussion:

"Comparisons with single-neuron recordings from able-bodied participants would validate this interpretation but may be difficult to acquire."

2) The stability of representations across learning is an interesting finding, but Reviewer #2 and #3 both made the valid request that you should more clearly discuss possible reasons for why the basic manifold of neural activity was preserved, by comparing it to previous studies that have found "off-manifold" learning (delayed feedback, recalibration, training protocol – see reviews). In this context you should also discuss the limitations of having N=1. Although the reviewers agree that case studies in this scenario can be valuable, the limitations in generalization should be noted and discussed adequately.

We agree with the reviewers that a different experimental paradigm could encourage off-manifold learning. We have expanded the discussion to clarify this point.

Revised text in the Discussion:

"The stability of finger representations here suggests that BCIs can benefit from the pre-existing, natural repertoire (Hwang et al., 2013), although learning can play an important role under different experimental constraints. In our study, the participant received only a delayed, discrete feedback signal after classification (Figure 1a). Because we were interested in understanding participant NS’s natural finger representation, we did not artificially perturb the BCI mapping. However, a previous study showed that BCI users can slowly learn to generate new, off-manifold neural activity patterns when the user receives continuous feedback and the BCI decoder is perturbed incrementally (Oby et al., 2019). Notably, learning was inconsistent when perturbations were sudden, indicating that learning is sensitive to specific guiding procedures. To further understand how much finger representations can be actively modified, future studies could benefit from perturbations (Kieliba et al., 2021; Oby et al., 2019), richer neurofeedback, and additional participants. "

2) The stability of representations across learning is an interesting finding, but Reviewer #2 and #3 both made the valid request that you should more clearly discuss possible reasons for why the basic manifold of neural activity was preserved, by comparing it to previous studies that have found "off-manifold" learning (delayed feedback, recalibration, training protocol – see reviews). In this context you should also discuss the limitations of having N=1. Although the reviewers agree that case studies in this scenario can be valuable, the limitations in generalization should be noted and discussed adequately.

3) All reviewers agreed also that the temporal change of the representational structure was a very interesting point. Please address the specific comments raised in the reviews for this analysis. As this is likely the most novel contribution of the paper, we believe it should be more explicitly mentioned in the abstract / introduction and more clearly discussed at the end of the paper.

We have expanded the abstract and introduction to further highlight the temporal analysis.

Revised text in the Abstract:

"Within individual BCI movements, the representational structure was dynamic, first resembling muscle activation patterns and then resembling the anticipated sensory consequences."

Revised text in the Introduction:

"Temporal dynamics provide another lens to investigate neural organization and its changes after paralysis. Temporal signatures can improve BCI classification (Willett et al., 2021) or provide a baseline for motor adaptation studies (Stavisky et al., 2017; Vyas et al., 2018). Notably, motor cortex activity exhibits quasi-oscillatory dynamics during arm reaching (Churchland et al., 2012). More generally, the temporal structure can depend on the movement type (Suresh et al., 2020) and the recorded brain region (Schaffelhofer and Scherberger, 2016). In this study, we recorded from the posterior parietal cortex (PPC), which is thought to compute an internal forward model for sensorimotor control (Desmurget and Grafton, 2000; Li et al., 2022; Mulliken et al., 2008; Wolpert et al., 1998). A forward model overcomes inherent sensory delays to enable fast control by predicting the upcoming states. If PPC activity resembles a forward model even after paralysis, this would suggest that even the temporal details of movement are preserved after injury."

“Furthermore, the neural representational dynamics matched the temporal profile expected of a forward model in able-bodied individuals, first resembling muscle activation patterns and then resembling expected sensory outcomes.”

Reviewer #1 (Recommendations for the authors):

I found the investigation very interesting and illuminating. Personally, I believe that the third result (i.e., the time-dependency of the representational structure) is one of the most interesting results of the paper, and should be mentioned in the abstract / introduction.

We have expanded the abstract and introduction to include this, as mentioned above in the response to the review summary.

My only main comment concerns the first result – the comparison of the representational structure with previous fMRI results. Judging from Figure 1 – supp 1, the implantation site is located right on the boundary between the postcentral gyrus and posterior parietal cortex, seemingly ~5mm anterior to the peak of activation observed during grasping (AIP). Thus, based on the ROI definition used in both Ejaz et al., (2015) and Kieliba at al. (2021), it would be right on the boundary between SPLa and S1 – and cycto-architectonically, right on the boundary between BA2 and BA5. The SPLa ROI in Kieliba at al. extends substantially posterior to this. Thus, it is not clear whether the observed mismatch between neuronal recordings and fMRI reveals a real discrepancy (see Arbuckle et al., 2020), a function of the low reliability of the SPLa RDM (as discussed), or the mismatch in location. To improve this analysis, I would be happy to provide an RD from a more tightly matched ROI from the Ejaz et al., (2015) data, if needed.

We have updated the analyses using anatomically matched fMRI data from Reviewer #1 (Jorn Diedrichsen) and from our original collaborator (Elena Amoruso).

Reviewer #2 (Recommendations for the authors):

Suggestions:

1) As described in the public review, the claim that motor representations are preserved after paralysis needs to be revised and described with more nuance to reflect the complex results that were presented in the manuscript. These revisions should be made to the title and throughout the discussion.

We have updated the manuscript as per the suggestion and as described above, under essential revisions.

2) The third hypothesis (that pre-injury motor representations could have de-specialized) references Figure 2C, but that same figure panel was referenced for the second hypothesis (that the PC-IP representation would match that from SPLa in the imaging studies). It seems that the third hypothesis is not represented in Figure 2 so it should be added.

This was an error in the subfigure indexing and has now been updated.

3) In figure 3a, it does appear that the single unit model became slightly less similar to the M1 model over time and in Figure 3c, the MDS plot suggests that the activity associated with each finger became slightly more variable (larger standard deviation) in the later sessions. While not a dramatic effect, could this be indicative of learning?

Overall, the decreases in distances contradict the expected learning hypothesis that participant NS would increase distances to perform the BCI task better. It appears that the decrease in inter-finger distances, while significant, was small enough to not affect classification.

Looking at the specific figures mentioned, the model fit metric (WUC similarity) used in Figure 3a requires a baseline to be compared against (Diedrichsen et al., 2021). Thus, the appropriate statistical test is to compare the slope of the fMRI model’s fit-values against the slope of the Unstructured model’s fit-values. The statistical tests do not find a difference in slopes, although we acknowledge that a small, significant trend could emerge if participant NS had participated in even more sessions.

The inter-finger distances decrease slightly in later sessions (Figure 3c, closer points in the MDS plot; Figure 3-—figure supplement 1, negative slope), although the effect size is rather small. As lower distances would generally indicate worse classification, this also runs counter to the learning hypothesis, which expects that participant NS should increase the neural distances to perform the BCI control task better.

Is it possible that learning took some form other than increased neural distances? One option would be that participant NS learned to control her neural variance; i.e. her attempted movements were more stereotyped, even, for example, as the means of her ring and pinky fingers moved closer together. However, the crossnobis (cross-validated Mahalanobis) distance already normalizes by the neural variance. In other words, lower-variance distributions appear farther away from each other.

The BCI performance should be shown for each session. Was there any trend in the overall success rate or error pattern over time?

We have added the BCI performance for each session in the new Figure 1—figure supplement 3. The evidence was not strong for an across-session trend in the online BCI accuracy (P = 0.060, BF = 1.7, slope = 0.009 / session).

4) The somatotopic model fit appeared to be stronger (in terms of model coefficient) than the muscle or usage model. Is this surprising given the interpretation that the somatotopic representation is primarily driven by sensory feedback, and this is a covert task. Might one expect the feedforward efferent signal to be stronger?

This is an interesting question. First, to clarify, we do not think that the somatotopic representation is driven by sensory feedback, but instead, represents the internal state estimate of the limb. State estimates, or internal abstractions that encode the brain’s estimate of the current position and velocity of the limb, are needed as sensory feedback is delayed, noisy, and multimodal. In motor intact individuals, this state estimate would combine sensory inputs with an internally produced estimate of the state as derived from known motor commands and a model of the body. As you mention, a loss of sensory inputs could potentially lead to a massive loss of neural representation. However, we know from the phantom limb literature that individuals can maintain internal representations of their lost limbs that they can still manipulate in their minds eye. This supports the idea that individuals maintain strong internal models of the body, even after catastrophic injury. Based on this prior evidence, it is not overly surprising that there exists a strong somatotopy representation. In terms of the relative strength of the two signals, we speculate that this will largely depend on brain region, e.g., we can imagine a region that represents the computed state estimate without any muscle-like model.

Further, the discussion mentions that there is no topography in PPC, These seem to be in conflict. This may require a more nuanced discussion of the limitations of fMRI temporal resolution that may impact the ability to measure a consistent pattern. Further, the authors reference a previous study, Chivukula et al., as evidence of touch representations in PPC.

Here it is important to distinguish cortical topography (e.g. a spatially organized mapping of the body onto the 2D cortical sheet, such as the homunculus) from a somatotopic representation (e.g., firing of single neurons as a smooth parametric function of body location.) Somatotopic representations can occur without a cortical topography if the somatotopic neurons are randomly distributed across the 2D cortical surface. This is consistent with our data in that representations of individual digits are all randomly intermingled within the same 4x4 mm patch of cortex that is sampled by our electrode array. Likewise, in Chivukula et al., tactile fields for all sensate parts of the body were also randomly intermingled in the 4x4 mm patch of cortex. Thus, even at the single neuron scale, we see no cortical topography, which may help to explain why no topography exists in the fMRI literature within PPC.

5) The discussion about limitations in the learning aspects of the study should be expanded as described in the weaknesses section of the public review.

This has been updated, as described above in essential revisions.

In addition, please comment on how the choice of the decoding window and features should impact the results given that there appears to be two temporally distinct representations measured in parietal cortex. Perhaps the participant converged on a strategy that was essentially a compromise between these two patterns, which both would have been included during the decoding window.

The input features were 1 firing rate/electrode, averaged over the 1-second decoding window. Knowing now that there appear to be two temporally distinct representations, it might have been interesting to split up the BCI decoding window into separate features through time. Perhaps this could help facilitate separate learning in the early and late components.

Revised text in the Discussion:

“Finally, given that finger representations were dynamic (Figure 4e), learning might occur separately in the early and late dynamics. Time-variant BCI decoding algorithms, such as recurrent neural networks (Sussillo et al., 2012; Willett et al., 2021), could help facilitate learning specific to different time windows of movement.”

Nevertheless, that no learning was observed in the current study is an important finding worthy of dissemination.

Reviewer #3 (Recommendations for the authors):

I'll avoid repeating concerns from the public review; most of those can be addressed by adding additional caveats or changing the introduction framing. As just one specific example, the manuscript would thus be substantially improved by appropriately narrowing the scope of claims. E.g., the last line of the Introduction would be only marginally longer, but would much better reflect the overall state of knowledge, if it were edited to say "Our results reveal that adult motor representations are preserved **in PPC** even after years without use".

This has been updated, as described above in essential revisions.

Page 26 line 481: I really didn't understand this paragraph about estimating noise ceiling (and thus also I struggled to assess the claim that the model fit was "close to the theoretical maximum" on page 8 line 125). I understand at a high level what you're trying to do, but there are a lot of ways to estimate what even a perfect fit would look like given the single trial / spike measurement noise present. Can you explain the motivation/specific approach more (what do you mean by noise of the RDM?) and in more detail explain this calculation? Both more explanatory text, and adding the equations to Methods, would help.

We agree that many different sources that could contribute to variability in neuroscience results. RDM noise ceilings usually only measure variability across RDMs (i.e., variability across sessions or subjects), not single-trial/spike measurement noise. Trial-to-trial variability (including spike measurement noise) is accounted for by the use of the cross-validated Mahalanobis distance, which normalizes single-neuron firing rates by an estimate of their trial-to-trial variability (i.e., noise covariance). Thus, we are mainly concerned with the variability across repeated sessions or subjects.

Revised Methods section:

“Noise ceiling for model fits

Measurement noise and behavioral variability cause data RDMs to vary across repetitions, so even a perfect model RDM would not achieve a WUC similarity of 1. To estimate the noise ceiling(Nili et al., 2014) (the maximum similarity possible given the observed variability between data RDMs), we assume that the unknown, perfect model resembles the average RDM. Specifically, we calculated the average similarity of each individual-session RDM (Figure 2—figure supplement 2) with the mean RDM across all other sessions (i.e., excluding that session):

where is the WUC similarity function, is the number of RDMs, refers to a single RDM from an individual session, and  is the “lower” noise ceiling. This noise ceiling is analogous to leave-one-out-cross-validation. If a model achieves the noise ceiling, the model fits the data well (Nili et al., 2014).”

In general I found that the Methods could be more clear throughout and would recommend asking a colleague in a different systems neuroscience lab to read the Methods and explain what you did back to you (to reveal places where a reader would not be able to replicate the steps).

Great suggestion. We followed this advice and hope the revised Methods are more clear.

Page 26 line 491: Can you first briefly explain what the "directional hypothesis" is? Or else the reader would have had to read the Ejaz et al., methods to follow this section.

We originally used the term “directional hypothesis” to refer to one-sided t-tests; we have revised our terminology to make this explicit.

Revised text in Methods:

“To calculate Bayes factors for one-sided t-tests (for example, ), we sampled (N = 106) from the posterior of the corresponding two-sided t-test (), calculated the proportion of samples that satisfied the one-sided inequality, and divided by the prior odds (Morey and Wagenmakers, 2014) ( ).”

Can authors explain how hypothesis #1 (RSA structure matches natural statistics of use, as per ref 14) is distinct from hypothesis #2 (the BCI finger representation in PC-IP might instead match the representation of able-bodied individuals in the same brain area”)? Couldn’t #1 and #2 describe the same RSA structure (it seems that in fact they are highly correlated)? Was there past work showing they are different? Is there a reference to help unpack hypothesis #2? Right now hypothesis #1 vs #2 are really just illustrated by the two different RSA matrices in Figure 2b vs 2c, but that’s asking the reader to squint and judge how different those matrices are. Some more hand-holding up front of where these hypotheses differ (and why) would help set the stage for this key question and result.

We agree that, a priori, it is possible hypotheses #1 and #2 could describe the same RSA structure. Indeed, at a high level, the PC-IP and motor cortex (MC) fMRI finger RDMs are more similar than they are different. To make this less confusing to readers, we have revised the way we introduce the hypotheses, deferring the previous hypothesis #1 to later.

Revised text in the Results:

“We used RSA to test three hypotheses: (1) the BCI finger representational structure could match that of able-bodied individuals (Ejaz et al., 2015; Kieliba et al., 2021) (Figure 2b and Figure2—figure supplement 1), which would imply that motor representations did not reorganize after paralysis… We note that our able-bodied model was recorded from human PC-IP using fMRI, which measures fundamentally different features (millimeter-scale blood oxygenation) than microelectrode arrays (sparse sampling of single neurons). Another possibility is that (2) the participant’s pre-injury motor representations had de-specialized after paralysis, such that finger activity patterns are unstructured and pairwise-independent…”

Only later do we compare to fMRI recordings from another region (i.e., MC). We also have revised the text to explain the quantitative differences between the PC-IP and MC fMRI finger RDMs.

Revised text in the Results:

“We also compared the PC-IP BCI RDM with able-bodied fMRI motor cortex (MC) RDMs, which have been previously shown to match the patterns of natural hand use (Ejaz et al., 2015). The able-bodied MC and PC-IP fMRI finger organizations are similar in that they represent the thumb distinctly from the other fingers, but PC-IP represents each of the non-thumb fingers similarly while MC distinguishes between all five fingers (Figure 2—figure supplement 1).”

Revised text in the Discussion:

“Compared to the PC-IP fMRI finger representation, MC represents the non-thumb fingers more distinctly from each other (Figure 2—figure supplement 1).”

“In the first human BCI study of neural-population plasticity” (page 18, line 330) is a questionable claim (and unnecessarily provocative, in my view). For example, Wodlinger et al., (2015) examined BCI control of a robot arm over many sessions across ~250 days of use (Figure 5, and "showing that learning took place consistently over a long period of time") – this would seem to be an arm BCI learning study (asking, does neural activity change in a way that would lead to improved performance over time, by examining the decoder output which is, by definition, the task-relevant readout of the motor cortical population). That is not that different from how this study examines finger BCI learning. Another example is Eichenlaub, Jarosiewicz et al., (2020), which looked at learning-associated replay events in a BCI-controlled sequence task.

Removed as recommended.

The page 18 line 331-332 sentence (about handwriting) also puzzled me: didn't that study explicitly find that yes, handwriting remains preserved through at least a decade of injury? Perhaps the present manuscript can be more specific about what future studies should ask about these complex motor skills

We have updated the discussion to clarify the point we were trying to make:

“As BCIs enable more complex motor skills, such as handwriting (Willett et al., 2021), future studies could investigate whether these complex skills also retain their pre-injury representational structure. For example, does a tetraplegic participant’s BCI handwriting look like their physical, pre-injury handwriting? These results will have important implications for the design of future neural prosthetics.”