The Neural Representation of Force across Grasp Types in Motor Cortex of Humans with Tetraplegia

Abstract Intracortical brain-computer interfaces (iBCIs) have the potential to restore hand grasping and object interaction to individuals with tetraplegia. Optimal grasping and object interaction require simultaneous production of both force and grasp outputs. However, since overlapping neural populations are modulated by both parameters, grasp type could affect how well forces are decoded from motor cortex in a closed-loop force iBCI. Therefore, this work quantified the neural representation and offline decoding performance of discrete hand grasps and force levels in two human participants with tetraplegia. Participants attempted to produce three discrete forces (light, medium, hard) using up to five hand grasp configurations. A two-way Welch ANOVA was implemented on multiunit neural features to assess their modulation to force and grasp. Demixed principal component analysis (dPCA) was used to assess for population-level tuning to force and grasp and to predict these parameters from neural activity. Three major findings emerged from this work: (1) force information was neurally represented and could be decoded across multiple hand grasps (and, in one participant, across attempted elbow extension as well); (2) grasp type affected force representation within multiunit neural features and offline force classification accuracy; and (3) grasp was classified more accurately and had greater population-level representation than force. These findings suggest that force and grasp have both independent and interacting representations within cortex, and that incorporating force control into real-time iBCI systems is feasible across multiple hand grasps if the decoder also accounts for grasp type.

While kinematic iBCIs can restore basic reaching and grasping movements, restoring the ability to grasp and interact with objects requires both kinematic and kinetic (force-related) information (Chib et al., 2009;Flint et al., 2014;Casadio et al., 2015). Specifically, sufficient contact force is required to prevent object slippage; however, excessive force may cause mechanical damage to objects (Westling and Johansson, 1984). Therefore, introducing force calibration capabilities during grasp control would enable iBCI users to perform more functional tasks.
Intended forces are usually produced in the context of task-related factors, including grasp postures used to generate forces (Murphy et al., 2016). The representation and decoding of grasps, independent of forces, has been studied extensively in primates Vargas-Irwin et al., 2010;Carpaneto et al., 2011;Townsend et al., 2011;Hao et al., 2014;Schaffelhofer et al., 2015) and humans (Pistohl et al., 2012;Chestek et al., 2013;Bleichner et al., 2014Bleichner et al., , 2016Klaes et al., 2015;Leo et al., 2016;Branco et al., 2017). Importantly, previous studies suggest that force and grasp are encoded by overlapping populations of neural activity (Sergio and Kalaska, 1998;Carmena et al., 2003;Sergio et al., 2005;Milekovic et al., 2015;Sburlea and Müller-Putz, 2018). While some studies suggest that force is encoded at a high level independent of motion and grasp (Chib et al., 2009;Hendrix et al., 2009;Pistohl et al., 2012;Intveld et al., 2018), others suggest that it is encoded at a low level intertwined with grasp (Hepp-Reymond et al., 1999;Degenhart et al., 2011). Thus, the degree to which intended hand grasps and forces interact within the neural space, and how such interactions affect force decoding performance, remain unclear. To our knowledge, these scientific questions have not been explored in individuals with tetraplegia, who constitute a target population for iBCI technologies.
To answer these questions, we characterized the extent to which three discrete, attempted forces were neurally represented and offline-decoded across up to five hand grasp configurations in two individuals with tetraplegia. Our results suggest that force has both grasp-independent and grasp-dependent (interacting) representation in motor cortex. Additionally, while this study demonstrates the feasibility of incorporating discrete force control into human-operated iBCIs, these systems will likely need to incorporate grasp and other task parameters to achieve optimal performance.

Study permissions and participants
Study procedures were approved by the United States Food and Drug Administration (Investigational Device Exemption #G090003) and the Institutional Review Boards of University Hospitals Case Medical Center (protocol #04-12-17), Massachusetts General Hospital (2011P001036), the Providence VA Medical Center (2011-009), Brown University (0809992560), and Stanford University (protocol #20 804). Human participants were enrolled in the BrainGate2 Pilot Clinical Trial (ClinicalTrials.gov number NCT00912041). Informed consent, including consent to publish, was obtained from the participants before their enrollment in the study.
This study includes data from two participants with chronic tetraplegia. Both participants had two, 96-channel microelectrode intracortical arrays (1.5-mm electrode length, Blackrock Microsystems) implanted in the hand and arm area ("hand knob"; Yousry et al., 1997) of dominant motor cortex. Participant T8 was a 53-year-old righthanded male with C4-level AIS-A spinal cord injury eight years before implant, and T5 was a 63-year-old right-handed male with C4-level AIS-C spinal cord injury. More surgical details can be found at (Ajiboye et al., 2017) for T8 and (Nuyujukian et al., 2018) for T5.

Participant task
The goal of this study was to measure the degree to which various hand grasps affect decoding of grasp force from motor cortical spiking activity. To this end, participants T8 and T5 took part in several research sessions in which they attempted to produce three discrete squeeze forces (light, medium, hard) using one of four designated hand grasps (closed pinch, open pinch, ring pinch, power). Squeeze force, defined here as the amount of force needed to deform an object, is distinct from grip force, which is the amount of force needed to grasp an object of particular weight and friction (Westling and Johansson, 1984). In this study, participants were instructed to produce squeeze forces as opposed to grip forces. This was because participants could not receive somatosensory feedback about the object properties that usually inform grip force production, yet retained the capacity to emulate squeeze forces in response to audio and visual cues.
The four hand grasps used to emulate squeeze forces were chosen to study force representation within multiple grasp-related contexts. The open and closed pinch grasps were included to determine how forces were represented when emulated with grasps of similar function (thumb-index precision grasp) but different posture. The ring pinch grasp was included to determine the effects of using different fingers to produce similar forces. Finally, the power grasp was included to determine the influence of power versus precision grasping on force representation. Participant T8 completed six sessions between trial days 735-956 relative to the date of his microelectrode array placement surgery; and T5 completed one session on trial day 390. During session 5, participant T8 emulated discrete forces using attempted elbow extension in addition to the four distal hand grasps. This enabled the study of force representation across the entire upper limb. Table  1 lists all relevant sessions and their associated task parameters.
Each research session consisted of multiple 4-min data collection blocks, which were each assigned to a particular hand grasp or elbow movement, as illustrated in Figure  1B. Blocks were presented in a pseudorandom order, in which hand grasps were assigned randomly to each set of two (session 1), four (sessions 2-4, 6-7), or five (session 5) blocks. This allowed for an equal number of blocks per hand grasp, distributed evenly across the entire research session.
All blocks consisted of ;20 trials, which were presented in a pseudorandom order by repeatedly cycling through a complete, randomized set of force levels until the end of the block. During each trial, participants used  Rastogi et al., 2020). Participants had two 96-channel microelectrode arrays placed chronically in motor cortex, which recorded neural activity while participants completed a force task. TC and SBP features were extracted from these recordings. Figure 1A is reprinted by permission from Springer Nature as indicated in the Terms and Conditions of a Creative Commons Attribution 4.0 International license (https://www.nature. com/srep/). B, Research session architecture. Each session consisted of 12-21 blocks, each of which contained ;20 trials (see Table 1). In each trial, participants attempted to generate one of three visually-cued forces with one of four grasps: power, closed pinch, open pinch, ring pinch. During session 5, participant T8 also attempted force production using elbow extension. Each trial contained a preparatory (prep) phase, a go phase where forces were actively embodied, and a stop phase where neural activity was allowed to return to baseline. Participants were prompted with both audio and visual cues, in which a researcher squeezed or lifted an object associated with each force level. During pinch blocks, the researcher squeezed the pinchable objects (cotton ball, eraser, nasal aspirator tip) using the particular pinch grip dictated by the block (ring pinch, open pinch, closed pinch). Here, only closed pinches of objects are shown. kinesthetic imagery (Stevens, 2005;Mizuguchi et al., 2017) to internally emulate one of three discrete force levels, or rest, with the dominant hand. Participants received simultaneous audio and visual cues indicating which force to produce, when to produce it, and when to relax. Participants were visually cued by observing a researcher squeeze one of nine graspable objects corresponding to light, medium, and hard forces (no object was squeezed during "rest" trials), as shown in Figure 1B. The participants were asked to "follow along" and attempt the same movements that the researcher was demonstrating. The graspable objects were grouped into three sets of three, corresponding to forces embodied using a power grasp (sponge = light, stress ball = medium, tennis ball = hard), a pincer grasp (cotton ball = light, nasal aspirator tip = medium, eraser = hard), or elbow extension (5-lb dumbbell = light, 10-lb dumbbell = medium, 15-lb dumbbell = hard). These visual cues, which were included to make the concept of light, medium, and hard forces seem less abstract to participants after years of deafferentation, were deemed unlikely to be reflected within the go-phase neural response based on previous investigations (Rastogi et al., 2020). Objects were chosen to be of similar weight, size, and shape to minimize the effects of visual confounds within the neural data.
During the prep phase, which lasted a pseudo-randomly determined period between 2.7 and 3.3 s to reduce confounding effects from anticipatory activity, the researcher presented an object indicating the force level to be attempted. The researcher then squeezed the object (or lifted the object, in the case of elbow extension) during the go phase (3-5 s), and finally released the object at the beginning of the stop phase (5 s). When squeezing (or lifting) objects, the researcher used the grasp type dictated by the block. For example, to visually cue hard forces, the researcher used a ring pinch to squeeze the eraser during ring pinch blocks, but used an open pinch grasp to squeeze the eraser during open pinch blocks.

Neural recordings Preprocessing
In both participants, each intracortical microelectrode array was attached to a percutaneous pedestal connector on the head. A Blackrock shielded Patient Cable connected the pedestals to front-end amplifiers and a NeuroPort System (Blackrock Microsystems) that bandpass filtered (0.3 Hz to 7.5 kHz) and digitized (30 kHz) the neural signals from each channel on the microelectrode array. These digitized signals were preprocessed in Simulink using the xPC real-time operating system (The MathWorks Inc.). Each channel was bandpass (BP) filtered (250-5000 Hz), common average referenced (CAR), and down-sampled to 15 kHz in real time. CAR was implemented by selecting 60 channels from each microelectrode array that exhibited the lowest variance, and then averaging these channels together to yield an array-specific CAR. This reference signal was subtracted from the signals from all channels within each of the arrays.

Extraction of neural features
From each filtered, CAR channel, two neural features were extracted in real time using the xPC operating system from non-overlapping 20-ms time bins. These features, as illustrated in Figure 1A, included unsorted threshold crossing (TC) rate and spike band power (SBP) features. Each TC feature, which was equivalent to multiunit activity  was defined as the number of times a channel's recorded voltage time series crossed a predefined noise threshold [À4.5 Â root mean square (RMS) voltage], divided by the width of the time bin (Christie et al., 2015) . The RMS voltage on each channel was calculated from 1 min of neural data recorded at the beginning of each research session. Additionally, each SBP feature was computed as the average signal power of the spike band (250-5000 Hz) within each time bin. Thus, SBP features were computed in the same manner as local field potentials (LFPs) and EEG signal power bands.
These calculations yielded 384 neural features per participant, which were used for offline analysis without spike sorting (Trautmann et al., 2019). TC features were labeled from 1 to 192 according to the recording electrodes from which they were extracted. Corresponding SBP features were labeled from 193 to 384. All features were normalized by subtracting the block-specific mean activity of the features within each recording block, to minimize nonstationarities in the data.
Unless otherwise stated, all subsequent offline analyses of neural data were performed using MATLAB software within a Windows 64-bit operating system.

Characterization of individual neural feature tuning
The first goal of this study was to determine the degree to which force-related and grasp-related information are represented within individual TC and SBP neural features. Specifically, neural activity resulting from three discrete forces and two (session 1), four (sessions 2-4, 6-7), or five (session 5) grasps, resulted in 6, 12, or 15 conditions of interest per session, respectively. See Table 1 for a list of grasps included for each individual research session. To visualize individual feature responses to force and grasp, each feature's peristimulus time histogram (PSTH) was computed for each of these conditions by averaging the neural activity over go-cue-aligned trials. These trials were temporally smoothed with a Gaussian kernel (100ms SD) to aid in visualization.
To determine how many of these individual features were tuned to force and/or grasp, statistical analyses were implemented in MATLAB and with the WRS2 library in the R programming language (Wilcox, 2017), as in Rastogi et al. (2020). Briefly, features were preprocessed in MATLAB to compute each feature's mean go-phase deviation from baseline during each trial. Baseline activity was computed by averaging neural activity across multiple rest trials.
In R, the distribution of go-phase neural deviations was found to be normal (analysis of Q-Q plots and Shapiro-Wilk tests, p , 0.05) but heteroskedastic (Levene's test, p , 0.05), necessitating a two-way Welch ANOVA analysis to determine neural tuning to force, grasp, and their interaction (p , 0.05). Features exhibiting an interaction between force and grasp were further separated into individual grasp conditions (closed pinch, open pinch, ring pinch, power, elbow), within which one-way Welch-ANOVA tests were implemented to find interacting features that were tuned to force. All p values were corrected for multiple comparisons using the Benjamini-Hochberg procedure (Benjamini and Hochberg, 1995).

Neural population analysis and decoding
The second goal of this study was to determine the degree to which force and grasp are represented within, and can be decoded from, the neural population. Here, the neural population was represented using both traditional and demixed principal component analysis (dPCA).

Visualizing force representation with traditional PCA
In order to visualize how consistently forces were represented across different grasps, neural activity collected during individual sessions were graphically represented within a low-dimensional space found using PCA. Notably, during session 5, participant T8 attempted to produce three discrete forces not only with several grasps, but also with an elbow extension movement. Therefore, two sets of PCA analyses were implemented on the data. The first, which was applied to all sessions, performed PCA on all force and grasp conditions within the session. In the second analysis specific to session 5 only, PCA was applied solely on power grasping and elbow extension trials to elucidate whether forces were represented in a consistent way across the entire upper limb. For both analyses, the PCA algorithm was applied to neural feature activity that was averaged over multiple trials and across the go phase of the task.
The results of each decomposition were plotted in a low-dimensional space defined by the first two principal components (PCs). The force axis within this space, given by Equation 1, was estimated by applying multiclass linear discriminant analysis (LDA; Juric, 2020) to the centered, force-labeled PCA data, and then using the largest LDA eigenvector as the multidimensional slope m of the force axis. Here, PC score is the principal component score, or representation of the neural data in PCA space, and f is the intended force level. A consistent force axis across multiple grasps within PCA space would suggest that forces are represented in an abstract (and thus grasp-independent) manner.
PC score ¼ mf: (1) dPCA The remainder of population-level analysis was implemented using dPCA. dPCA is a dimensionality reduction technique that, similarly to traditional PCA, compresses neural population activity into a few components that capture the majority of variance in the source data (Kobak et al., 2016). Unlike traditional PCA, which yields PCs that each capture signal variance from multiple parameters of interest, dPCA performs an ANOVA-like decomposition of data into task-dependent dimensions of neural activity. That is, the resulting dPCs are tuned to individual task parameters; thus, they are much easier to interpret than traditional PCs. Additionally, because dPCA performs an ANOVA-like decomposition of data, it serves as a population-level analog to the two-way Welch ANOVA analysis implemented on individual neural features.
Briefly, the matrix X of neural data is decomposed into trial-averaged neural activity explained by time (t), various task parameters (p 1 , p 2 ), their interaction (p 1 p 2 ), and noise, according to Equation 2. Next, dPCA finds separate decoder (D) and encoder (E) matrices for each marginalization M by minimizing the loss function L exhibited in Equation 3.
The resulting dPCs, obtained by multiplying the neural data X by the rows of each decoder matrix D M , are, in theory, de-mixed, in that the variance explained by each component is because of a single, specific task parameter M. These dimensions of neural activity not only reveal population-level trends in neural data, but they can also be used to decode task parameters of interest. Critically, dPCA can be used to decode task parameters from the data while still preserving its original geometry. Thus, a single technique can be used to analyze the underlying structure of neural data as it relates to the encoding of task parameters, and to simultaneously quantify how well these parameters can be decoded for use in an iBCI system (Kobak et al., 2016).

Single dPCA component implementation
In the present study, the task parameters of interest were force and grasp. Here, one goal was to use variance as a metric to quantify the degree to which force and grasp were represented within the neural population as a whole. Therefore, for each research session listed in Table 1, the neural data X was temporally smoothed using a Gaussian filter (100-ms SD) and decomposed into neural activity that varied with four marginalizations X M , as per Equation 2: time (condition independent), force, grasp, and an interaction between force and grasp. The variance that each marginalization accounted for was computed as the sum of squares of the mean-centered neural data contained within the marginalization.
An additional goal was to isolate neural components that contained useful information about force and grasp, i.e., components that would enable discrimination between individual force levels and grasp types. First, dPCA was used to reduce each of the four, 384-dimensional, mean-centered marginalizations X M into 20 dPCs, as described by Equation 3. This yielded 80 dPCs across all four marginalizations. Second, the variances accounted for by each of the 80 components were computed as the sum of squares. Third, the top 20 out of 80 components with the highest variance were selected as representing the majority of variance in the neural dataset and were assembled into a decoder matrix D. Finally, each of these top 20 components was assigned to one of the four marginalizations of interest according to the marginalization from which it was extracted. For example, dPCs that were extracted from the force marginalization X force were deemed as force-tuned dPCs; those extracted from the grasp marginalization X grasp were deemed as grasp tuned dPCs, and those extracted from the marginalization X F/G representing an interaction between force and grasp were deemed as interacting dPCs.
Each dPC's information content was further quantified in two ways. First, to assess the degree to which dPCs were demixed, each dPC's variance was subdivided into four sources of variance corresponding to each of the four marginalizations of interest, as per Equation 2. Second, the decoder axis associated with each dPC was used as a linear classifier to decode intended parameters of interest. Specifically, each force-tuned dPC was used to decode force at every time point of the behavioral task, while each grasp-tuned dPC was used to decode grasp, but not force. Likewise, components that exhibited an interaction between force and grasp were used to decode forcegrasp pairs. Condition-independent dPCs, which were tuned to time, were not used to decode force or grasp from the neural activity.
Linear classification was implemented using 100 iterations of stratified Monte Carlo leave-group-out cross-validation (Kobak et al., 2016). During each iteration, one random group of F x G test "pseudo-trials," each corresponding to one of the several force-grasp conditions, was set aside during each time point (F = number of intended forces, G = number of intended grasps). Next, dPCA was implemented on the remaining trials, and the decoder axes of the resulting dPCs were used to predict the intended forces or intended grasps indicated by the test set of pseudo-trials at each time point. This was accomplished by first computing mean dPC values for each force-grasp condition, separately for each time point; projecting the F x G pseudo-trials onto the decoder axes of the dPCs at each time point; and then classifying the pseudo-trials according to the closest class mean (Kobak et al., 2016). The proportion of F x G pseudo-trials correctly classified across 100 iterations at each time point constituted a time-dependent classification accuracy. Chance performance was computed by performing 100 shuffles of all available trials, randomly assigning force or grasp conditions to the shuffled data, and then performing the same cross-validated classification procedure within each of the 100 shuffles. Classification accuracies that exceeded the upper range of chance performance were deemed significant.

Force and grasp decoding using multiple dPCs
Two additional goals of this study were to determine whether intended forces could be accurately predicted from neural population data and whether these predictions depended on hand grasp configuration. To this end, dPCs that were tuned to force, grasp, and an interaction between force and grasp were used to construct multidimensional force and grasp decoders within each session. Specifically, the force decoder was constructed by combining the decoding axes of force-tuned and interacting components into a single, multidimensional decoder D F ; likewise, the grasp decoder D G was constructed by combining the decoding axes of grasp-tuned and interacting components.
Each of these decoders was used to perform 40 runs of linear force and grasp classification for each of S research sessions per participant, implemented using the aforementioned stratified Monte Carlo leave-group-out crossvalidation procedure (S = 6 for T8; S = 1 for T5). As in the single component implementation (Kobak et al., 2016), each run was accomplished in multiple steps. First, the mean values of all dPCs included within the multidimensional decoder were computed for each force-grasp condition, separately for each time point. Second, at each time point, the F x G pseudo-trials were projected onto the multidimensional decoder axis and classified according to the closest class mean. The proportion of test trials correctly classified at each time point over 100 iterations constituted a time-dependent force or grasp classification accuracy.
The aforementioned computations yielded 40 Â S time-dependent force and grasp classification accuracies per participant. Session-averaged, time-dependent force and grasp classification accuracies were computed by averaging the performance over 240 sessionruns for participant T8 (40 runs Â six sessions) and 40 session-runs for participant T5 (40 runs Â one session). These averages were compared with chance performance, which was computed by performing 100 shuffles of all trials, randomly assigning force or grasp conditions to the shuffled data, and then performing force and grasp classification on each of the shuffled datasets using the multidimensional decoders D F and D G . Time points when force or grasp classification exceeded the upper bound of chance were deemed to contain significant force-related or grasp-related information.
To visualize the degree to which individual forces and grasps could be discriminated, confusion matrices were computed over go-phase time windows during which the neural population contained significant force-related and grasp-related information. The time window began when session-averaged, time-dependent classification accuracy exceeded 90% of maximum achieved performance within the go phase, and ended at the end of the go phase. First, classification accuracies for each of the S Â 40 sessionruns were approximated by averaging classification performance across the prespecified go-phase time window. These time-averaged accuracies, which are henceforth referred to as mean force and grasp accuracies, were next averaged over all S Â 40 session-runs to yield confusion matrix data. In this way, confusion matrices were computed to visualize force-related discriminability across all trials, force-related discriminability within individual grasp types, and grasp-related discriminability across all trials.
Classification performances for individual forces and individual grasps were statistically compared using parametric tests implemented on mean force and grasp accuracies. Specifically, for each participant, mean classification accuracies for light, medium, and hard forces were compared by implementing one-way ANOVA across mean force accuracies from all S Â 40 session runs. The resulting p values were corrected for multiple comparisons using the Benjamini-Hochberg procedure (Benjamini and Hochberg, 1995). Likewise, mean classification accuracies for closed pinch, open pinch, ring pinch, power, and elbow "grasps" were compared by implementing one-way ANOVA across all mean grasp accuracies. These comparisons were implemented to determine whether offline force and grasp decoding yielded similar versus different classification results across multiple forces and multiple grasps.
Statistical analysis was also used to determine the degree to which grasp affected force decoding accuracy. This was achieved by implementing two-way ANOVA on mean force accuracies that were labeled with the grasps used to emulate these forces. The results of the two-way ANOVA showed a statistically significant interaction between force and grasp. Therefore, the presence of simple main effects was assessed within each force level and within each grasp type. Specifically, one-way ANOVA was implemented on mean accuracies within individual force levels to determine whether light forces, for example, were classified with similar degrees of accuracies across all grasp types. Similarly, one-way ANOVA was implemented on mean accuracies within individual grasps to see whether intended forces affected grasp classification accuracy; p values resulting from these analyses were corrected for multiple comparisons using the Benjamini-Hochberg procedure.
Finally, this study evaluated how well dPCA force decoders could generalize to novel grasp datasets in T8 session 5 and T5 session 7. Specifically, within each session, a multidimensional force decoder D F was trained on neural data generated during all but one grasp type, and then its performance was evaluated on the attempted forces emulated using the left-out "novel" grasp. To establish the generalizability of force decoding performance across many novel grasps, this analysis cycled through all available grasps attempted during session 5 (closed pinch, open pinch, ring pinch, power, elbow extension) and session 7 (closed pinch, open pinch, ring pinch, power). For each novel grasp, the trained decoder D F was used to perform 40 runs of stratified Monte Carlo leave-group-out cross-validated linear force classification on two sets of test data: the "initial grasp" dataset, which originated from the grasps on which the force decoder was trained; the novel grasp dataset, which originated from the leave-out test grasp. The resulting time-dependent, "initial grasp" and "novel grasp" decoding performances from the go-phase time window during above-90% maximum classification accuracy were averaged over 40 runs, and then compared using a standard t test. P values resulting from the statistical analysis were corrected for multiple comparisons across forces and test grasps using the Benjamini-Hochberg procedure.

Comparison of force encoding models
The overarching goal of this study, which is to determine the extent to which force representation within motor cortex depends on grasp, arose from two conflicting hypotheses indicating that force representation is either grasp-independent (Chib et al., 2009;Hendrix et al., 2009;Pistohl et al., 2012;Intveld et al., 2018) or grasp-dependent (Hepp-Reymond et al., 1999;Degenhart et al., 2011). The grasp-independent and grasp-dependent force encoding hypotheses can be mathematically modeled as per Equations 4, 5, respectively: In these equations, x ij is an N Â T Â TR matrix of neural activity generated within N neural features over T time points, during TR trials of a particular grasp i and force j. The term g i is an N Â T Â TR matrix of baseline feature activity during the grasp i, f is an N Â T Â TR matrix of baseline activity feature activity during force generation, and s j is a discrete scalar force level. Finally, the coefficients a, b, c, and d are constants. In Equation 4, the overall neural activity x ij consists of an addition of independent force-related and grasp-related terms, as is thus referred to as the additive model of force encoding. In contrast, Equation 5 models the neural activity x ij as a multiplication of the force level s j with baseline grasp activity g i , and is hence referred to as the scalar model of force encoding.
An additional model, indicated by Equation 6, incorporates terms from both the additive and scalar models of force encoding, and is thus referred to as the combined model: x ij ¼ ag i 1 bfs j 1 cs j g i 1 d The additive (grasp-independent) and scalar (grasp-dependent) hypotheses of force encoding were graphically illustrated with a toy example of expected grasp-independent versus grasp-dependent (interacting) representations of force within the neural space. In the toy example, the model coefficients a, b, and c were set to one, and the model coefficient d was set to zero. The neural activity x ij was a vector of trial-averaged activity from 100 simulated neural features during a single time point, generated during a particular grasp i and force j. The variable g i was a 100 Â 1 vector of normalized baseline feature activity during the grasp i, f was a 100 Â 1 vector of normalized baseline neural feature activity during force generation, and s j was a discrete, scalar force level (1, 2, or 5). The vectors g i and f contained values drawn from the standard normal distribution.
Additionally, cross-validated ordinary least squares regression was used to quantify the degree to which the additive, scalar, and combined models explained the neural data x ij recorded from participants T8 and T5. Here, the neural data x ij consisted solely of force, grasp, and interacting components; condition-independent components of x ij were omitted. Thus, the matrix x ij was computed by compressing the 384-dimensional neural feature data using the dPCA decoder matrix D, eliminating CI-tuned dPCs, and then transforming the data back to feature space using the encoder matrix E (see Eq. 3). Baseline grasp activity g i was estimated by isolating grasp-tuned components from the neural data, transforming these components back to feature space using the encoder matrix E grasp , and then averaging the resulting activity over force conditions. Similarly, the baseline force activity f was estimated by isolating force-tuned dPCs, transforming these components back to feature space using the encoder matrix E force , and then averaging the resulting data over all force-grasp conditions. All three neural activity variables consisted of 384 Â T Â TR matrices, where T was the number of go-phase time points and TR was the number of trials emulated with an individual force-grasp combination. As in the toy example, s j was a discrete scalar force level (1, 2, or 5).
Cross-validated regression analysis for each model was performed using 100 iterations of a stratified Monte Carlo leave-group-out scheme. Notably, the regression was performed on the data x generated during all combinations of forces and grasps, as opposed to x ij , generated during a particular grasp i and force j. During each iteration of cross-validation, one random group of F x G pseudo-trials, each corresponding to one of the several force-grasp conditions, was set aside as a test dataset. Next, model coefficients were trained via ordinary least squares regression on the remaining data. Finally, the trained model was used to predict the neural activity generated during the emulated pseudo-trials, resulting in an R 2 value for each iteration. The distributions of R 2 values generated from each model were statistically compared by implementing a multiple comparisons test (Tukey method) on the results of a one-way ANOVA analysis.

Data and code accessibility
This study made use of several computational algorithms implemented using publicly available source code packages. Code for the WRS2 R package, which was used to characterize single features, is available at https:// CRAN.R-project.org/package=WRS2. Source code for the dPCA algorithm can be implemented either in MATLAB or Python and is available at https://github.com/machenslab/ dPCA. The dPCA source code was modified to perform multidimensional decoding of force and grasp; these modified scripts can be made available on reasonable request by contacting the lead or senior authors. Finally, MATLAB code for the multiclass LDA algorithm used to compute low-dimensional force axes within PCA space is available on the MATLAB file exchange at https://www.mathworks. com/matlabcentral/fileexchange/31760-multiclass-lda.
The data presented in this study can be made available on reasonable request by contacting the lead or senior authors. Figure 2 shows the activity of four exemplary features from session 5 chosen to illustrate tuning to force, grasp, both force and grasp independently, and an interaction between force and grasp, as evaluated with two-way Welch-ANOVA (corrected p , 0.05, Benjamini-Hochberg procedure). These features demonstrate neural modulation to forces that T8 attempted to produce using all five grasp conditions: closed pinch, open pinch, ring pinch, power grasp, and elbow extension. Extended Data Figure  2-1 shows the activity of four additional features from participant T5. TC features are labeled from 1 to 192 according to the recording electrodes from which they were extracted. Corresponding SBP features are labeled from 193 to 384.

Characterization of individual neural features
For each feature, column 1 shows neural activity that was averaged across grasp types (within force levels), resulting in trial-averaged feature traces whose differences in modulation were due to force alone. Similarly, column 2 shows neural activity averaged within individual hand grasps. Here, SBP feature 302 exhibits modulation to force only (row 1), as indicated by statistically significant go-phase differentiation in activity across multiple force levels, but not across multiple grasp levels. This forceonly tuning is what might be expected for a "high-level" coding of force that is independent of grasp type. Similarly, TC feature 190 is statistically tuned to grasp only, in that it exhibits go-phase differentiation across multiple grasps, but not across multiple forces. SBP feature 201, in which multiple forces and multiple grasps are statistically discriminable, is tuned to both force and grasp. Figure 2, column 3, displays a graphical representation of the simple main effects of the two-way Welch-ANOVA analysis, as shown by mean go-phase neural deviations from baseline feature activity during the production of each individual force level using each individual grasp type. Here, SBP features 302 and 201, which were both tuned to force independent of grasp, showed similar patterns in modulation to light, medium, and hard forces within individual grasp types. In contrast, TC feature 83 was tuned to an interaction between force and grasp; accordingly, its modulation to light, medium, and hard forces varied according which grasp type the participant used to emulate these forces. This type of interaction is what might be expected for a more "motoric" encoding of force and grasp type. If each grasp requires a different set of muscles and joints to be active, then a motoric encoding of joint or muscle motion would end up representing force differently depending on the grasp. Figure 3 summarizes the tuning properties of all 384 TC and SBP neural features in participants T8 and T5, as evaluated with robust two-way Welch-ANOVA. Specifically, Figure 3A shows the fraction of neural features tuned to force, grasp, both force and grasp, and an interaction between force and grasp. Features belonging to the former three groups (i.e., those that exhibited no interactions between force and grasp tuning) were deemed as independently tuned to force and/or grasp. As shown in row 1, the proportion of features belonging to each of these groups varied considerably across experimental sessions. However, during all sessions in both participants, a substantial proportion of features (ranging from 15.4% to 54.7% of the feature population across sessions) were  Fig. 2-1). Rows indicate average per-condition activity (PSTH) of four exemplary features tuned to force, grasp, both factors, and an interaction between force and grasp, recorded during session 5 from participant T8 (two-way Welch-ANOVA, corrected p , 0.05, Benjamini-Hochberg method). Bolded, starred p values indicate significant tuning to force (Rows 1 and 3), grasp (Rows 2 and 3), or a force-grasp interaction (Row 4). Neural activity was normalized by subtracting block-specific mean feature activity within each recording block, and then smoothed with a 100-ms Gaussian kernel to aid in visualization. Column 1 contains PSTHs averaged within individual force levels (across multiple grasps), such that observable differences between data traces are because of force alone. Similarly, column 2 shows PSTHs averaged within individual grasps (across multiple forces). Column 3 shows a graphical representation of the simple main effects as normalized mean neural deviations from baseline activity during force trials within each of the five grasps. (cp, cpinch = closed pinch; op, o-pinch = open pinch; rp, r-pinch = ring pinch, pow = power, elb = elbow extension). Mean neural deviations were computed over the go phase of each trial and subsequently averaged within each force-grasp pair. Error bars indicate 95% confidence intervals. tuned to force, independent of grasp. In other words, a substantial portion of the measured neural population represented force and grasp independently.
A smaller subset of features exhibited an interaction between force and grasp in both T8 (5.2 6 4.2%) and T5 (13.8%). Figure 3, row 2, further separates these interacting features into those that exhibited force tuning within each individual grasp type, as evaluated by one-way Welch-ANOVA (corrected p , 0.05). Here, the proportion of interacting features tuned to force appeared to depend on grasp type, particularly during sessions 2, 4, 5, 6, and 7, in a session-specific manner. In other words, within a small contingent of the neural feature population, force representation showed some dependence on intended grasp. Taken together, Figure 3 suggests that force and grasp are represented both independently and dependently within motor cortex at the level of individual neural features.

Neural population analysis and decoding Simulated force encoding models
The goal of this study was to clarify the degree to which hand grasps affect neural force representation and decoding performance, in light of conflicting evidence of grasp-independent (Chib et al., 2009;Hendrix et al., 2009;Pistohl et al., 2012;Intveld et al., 2018) versus grasp-dependent (Hepp-Reymond et al., 1999;Degenhart et al., 2011) force representation in the literature. Before visualizing population-level representation of force, we first illustrate these differing hypotheses with a toy example of expected grasp-independent versus grasp-dependent (interacting) representations of force within the neural space. Figure 4 simulates grasp-independent force encoding with an additive model (Eq. 4), and grasp-dependent force encoding with a scalar model (Eq. 5; reproduced in Fig. 4, row 1). Within the additive model, the overall neural activity x ij generated during a grasp i and force j is represented as a summation of independent force-related and grasp-related contributions. Thus, the additive model simulates independent neural force representation, in which force is represented at a high level independent of grasp. In contrast, the scalar encoding model simulates the neural activity x ij as resulting from a multiplication of the force level s j and the baseline grasping activity g i . Such an effect might be expected if force were encoded as low-level tuning to muscle activity. In this case, different force levels would result in the same pattern of muscle activity being activated to a lesser or greater degree, thus scaling the neural activity associated with that grasp, resulting in a coupling between force and grasp. Therefore, the scalar model simulates an interacting (grasp-dependent) neural force representation. Figure 4, row 2, shows simulated neural activity resulting from the additive and scalar encoding models within two-dimensional PCA space. In the independent model, force is represented in a consistent way across multiple simulated grasps, as indicated by the force axis. In contrast, within the interacting model, force representation differs according to grasp. These differences are further highlighted in Figure 4, row 3, in which dPCA was applied to the simulated neural data (over 20 simulated trials) resulting from each model. While the additive model exhibited no interaction-related neural variance, the scalar model yielded a substantial proportion of force, grasp, and interaction-related variance. Note that within these toy models, the simulated neural activity did not vary over its time course and, thus, exhibited no condition-independent (time-related) variance.
Neural population analysis Figure 5 shows neural population-level activity patterns during sessions 5 and 7 from participants T8 and T5, Figure 3. Summary of neural feature population tuning to force and grasp. Row 1, Fraction of neural features significantly tuned to force, grasp, both force, and grasp and an interaction between force and grasp in participants T8 and T5 (two-way Welch-ANOVA, corrected p , 0.05). Row 2, Fraction of neural features significantly tuned to an interaction between force and grasp, subdivided into force-tuned features within each individual grasp (c-pinch = closed pinch, o-pinch = open pinch, r-pinch = ring pinch). Note that the number of grasp types differed between sessions (see Table 1).
respectively. Here, session 5 data were shown to illustrate the neural population response to forces emulated using all five grasp conditions. Additionally, session 7 data were shown as the representative dataset from participant T5. In the first two columns of Figure 5, dPCA and traditional PCA were applied to all force-grasp conditions in both participants. In the third column, these dimensionality reduction techniques were applied solely to force trials attempted using the power grasping and elbow extension, to further quantify force representation across the entire upper limb. Population-level activity patterns for additional sessions are shown in Extended Data Figures 5-1, 5-2.
The 12 dPCs shown in Figure 5A explain the highest amount of variance within each of the four marginalizations of interest, for each participant. For example, participant T8's component #4 (row 2, column 1) is the largest force-tuned component in the dataset and explains 3.3% of the neural data's overall variance. Similarly, T8's component #2 (row 3, column 1), which captures grasp-related activity, explains 8.1% of neural variance. Horizontal black bars on each panel indicate time points at which individual dPC decoding axes predict intended forces (row 2), grasps (row 3), and forcegrasp pairs (row 4) more accurately than chance performance. In both participants, single components were able to offline-decode intended forces at abovechance levels solely during the active "go" phase of the trial, indicated by the vertical gray lines. However, grasptuned components were able to accurately predict intended grasps at nearly all time points during the trial, including the prep and stop phases. These trends were observed when dPCA was applied across all force-grasp conditions (columns 1 and 2) and across solely power and elbow trials in participant T8 (column 3). Figure 5B summarizes the variance accounted for by the entire set of dPCs extracted from each dataset. Specifically, the first row shows the cumulative variance captured by the dPCs (red), as compared with components extracted with traditional PCA (black). Here, dPCs extracted from different marginalizations were not necessarily orthogonal and accounted for less cumulative variance than traditional PCs because the axes were optimized for demixing in addition to capturing maximum variance. However, the cumulative dPC variance approached total signal variance, as indicated by the dashed horizontal lines in each panel, and were thus deemed as a faithful representation of the neural population data. . Simulated models of independent and interacting (grasp-dependent) neural representations of force. Row 1, Equations corresponding to the independent and interacting models of force representation. Here, x ij represents neural feature activity generated during a particular grasp i and force j, g i represents baseline feature activity during grasp i, f represents force-related neural feature activity, and s j is a discrete force level. Row 2, Simulated population neural activity projected into a two-dimensional PCA space. Estimated force axes within the low-dimensional spaces are shown as blue lines. Row 3, Summary of variances accounted for by the top 20 dPCs extracted from the simulated neural data from each model. Here, the variance of each individual component is separated by marginalization (force, grasp, and interaction between force and grasp). Pie charts indicate the percentage of total signal variance due to these marginalizations. Figure 5. Neural population-level activity patterns. A, Demixed principal components (dPCs) isolated from all force-grasp conditions from T8 session 5, all force-grasp conditions from T5 session 7, and power versus elbow conditions from T8 session 5 neural data. Figure 5B, second row, further subdivides the variances of individual dPCs into per-marginalization variances.
Here, most of the variance in each extracted component can be attributed to one primary marginalization, indicating that the extracted components are fairly well demixed. Pie charts indicate the percentage of total signal variance (excluding noise) from force, grasp, force/grasp interaction, and condition-independent signal components. In both participants, condition-independent components accounted for the highest amount of neural signal variance, followed by grasp, then force, then force-grasp interactions. In other words, more variance could be attributed to putative grasp representation than force representation at the level of the neural population. Additionally, force-grasp interactions only accounted for a small amount of neural variance, even when dPCA was applied solely across power grasping and elbow extension trials (column 3). Session 5 contained a larger amount of interaction variance than other sessions, possibly because of the presence of elbow extension trials that were attempted over a larger range of forces than those emulated with distal hand grasps. However, interaction variance was nonetheless smaller than force-related and grasp-related variance. Figure 5C visualizes the trial-averaged, go-phase-averaged neural activity from each dataset within two-dimensional PCA space. Within these plots, each data point represents the average neural activity corresponding to an individual force-grasp condition. In all panels, light, medium, and hard forces, represented as different shapes within PCA space, are aligned to a consistent force axis (shown in blue) across multiple grasps, and also across power grasping and elbow extension movements.
Finally, Figure 5D quantifies how well the data can be explained by the additive (grasp-independent) and scalar (grasp-dependent) encoding models presented in Equations 4, 5 and illustrated in Figure 4. Fitted model coefficients obtained via cross-validated ordinary least squares regression are indicated within session-specific tables, while R 2 values for the trained models are indicated as bar plots. Here, the additive model significantly outperformed the scalar model for all sessions (p , 0.001, one-way ANOVA, Tukey method). In agreement with this result, Figure 5B,C resemble the simulation results from the additive force encoding model ( Fig. 4; Eq. 4), which would be expected for grasp-independent force representation. However, a small amount of interactionrelated variance was also present in Figure 5B, and the force activity patterns in Figure 5C deviated to a small degree from the force axis, indicating that the additive model may not fully explain the neural activity. Therefore, the neural data were also fit to a combined model (Eq. 6), which incorporated terms from both the additive and scalar models. When fitted to neural data recorded from all force-grasp conditions (columns 1 and 2), the combined model performed similarly to the additive model (p . 0.05), likely because the scalar term within the combined model was assigned a low weighting coefficient c. However, when applied solely to the power and elbow extension trials of session 5, the combined slightly outperformed the additive model (p , 0.01), in agreement with the slightly larger force/grasp interaction-related variance present within this subset of the data.
Time-dependent decoding performance Figure 6 summarizes the degree to which intended forces and grasps could be predicted from the neural activity using the aforementioned dPCs. Here, offline force decoding accuracies were computed by using a force decoder D F , created by assembling the decoding axes of multiple force-tuned and interacting components, to classify light, medium, and hard forces over multiple sessionruns of a 100-fold, stratified, leave-group-out Monte Carlo cross-validation scheme, as described in the Methods. Similarly, grasp decoding accuracies in row three were computed using a grasp decoder D G , created by assembling the decoding axes of grasp-tuned and interacting dPCs. Figure 6, row 1, shows time-dependent force decoding results, averaged over S Â 40 session-runs in participants T8 (S = 6) and T5 (S = 1). Row 2 further subdivides the results of row 1 into force decoding accuracies achieved during individual hand grasps. Finally, row 3, shows time-dependent grasp decoding results for both participants.
Here, intended forces were decoded at levels exceeding the upper bound of chance solely during the go phase across all sessions (Extended Data Fig. 6-3), regardless of the grasp used to emulate the force. The exception to this trend occurred during elbow extension trials, in which intended forces were decoded above chance during the stop phase. In contrast, intended grasps were decoded continued dPCs were tuned to four marginalizations of interest: Condition-Independent (CI) tuning (i.e., time), Force, Grasp, and an interaction between force and grasp (FxGrasp). dPCs that account for the highest amount of variance in the per-marginalization neural activity are shown here. These variances are included in brackets next to each component number. Vertical bars indicate the start and end of the go phase. Horizontal bars indicate time points at which the decoder axes of the pictured components classified forces (row 2), grasps (row 3), or force-grasp pairs (row 4) significantly above chance. B, Summary of variances accounted for by the top 20 dPCs and PCs from each session. Here, the variance accounted for by the dPCs approaches the variance accounted for by traditional PCs. Horizontal dashed lines indicate total signal variance, excluding noise. Row 2 shows the variance of each individual component, separated by marginalization. C, Go-phase activity within a two-dimensional PCA space. Estimated force axes within the low-dimensional PCA spaces are shown as blue lines. Here, c-pinch = closed pinch, o-pinch = open pinch, r-pinch = ring pinch. D, Encoding model performances. The task-dependent components of neural feature activity were fit to the additive, scalar, and combined encoding models via cross-validated ordinary least squares regression. Tables contain the fit model coefficients for each session. Bar graphs indicate mean R 2 values for each model over 100 iterations of Monte Carlo leave-group-out cross-validation. Error bars indicate SDs across iterations. Stars indicate statistically significant differences between model R 2 values; **p , 0.01 and ***p , 0.001. above chance during all trial phases, regardless of the number of grasps from which the decoder discriminated (Extended Data Fig. 6-2), although go-phase grasp decoding accuracies tended to exceed those achieved during other trial phases. In summary, both intended forces and grasps were decoded above chance during time periods when participants intended to produce these forces and grasps, and in some cases, during preparatory and stop periods. Session-averaged, time dependent decoding accuracies for individual force levels and grasp types are displayed in Extended Data Figure 6-1. Figure 7 summarizes go-phase force and grasp decoding accuracies as confusion matrices. Here, time-dependent classification accuracies for each force level and each grasp type were averaged over go-phase time windows (see Fig. 6) that commenced when overall classification performance exceeded 90% of their maximum, and ended with the end of the go phase. This time period was selected to exclude the rise time in classification accuracy at the beginning of the go phase, so that the resulting mean trial accuracies reflected stable values. The mean trial accuracies were then averaged over all session-runs in each participant to yield confusion matrices of true versus predicted forces and grasps. Figure 7B further subdivides overall three-force classification accuracies into force classification accuracies achieved during each individual grasp type (columns) in both participants (rows). The confusion matrices in Figure 7 represent cumulative data across multiple sessions in participant T8, and one session in participant T5. Extended Data Figures 7-1, 7-2, 7-3 statistically compare decoding accuracies between individual force levels and grasp types within each individual session.

Go-phase decoding performance
In Figures 6A, 7A and Extended Data Figure 6-1, overall, three-force classification accuracies exceeded the upper limit of chance in both participants. However, the decoding accuracies of individual force levels were statistically different. For almost all sessions, hard forces were classified more accurately than light forces (with the exception of session 4, during which light and hard force classification accuracy was statistically similar), and both light and hard forces were always classified more accurately than medium forces. More specifically, hard and light forces were decoded above chance across all sessions, while medium force classification accuracies often failed to exceed chance in both participants. In contrast, both overall and individual grasp decoding accuracies always exceeded the upper limit of chance. According to Figure 7A and Extended Data Figure 7-1B, certain grasps were decoded more accurately than others. Specifically, in participant T8, the power and ring pincer grasps were often classified more accurately than the open and closed pincer grasps across multiple sessions (corrected p ( 0.05, one-way ANOVA). Elbow  Fig. 6) were averaged over go-phase time windows that commenced when performance exceeded 90% of maximum and ended with the end of the go phase. These yielded mean trial accuracies, which were then averaged over all session-runs in each participant. Overall force and grasp classification accuracies are indicated above each confusion matrix. SDs across multiple session-runs are indicated next to mean accuracies (cp = closed pinch, op = open pinch, rp = ring pinch, pow = power, elb = elbow extension). Statistical comparisons between the achieved classification accuracies are shown in Extended Data Figure 7-1. B, Confusion matrices, now separated by the grasps (c-pinch = closed pinch, opinch = open pinch, r-pinch = ring pinch, power, elbow) that participants T8 (row 1) and T5 (row 2) used to attempt producing forces. Statistical comparisons between the achieved force accuracies are shown in Extended Data Figures 7-2 , 7-3. extension, which required the participants to attempt force production in the upper limb in addition to the hand, was classified more accurately than any of the grasping forces during session 5 (corrected p ( 0.05). In participant T5, grasp classification accuracies, in order from greatest to least, were ring pincer . open pincer . power . closed pincer. Regardless, grasp decoding performance always exceeded force decoding performance in both participants, as seen in Figures 6, 7.
In Figure 7 and Extended Data Figure 7-3, overall and individual force classification accuracies varied depending on the hand grasps used to attempt these forces. Specifically, classification accuracies for forces attempted with different grasps were, with few exceptions, statistically different (corrected p ( 0.05, one-way ANOVA). For example, in Figure 7B and Extended Data Figure 7-3, hard forces attempted using the open pincer grasp were always classified more accurately than hard forces attempted using the ring pincer grasp in both participants. In other words, grasp type affected how accurately forces were decoded.
Finally, Figure 8 summarizes how well force decoders trained on one set of grasps generalized to novel grasp types in T8 session 5 (row 1) and T5 session 7 (row 2). A force decoder was used to discriminate forces among a set of grasps used for training ("left-in," gray bars) or a leave-out novel grasp (white bars). Here, the force decoding performance between the leave-in and leave-out grasps was significantly different in seven out of nine comparisons, suggesting that grasp affects how well forces are decoded from neural activity. However, for all sets of grasps, force decoding performance always exceeded chance. This was even true when, during T8 session 5, the force decoder was trained on four hand grasps and evaluated on elbow extension data. This is consistent with the previous population-level analyses that show that components of force representation in motor cortex are conserved across grasps and even arm movements.

Discussion
The current study sought to determine how human motor cortex encodes hand grasps and discrete forces, how much these representations interacted, and how well forces and grasps could be decoded. Three major findings emerged from this work. First, force information was present in, and could be decoded from, intracortical neural activity in a consistent way across multiple hand grasps. This suggests that force is, to some extent, represented at a high level in individuals with tetraplegia, independent of motion and grasp. However, as a second finding, grasp affected force representation and classification accuracy, suggesting a simultaneous, low-level, motoric representation of force in individuals with tetraplegia. Finally, hand grasps were classified more accurately and explained more neural variance than forces. These three findings and their implications for future online force decoding efforts are discussed here.

Force and grasp representation in motor cortex
Force information persists across multiple hand grasps in individuals with tetraplegia.
Overall force representation. Force was represented in a consistent way across multiple hand grasps within the neural activity. In particular, a substantial contingent of neural features was tuned to force independent of grasp ( Fig. 3), force-tuned components explained more population-level variance than components tuned to force-grasp interactions (Fig. 5), and intended forces were accurately predicted from population-level activity across multiple grasps (Figs. 6-8). The study results suggest that in individuals with tetraplegia, to a large extent, force is represented at a high level within motor cortex, distinct from grasp, in accordance with the grasp-independent force encoding model described by Equation 4 (Figs. 4, 5D). This conclusion agrees with previous motor control studies (Mason et al., 2004;Chib et al., 2009;Casadio et al., 2015), which suggest that at the macroscopic level, force Figure 8. Go-phase force classification accuracy for novel (test) grasps. Within each session (rows), dPCA force decoders were trained on neural data generated during all grasps, excluding a single leave-out grasp type (columns). The force decoder was then evaluated over the set of training grasps (gray bars), as well as the novel leave-out grasp type (white bars). Stars indicate statistically significant differences in performance between training and novel grasps; **p , 0.01, ***p , 0.001. Error bars indicate the 95% confidence intervals. The horizontal dotted line in each panel indicates upper bound of the empirical chance distribution for force classification. Here, c-pinch = closed pinch, o-pinch = open pinch, r-pinch = ring pinch. and motion may be represented independently. In particular, Chib and colleagues showed that descending commands pertaining to force and motion could be independently disrupted via transcranial magnetic stimulation (TMS), and that these commands obeyed simple linear superposition laws when force and motion tasks were combined.
Furthermore, intracortical non-human primate studies (Mason et al., 2006;Hendrix et al., 2009;Intveld et al., 2018) suggest that forces are encoded largely independently of the grasps used to produce them. However, in these studies and within the present work, hand grasps likely recruited overlapping sets of muscle activations. Thus, the relatively low degree of interactions observed here and in the literature could actually be because of overlapping muscle activations rather than truly grasp-independent force representation. For this reason, participant T8 emulated forces using elbow extension in addition to the other hand grasps during session 5. The elbow extension task, which recruited both proximal and distal muscle activations, was chosen to overlap less with the other hand grasps, which recruited distal muscle activations only. In Figure 5, column 3, dPCA was implemented solely on force trials emulated using elbow extension and power grasping. The resulting dPCA composition yielded a slightly larger interaction variance (4%) that was nonetheless smaller than variance due to force (;12%) or grasp (;35%). Furthermore, discrete force data, when represented within two-dimensional PCA space, aligned closely with a force axis that was conserved over both power grasping and elbow extension movements, providing further evidence that force may be encoded independently of movements and grasps.
Representation of discrete forces. While overall force accuracies exceeded chance performance (Fig. 6), hard and light forces were classified more accurately than medium forces across all hand grasps, sessions, and participants. Medium forces often failed to exceed chance classification performance ( Fig. 7A; Extended Data Figs. 6-1B, 6-3). Notably, classification performance depended on participants' ability to kinesthetically attempt various force levels and grasps without feedback, despite having tetraplegia for several years before study enrollment. Anecdotally, participant T8 reported that light and hard forces were easier to attempt than medium forces, because they fell at the extremes of the force spectrum and could thus be reproduced consistently. Although his confidence with reproducing all forces improved with training, it is conceivable that without sensory feedback, medium forces were simply more difficult to emulate, and thus yielded neural activity patterns that were less consistent and more difficult to discriminate.
Additionally, prior studies suggest that neural activity increases monotonically with increasing force magnitude (Evarts, 1969;Thach, 1978;Cheney and Fetz, 1980;Wannier et al., 1991;Ashe, 1997;Cramer et al., 2002). Therefore, by virtue of being intermediate to light and hard forces, medium forces may be represented intermediate to light and hard forces in the neural space, and may thus be more easily confused with forces at the extremes of the range evaluated (Murphy et al., 2016;Downey et al., 2018). To this point, population-level activity during medium and light forces exhibited similarities ( Fig. 5; Extended Data Fig. 5-1); accordingly, medium forces were most often confused with light forces during offline classification (Fig. 7).

Hand posture affects force representation and force classification accuracy
Single-feature versus population interactions between force and grasp. As previously stated, force information was neurally represented, and could be decoded, across multiple hand grasps (Figs. 3, 5-8). However, hand grasp also influenced how force information was represented within (Fig. 5) and decoded from ( Fig. 7; Extended Data Fig. 7-3) motor cortex. Furthermore, despite small force-grasp interaction population-level variance ( Fig. 5B; Extended Data Fig. 5-1B), as many as 12.0% and 13.8% of neural features exhibited tuning to these interaction effects in participants T8 and T5, respectively (Fig. 3), providing further evidence that the force and grasp representation are not entirely independent.
When considering the relatively large number of interacting features and the small population-level interaction variance, one might initially conclude that a discrepancy exists between feature-level and population-level representation of forces and grasps. However, the amount of variance explained by a parameter of interest may not always correspond directly to the percentage of features tuned to this parameter. Here, the interaction effects within individual features likely reached statistical significance with small effect size. In other words, while real interaction effects were present within the feature data (Fig. 3), the overall effect was small, as exhibited within the population activity (Fig. 5). From this perspective, the seemingly incongruous feature-level and population-level results actually complement one another and inform our understanding of how forces are represented in motor cortex in individuals with tetraplegia.
Force and grasp have both abstract (independent) and motoric (interacting) representations in cortex. Thus far, studies of force versus grasp representation have largely fallen into two opposing groups. The first proposes that motor parameters are represented independently (Carmena et al., 2003;Mason et al., 2006;Hendrix et al., 2009;Intveld et al., 2018). Such representation implies that the motor cortex encodes an action separately from its intensity, then combines these two events downstream to compute the EMG patterns necessary to realize actions in physical space.
In contrast, the second group suggests that force, grasp, and other motor parameters interact (Hepp-Reymond et al., 1999;Degenhart et al., 2011). They propose that motor parameters cannot be fully de-coupled (Kalaska, 2009;Branco et al., 2019) and that it may be more effective to use the entire motor output to develop a comprehensive mechanical model, rather than trying to extract single parameters such as force and grasp .
The current study presents evidence supporting both independent and interacting representations of force and grasp in individuals with tetraplegia. These seemingly contradictory results actually agree with a previous nonhuman primate study that recorded from motor areas during six combinations of forces and grasps (Intveld et al., 2018). Intveld and colleagues found that, while forcegrasp interactions explained only 0-3% of population variance, roughly 10-20% of recorded neurons exhibited such interactions, which is highly consistent with the present study results (Figs. 3, 5). Thus, in individuals with tetraplegia, the neural space could consist of two contingents: one that encodes force at a high level independent of grasp and motion, and another that encodes force as low-level tuning to muscle activity, resulting in interactions between force and grasp. The second contingent, however small, significantly impacts how accurately forces and grasps are decoded ( Fig. 7B; Extended Data Fig. 7-3) and should not be discounted.
Hand grasp is represented to a greater degree than force at the level of the neural population Go-phase grasp representation. In the present datasets, grasps were decoded more accurately (Figs. 6,7; and explained more signal variance ( Fig. 5B; Extended Data Fig. 5-1B) than forces. This suggests that within the sampled region of motor cortex, grasp is represented to a greater degree than force, which agrees with prior literature Milekovic et al., 2015;Intveld et al., 2018).
In the current work, force may be represented to a lesser degree than grasp for several reasons. First, force information may have stronger representation in caudal M1, particularly on the banks of the central sulcus (Kalaska and Hyde, 1985;Sergio et al., 2005; or within the depth of the sulcus (Rathelot and Strick, 2009), which cannot be accessed using planar microelectrode arrays. Second, force-tuned neurons in motor cortex respond more to the direction of applied force than its magnitude (Kalaska and Hyde, 1985;Kalaska et al., 1989;Taira et al., 1996). Finally, intracortical non-human primate studies (Georgopoulos et al., 1983(Georgopoulos et al., , 1992 and human fMRI studies (Branco et al., 2019) suggest that motor cortical neurons respond more to the dynamics of force than to static force tasks. The present work, which recorded from rostral motor cortex while study participants emulated static, non-directional forces, may therefore have detected weaker force representation than would have been possible from more caudallyplaced recording arrays during a dynamic, functional force task.
Additionally, both study participants were deafferented and received no sensory feedback regarding the forces and grasps they attempted. In individuals with tetraplegia, discrepancies may exist between the representation of kinematic parameters such as grasp, which remain relatively intact because of their reliance on visual feedback, and kinetic parameters such as force (Rastogi et al., 2020). Specifically, since force representation relies heavily on somatosensory feedback (Tan et al., 2014;Tabot et al., 2015;Schiefer et al., 2018), whose neural pathways are altered during tetraplegia (Solstrand Dahlberg et al., 2018), the current study may have yielded weaker forcerelated representation than if this feedback were present. Therefore, further investigations of force representation are needed in individuals with tetraplegia during naturalistic, dynamic tasks that incorporate sensory feedback, either from intact sensation or from intracortical microstimulation (Flesher et al., 2016), to determine the full extent of motor cortical force representation and to maximize decoding performance.
Grasp representation during prep and stop phases. Unlike forces, which were represented primarily during the go phase of the trial, grasps were represented throughout the entire task (Figs. 5,6), in agreement with previous literature (Milekovic et al., 2015). However, this ubiquitous grasp representation may be partially explained by the behavioral task. Research sessions consisted of multiple data collection blocks, each of which was assigned to a particular hand grasp, and cycled through three attempted force levels within each block (Fig. 1B). Thus, while attempted force varied from trial to trial, attempted hand grasps were constant over each block and known by participants in advance. When individuals have prior knowledge of one task parameter, but not another other, information about the known parameter can appear within the baseline activity (Vargas-Irwin et al., 2018). Therefore, grasp-related information may have been represented within the neural space during non-active phases of the trial, simply by virtue of being known in advance.
Additionally, the placement of the recording arrays could have influenced grasp representation in this study. In each participant, two microelectrode arrays were placed within the hand knob of motor cortex (Yousry et al., 1997). These arrays may have recorded from "visuomotor neurons," which modulate both to grasp execution and to the presence of graspable objects before active grasp (Carpaneto et al., 2011), or from neurons that are involved with motor planning of grasp (Schaffelhofer et al., 2015). These neurons have typically been attributed to area F5, a homolog of premotor cortex in non-human primates. Recent literature indicates that human precentral gyrus is actually part of the premotor cortex (Willett et al., 2020). Thus, the arrays in this study likely recorded from premotor neurons, which modulate to grasp during both visuomotor planning and grasp execution, as was observed here.

Implications for force decoding Hand grasp affects force decoding performance
Our decoding results demonstrate that, in individuals with tetraplegia, forces can be decoded offline from neural activity across multiple hand grasps (Figs. 6-8). These results agree with the largely independent force and grasp representation within single features (Fig. 3) and the neural population (Fig. 5). From a functional standpoint, this supports the feasibility of incorporating force control into real-time iBCI applications. On the other hand, grasp affects how accurately discrete forces are predicted from neural data ( Fig. 7B; Extended Data Fig. 7-3). Therefore, future robust force decoders may need to account for additional motor parameters, including hand grasp, to maximize performance.