A versatile robotic platform for the design of natural, three-dimensional reaching and grasping tasks in monkeys

The IIWA robotic arm features a large workspace (figure S1) allowing ranges of motion that are compatible with both human and monkey reaches. Additionally, the robotic arm is able to actively lift up to 7 Kg of weight, which makes it robust to manipulation by monkeys.

We developed a software package that implements a real-time closed-loop control (figure 2(A), https://doi.org/10.5281/zenodo.3234138) configured as a finite state machine. This allows fast configuration of tasks that proceed through several phases, where each phase requires a different behaviour of the robot.

Zoom In Zoom Out Reset image size

In our specific example application, at the beginning of an experimental session the robot lies in home position (figure 2(A), FS1), in which all joint coordinates are equal to 0°, resulting in a straight and vertical robotic arm configuration (supplementary figure S1(C) (stacks.iop.org/JNE/17/016004/mmedia)). The robotic arm maintains its position until the user triggers the start of a new trial by pressing a remote button. This brings the robot in a position control phase (figure 2(A), FS2) during which it moves the end effector to a determined position in space.

Once the sum of the errors in the positioning of all the joints decreases to a value of 0.01 radians, the robot is considered in target position and the control mode switches to impedance mode (figure 2(A), FS3). In this state the robot behaves as a mass-spring-damper system, trying to keep the target position while opposing a dynamic resistance to applied forces. In our specific examples we facilitated smooth movements of the end effector along the x-axis (towards the monkey), by imposing higher stiffnesses and damping parameters along the other orthogonal directions.

To interact with our platform we trained three monkeys to grasp a flexible instrumented object connected to the end effector, and pull it until a user-defined-position along the x-axis is crossed. Upon success, the robot quickly updates the next desired target position and moves towards it (figure 2(A), FS4, FS5). If the animal fails to pull the object across the threshold or the trial timeouts, the robot returns to target position. After three failed trials the robot switches to another target.

We designed an air-pressure-based gauge sensor that can measure the applied grip force independently from an object shape in the interval  −0.7 kPa–  40 kPa. We 3D-printed hollow objects of different shape and size (figure 2(B)) and connected to a steel enclosure hosting the sensor. In this configuration the grip strength is proportional to the air flux produced by compression of the object. We designed a programmable electronic circuit (figure 2(B)) that digitizes pressure measurements while power and data transfer are provided through connection ports located at the robot end effector thereby limiting external wires at the subject-robot interface.

The pressure sensor circuit converted air pressure into a voltage measurement using strain-gages (1620 measurement SPECIALTIES^™) directly connected to a front-end amplifier (FEA—LMP90100) of a full-Wheatstone bridge circuit. The footprint of the electronic circuit was maximally reduced in order to integrate it to the robot (diameter 3.2 cm). The embedded circuit converts the strain-gages signal into an analog output signal bounded between 0 and 4.1 V. The electronic assembly comprises a microcontroller (pic24FV08KM101) for the control of the FEA and the digital-to-analog converter (DAC8551) via a serial peripheral interface (SPI) protocol. Signals were amplified, digitized at 16 bits, and converted in analog signals after noise reduction. A 5 V regulator (EG113NA-5) adjusted the input power supply from 6 V to 9 V and a second regulator (LM4140ACM-4.1) outputted a precise reference voltage of 4.1 V used for signal conversions. Three LEDs indicated the state of the system and three input switches were used to impose gains ranging from 1 to 64. The in-circuit serial programming (ISCP) connector served to re-program the microcontroller. The system was also equipped from an ON/OFF button, a reset button and a power supply connector. The assembly was mechanically and electrically coupled to the flange of the robotic arm.

The objects were printed in a thermoplastic elastomer and were covered with two layers of silicon (DOWSIL^™ and Sil-Poxy^®) to ensure sealing. All the sensitive elements were characterized using a computer-controlled compression system (Zwick/Roell 1KN D0728165) to derive the internal pressure of each object as a function of the applied force.

All the surgical procedures were performed under full anaesthesia induced with midazolam (0.1 mg kg⁻¹) and ketamine (10 mg kg⁻¹, intramuscular injection) and maintained under continuous intravenous infusion of propofol (5 ml/kg/h) and fentanyl (0.2–1.7 ml/kg/h) using standard aseptic techniques. A certified neurosurgeon (Dr Jocelyne Bloch, CHUV, Lausanne, Switzerland) supervised all the surgical procedures. We implanted two 64-channel microelectrode arrays in Mk-Jo and Mk-Cs and two 48-channel microelectrode arrays in Mk-Ol (Blackrock Microsystems, 400 µm pitch and electrodes tip lengths 1.5 mm and 1 mm for M1 and S1, respectively). Mk-Cs was implanted in M1 and S1 area of the arm, and Mk-Jo and Mk-Ol were implanted in M1 and S1 area of the hand. Functional areas were identified with electrical stimulation delivered as biphasic pulses on the cortex surface at 3 mA and 300 Hz. A 20 mm diameter craniotomy was performed in order to span the brain areas of interest and the dura was incised. Implantation of the arrays was achieved using a pneumatic compressor system (Impactor System, Blackrock Microsystems). The pedestal was fixated to a compliant titanium mesh modelled to fit the skull shape for Mk-Jo and Mk-Cs. The pedestal was fixed directly to the skull in Mk-Ol. Surgical and post-operative care procedure are developed in details in Capogrosso et al (2018). Data presented in this paper were collected 3, 8 and 9 weeks post-implantation for Mk-Jo, Mk-Ol and Mk-Cs, respectively.

We built a custom primate chair (supplementary figure S1) in which only the neck of the animal was fixed with a metallic collar allowing wide range of voluntary arm and hand movements in the three-dimensional space. Two Plexiglas plates were used to restrain the access to the head and a third plate placed around the stomach of the animal served for placing a resting bar and elastic bands to immobilize the right hand. Monkeys were trained to maintain a resting position between trials and place the left hand on the resting bar a few centimetres in the front, at chest level. A trial started when the robot presented a graspable target in front of the animal, at a distance of approximatively 20 cm. As soon as the 'go' cue was played (1 s duration sound), monkeys reached for the target, grasped it, and pulled it towards themselves (x-direction). Once the robot end effector crossed a pre-set virtual spatial threshold (8 cm), a clicker sound played by the experimenter indicated the success of the trial encouraging the animal to release the grip, return to the resting bar and get a food reward. The reward was delivered manually by the experimenter for Mk-Jo and Mk-Cs and automatically for Mk-Ol. The robot returned to its vertical home position at the end of each trial. For Mk-Cs and Mk-Jo, cueing signals during the task were implemented as follow: a go cue was implemented as a high pitch sound coupled to a green light that was played when the robot was in position and reached the impedance control mode; a success cue consisted in a low pitch sound when the monkey crossed the spatial threshold by pulling on the robot. Sounds were played from the application controlling the robot while the light cue was triggered through a synchronisation module (Arduino Due, Arduino, Italy). For Mk-Ol, the visual cues were delivered on a computer screen. A start cue was displayed as a rectangular box on top of the screen. After 500 ms, the start box disappeared and a cursor was displayed on top of the screen with a rectangular target at the bottom. The cursor moved down vertically towards the target, proportionally to the displacement of the robot end effector towards the pre-set spatial threshold. Auditory cues were delivered automatically and consisted in a high pitch start cue, a medium pitch success cue and a low pitch reward cue. Liquid food was delivered together with the reward sound using a pump (masterflex consoledrive) triggered externally via the control software. The screen and the reward system were controlled using an in-house application (Matlab Mathworks^®) coupled to a synchronization board (National Instruments, US). All cues were identically delivered during training and experimental sessions.

In Task 1, Mk-Jo reached for, grasp and pull objects of different sizes and shapes all presented at the same location. A cylinder (length 8 cm; diameter 1.5 cm), a small sphere (diameter 1.5 cm) and a large sphere (diameter 3 cm) were used to induce cylindrical, small and wide spherical grips, respectively. CAD design of these objects are available at (https://doi.org/ 10.5281/zenodo.3234138).

Each of these targets were presented with various levels of resistance applied by the robot impedance controller (joint stiffness 200, 400 and 600 N m⁻¹ corresponding to 'low', 'medium' and 'high' stiffness levels). In addition, a small triangular pinch-like object (base 2 cm; height 1.5 cm) was used to prompt a lateral precision grip and was presented with the lowest level of resistance. Each session consisted in 12 to 20 trials per object per level of resistance. Two to four objects were presented during a single session.

In Task 2, Mk-Cs performed 3D centered-out task where the small spherical object was alternatively and randomly presented in three horizontal positions in the sagittal plane, center, left and right. The object was placed at  −40, 0 and 60 mm along the y -axis for the left, center and right position, respectively. The z distance was fixed for all conditions at approximately 180 mm above the animal seating height. Each session consisted in approximately 25 trials per position. Mk-Ol performed both Task 1 and Task 2. In Task 2, the left and right targets were placed approximately 3 cm higher than the central target along the z-coordinate. Each session consisted in at least 80 or 70 trials per condition, for position or object, respectively.

Three-dimensional spatial coordinates of arm and hand joints during upper limb movements were acquired using a 14-camera motion tracking system (figure 1, Vicon Motion Systems, Oxford, UK) at a 100 Hz-frame rate. The video system tracked the Cartesian position of up to 15 infrared reflective markers (6–9 mm in diameter each, Vicon Motion Systems, Oxford, UK). For each monkey, one marker was placed directly below the shoulder, three on the elbow (proximal, medial and distal) and two were placed on the wrist (lateral and medial) using elastic bands. For Mk-Jo and Mk-Ol, nine additional markers were positioned on the back of a customized viscose glove, on the metacarpal (MCP), proximal (PIP) and distal phalanxes (DIP) joints of the thumb, index and little finger (supplementary figure S2(A)). A model of each subject's marker placement was calibrated in Vicon's Nexus software.

The interaction force, measured as the force applied at the robotic joints, was sampled at 500 Hz and synchronized with triggers marking the beginning (go cue) and the end (spatial threshold crossing) of each trial together with the spatial position of the end effector. Values were also streamed in the BlackRock system using the single-board microcontroller. The pressure sensor signal was sampled at 1000 Hz with the VICON and Blackrock systems.

Neural signals were acquired with a Neural Signal Processor (Blackrock Microsystems, USA) using the Cereplex-E headstage with a sampling frequency of 30 kHz.

We evaluated the safety and usability of the robotic platform accuracy of positioning, stability of position over time and repeatability of behaviour across sessions (see section results).

We continuously recorded the three-dimensional position of the robot arm end effector, together with the three-dimensional force exerted on the object. We computed a positioning error, a time for the robot arm to converge on the target position (target positioning time) and a drop error (figure 3(A)). The precision error was computed as the distance between the targeted end effector position and the reached endpoint position, therefore illustrating the reproducibility of the robotic arm across trials. Since joint stiffness values influence the speed of the motion and the speed of convergence toward the commanded position we evaluated the variation of position error for several joint stiffness values ranging from 200 to 800 N m⁻¹. In addition, we enclosed the variation of the time needed for convergence toward the commanded position (time for target positioning). The drop error was computed as the three-dimensional drift of the end effector from the target position during the first 500 ms after placement, illustrating the stability of the robotic arm in holding the object in the target position. Subsequently, we verified the compliance of the robot movement upon interaction with a monkey previously trained to perform Task 2. We recorded three-dimensional coordinates of the joints and then computed a trajectory smoothness index (Hans-Leo Teulings et al 1997) computed as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle S=~\sqrt{\frac{{{T}^{5}}}{{{L}^{2}}}\sum\limits_{i=1}^{N}{{{\left\Vert {{\left( \frac{{{\partial }^{3}}\boldsymbol{p}}{\partial {{t}^{3}}} \right)}_{i}} \right\Vert}^{2}}dt}}\nonumber \end{align}$$

where p  is the end effector position, N is the number of samples considered for the measure, t is the time, T and L represents the total duration of the trajectory in seconds and the length of the trajectory in m, respectively. We computed a separated smoothness index over the reach and the pull.

Zoom In Zoom Out Reset image size

We analyzed the movement velocity by computing the maximum wrist speed during the reaching and pulling phase separately. We derived the wrist speed from the wrist kinematics of Mk-Jo using the cylinder for three levels of robotic joint resistances (200 N m⁻¹, 400 N m⁻¹ and 600 N m⁻¹).

All the grip pressure sensor-coupled objects were characterized using a computer-controlled compression system (Zwick/Roell 1KN D0728165) to derive the internal pressure of each object as a function of the applied force in Netwon. Calibration curves were acquired by applying dynamically load ranging from 0 to 10 N on each object.

Post processing of motion capture data was performed to ensure that all joint markers were labelled correctly. We converted the three-dimensional marker position data to rotational degrees of freedom. For Mk-Jo and Mk-Ol, we computed 12 joint angles (see figure S2): index finger PIP flexion/extension, pinkie finger PIP flexion/extension, thumb finger PIP flexion/extension, index finger MCP flexion/extension, pinkie finger MCP flexion/extension, thumb finger MCP flexion/extension, index finger abduction/adduction, pinkie finger abduction/adduction, thumb finger abduction/adduction, thumb opposition, wrist flexion, wrist ulnar deviation. For Mk-Cs, we computed 3 joint angles: shoulder adduction, elbow flexion/extension angle and wrist pronation/supination. Interaction force measurements were synchronized post-hoc to kinematics signals based on triggers marking the start (go cue) and end (threshold crossing) of each trial. For further analysis of arm and hand dynamics during natural reaching and grasping, we averaged the mesurements across each condition. Kinematics, force and grip pressure measurements were low-pass filtered at 10 Hz.

Principal component analysis (PCA) was performed on joint angle kinematics and kinetic variables to identify which dynamic features accounted for most of the variance in the data. We reconstructed the Y  ×  N_df matrix where N_df is the number of joint angles and Y  =  [Y₁, Y₂,..., Y_N] is the concatenated vector of all time points t  =  1,..., N. Average kinematic, force and pressure variables were normalized and mean-centered before singular values decomposition. The dynamic data projected along the first three eigenvectors in the PC space were averaged over each condition and smoothed using a 5-samples moving average filter before plotting.

We used PCA, as described in Santello et al (1998) to identify kinematic synergies as orthogonal axes of maximal correlated variance in the 3D joint coordinates of the hand. Briefly, the x, y and z joints Cartesian coordinates were first normalized and mean-centered. The principal components (PCs) were then computed from the eigenvalues and eigenvectors of the matrix of the covariance coefficients between each of the joint coordinates waveforms. As each of the 13 markers was represented by three coordinates, the PCA was computed over 39 waveforms. The first PC accounted for at least 70% of the variance. We represented the hand posture corresponding to the first synergy by plotting the 39 coefficients corresponding to the first PC in 3D (Mason et al 2001) and overlaid a sketch of a monkey hand for the sake of representation. The positions of the ring finger and middle finger were not measured in our experiment but were represented in a realistic inferred position with respect to the thumb index and pinkie to help visualize the hand posture in supplementary figure S2(D).

We acquired the multiunit spiking activity from each channel of intracortical neural recordings (2  ×  64 channels for Mk-Cs and Mk-Jo, 2  ×  48 channels for Mk-Ol) using the Neural Signal Processor (Capogrosso et al 2016). Specifically, a spiking event was defined on each channel if the band-pass filtered signal (250 Hz–5 kHz) exceeded 3.0–3.5 times its root-mean-square value calculated over a period of 5 s. Artefacts removal was achieved by eliminating all the spikes occurring within a time window of 0.5 ms after a spike event in at least 30 channels. We computed the firing rate of each channel as the number of spikes detected over non-overlapping bins of 10 ms.

In order to perform a neural population analysis, we identified a neural manifold (Gallego et al 2017) by conducting PCA on multi-unit firing rates ranging from movement onset to the end of the pulling phase. Movement onset happened in average 660 ms, 910 ms and 520 ms before the grasp for Mk-Jo, Mk-Cs and Mk-Ol, respectively. The pulling phase lasted at least 350 ms, 740 ms and 370 ms for Mk Jo, Mk-Cs and Mk-Ol, respectively. Events marking movement onset, grasp onset and release of the object were identified using video recordings. Neural features for M1 and S1 consisted in the firing rates computed over all intracortical neural channels of each array (64 channels for Mk-Cs and Mk-Jo, 48 channels for Mk-Ol). PCs were computed using at least 15 trials per object and 25 trials per position in Mk-Jo and Mk-Cs, respectively. In Mk-Ol, principal components during the objects-grasping task were computed using 70 trials per object (supplementary figure S3).

We manually spike-sorted the recordings from each electrode implanted in M1 and S1. We computed the average firing rate across the entire reach for each neuron to compare the two brain areas (figure 7(a)). Additionally, we inspected the tuning for each parameter by plotting the firing rate in each bin against the hand velocity and the magnitude of pulling force recorded at each time (figure 7(b)).

We then constructed encoding models to predict the spiking of each neuron using Generalized Linear Models (GLMs), adapting an analysis previously described by Lawlor et al (2018) and Perich et al (2018). First, we counted the number of spikes in non-overlapping 10 ms bins. In brief, GLMs generalize the idea of multilinear regression for non-Gaussian statistical distributions using a nonlinear link function. The neurons were assumed to have Poisson statistics, thus we used an exponential link function (Lawlor et al 2018). As inputs to the GLMs we provided: (1) full-limb kinematics, including the three-dimensional velocities and accelerations of the shoulder, elbow, and wrist; (2) three-dimensional pulling force applied to the robot by the monkey; (3) the time of the grab event, to capture either motor commands related to hand shaping or sensory feedback from the object contact. For the third input, the grab event was convolved with three raised cosine basis functions (Pillow et al 2008) spaced evenly up to 300 ms in the past or future for M1 and S1, respectively.

We quantified the performance of the GLMs using a pseudo-R² metric, which generalizes the notion of variance explained for the Poisson statistics of the model (Lawlor et al 2018, Perich et al 2018). This metric compares the log-likelihood of the tested model fit against a simpler model.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle p{{R}^{2}}=1-\frac{{\rm log}L\left( n \right)-{\rm log}L\left( \widehat{\lambda } \right)}{{\rm log}L\left( n \right)-{\rm log}L\left( {\bar{n}} \right)}.\nonumber \end{align}$$

Typically, the simpler model $\bar{n}$ is a mean-fit to the data. However, this formulation also allows us to test the relative contributions of different parameters. We compared the full model with all three types of inputs described above to a reduced model $\bar{n}$ where just one of the inputs is omitted. This relative pseudo-R² metric provides insights into how well specific parameters, such as object contact, helps to explain neural activity. We assessed significance for all of these using a Monte Carlo simulation resampling across the available trials. A model fit or parameter was assumed to be significant if the 95% confidence intervals on the parameter fit were greater than zero.

We implemented an approach based on a multiclass regularized linear discriminant analysis algorithm (mrLDA) to detect moments of movement onset and object grasp from the continuous neural recordings from either the motor or somatosensory cortices (Milekovic et al 2013a, Capogrosso et al 2016, Milekovic et al 2018). In brief, we synchronized the multiunit spike activity with the movement onset and object grasp events. The performance of decoders was evaluated using five-fold cross-validation (Hastie et al 2009). The movement onset and object grasp events, mo and og, were used to derive the respective classes of neural features, C_mo and C_og:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{C}_{mo}}=\left\{{{{\bar{a}}}_{mo}}\left| {{{\bar{a}}}_{mo}}\left( i \right) \right.=\left[ \begin{array}{@{}c@{}} {{x}_{1}}\left( {\rm mo}\left( i \right)-\Delta {{t}_{mo}} \right) \nonumber \\ {{x}_{1}}\left( {\rm mo}\left( i \right)-\Delta {{t}_{mo}}-\Delta t \right) \nonumber \\ \vdots \nonumber \\ {{x}_{1}}\left( {\rm mo}\left( i \right)-\Delta {{t}_{mo}}-\left( {{N}_{TP}}-1 \right)\cdot \Delta t \right) \nonumber \\ {{x}_{2}}\left( {\rm mo}\left( i \right)-\Delta {{t}_{mo}} \right) \nonumber \\ \vdots \nonumber \\ {{x}_{{{N}_{CH}}}}\left( {\rm mo}\left( i \right)-\Delta {{t}_{mo}}-\left( {{N}_{TP}}-1 \right)\cdot \Delta t \right) \nonumber \\ \end{array} \right] \right\}{{C}_{og}}=\left\{{{{\bar{a}}}_{og}}\left| {{{\bar{a}}}_{og}}\left( i \right) \right.=\left[ \begin{array}{@{}c@{}} {{x}_{1}}\left( {\rm og}\left( i \right)-\Delta {{t}_{og}} \right) \nonumber \\ {{x}_{1}}\left( {\rm og}\left( i \right)-\Delta {{t}_{og}}-\Delta t \right) \nonumber \\ \vdots \nonumber \\ {{x}_{1}}\left( {\rm og}\left( i \right)-\Delta {{t}_{og}}-\left( {{N}_{TP}}-1 \right)\cdot \Delta t \right) \nonumber \\ {{x}_{2}}\left( {\rm og}\left( i \right)-\Delta {{t}_{og}} \right) \nonumber \\ \vdots \nonumber \\ {{x}_{{{N}_{CH}}}}\left( {\rm og}\left( i \right)-\Delta {{t}_{og}}-\left( {{N}_{TP}}-1 \right)\cdot \Delta t \right) \nonumber \\ \end{array} \right] \right\}\nonumber \end{align}$$

where ${{\bar{a}}_{mo}}$ and ${{\bar{a}}_{og}}$ are feature vector members of classes C_mo and C_og, respectively; N_TP is the number of multiunit spike rate measurements taken from the same neural channel; Δt is the temporal difference between the two consecutive spike rate measurements taken from the same neural channel; N_CH  =  64 is the number of neural channels; and Δt_mo and Δt_og are the temporal offsets for each type of event. N_TP, N_TP · Δt (history used to sample neural features), Δt_mo and Δt_og were used as decoder parameters—their values were selected from a following list of values: N_TP: 3 and 5; N_TP · Δt: 0.3 s and 0.5 s; Δt_mo and Δt_og: any value from  −200 ms to 200 ms with 0.5 ms steps. We additionally formed another class C_OTHER representing states at least 10 ms away from all mo and og events.

We then calibrated a set of mrLDA decoding models using C_mo, C_og and C_OTHER and all possible combinations of parameter values (Milekovic et al 2013a, Capogrosso et al 2016, Milekovic et al 2018). We used the mrLDA regularization parameter as an additional parameter with values of 0, 0.001, 0.1, 0.3, 0.5, 0.7, 0.9 and 0.99. The performance of each of the models was validated using four-fold crossvalidation on the training dataset as follows. A model was calibrated on three quarters of the training dataset and tested on the remaining part of the training dataset. For each time point of this remaining part, the mrLDA model calculated the probability of observing neural activity belonging to C_mo and C_og classes. When one of these two probabilities crossed a threshold of 0.9, the decoder detected movement onset or object grasp, mo_DET or og_DET, respectively. To reproduce the sparsity of these events, we ignored the probability values of the detected event for 1 s after the detection. We pooled the time series of actual and detected events across all four folds and calculated normalized mutual information using a tolerance window of 200 ms (Milekovic et al 2013a, Capogrosso et al 2016, Milekovic et al 2018). We then calibrated another decoding model on the left-out testing dataset using the parameter values that resulted with the maximum normalized mutual information. This decoder was then used on the testing part of the dataset to detect the mo_DET or og_DET events. We again pooled actual and detected events across all five testing folds and measured the decoding performance using temporal accuracy—the ratio of actual events that have one and only one event of the same type within the 200 ms neighbourhood. We estimated the standard error of the temporal accuracy using bootstrapping with 10 000 resamples (Moore et al 2009).

We applied a similar approach to classify grasp types (Mk-Jo) and reach target locations (Mk-Cs) using mrLDA from the neural recordings from either the motor or somatosensory cortex (Milekovic et al 2015). We measured the classification performance using leave-one-out cross-validation (Hastie et al 2009)—in each validation fold ('e'), a different trial e is selected as a test dataset and all other trials are used as a training dataset to calibrate a mrLDA model. This calibration involves a procedure to select the regularization parameter of the mrLDA model from the values of 0, 0.001, 0.1, 0.3, 0.5, 0.7, 0.9 and 0.99. The training dataset is used to derive the classes of neural features, C₁(e,t), ..., C_N(e,t) (N  =  4 for four objects in Mk-Jo, and N  =  3 object positions for Mk-Cs) as described above for event timing using leave-one-out validation on the training dataset. We then calibrated a mrLDA classifier on the complete testing dataset using the regularization coefficient that resulted with the maximum decoding accuracy. We pooled the classification outcomes across all folds and measured the decoding accuracy. This procedure resulted in a confusion matrix and decoding accuracy value for every value of the latency t around both movement onset and object grasp events. We estimated the standard error of the decoding accuracy using bootstrapping with 10 000 resamples (Moore et al 2009).

We assessed continuous predictions of limb kinematics using a Kalman Filter (Wu et al 2006). The Kalman Filter provides a probabilistic framework to predict the state of the limb during the reach and grasp task based on the neural recordings. The output of the filter was the state of the limb. For Mk-Cs, this comprised the position and speed of the elbow and wrist. For Mk-Jo, it comprised the position and speed of the thumb, index and pinkie finger joints (distal phalanx, PIP joint and MCP joint), as well as of the wrist, elbow and shoulder joints. We computed the limb kinematics and the instantaneous firing rate of each neuron at 50 ms intervals. We used the multiunit firing rate obtained from thresholding on M1 and S1 arrays as inputs to the decoder. The neural signals were shifted relative to the kinematics by a static value of 100 ms in Mk-Cs and 70 ms in Mk-Jo. These lags were determined by testing the models on various delays to optimize the model performance. We trained and tested the models using the leave-one out cross-validation method: we iteratively set aside one trial for testing and trained the model using the remaining trials. We pulled together kinematic and neural data over different objects (small sphere and large sphere) for Mk-Jo and over different positions (left, center and right) for Mk-Cs. The performance was assessed using the coefficient of determination R² for which we computed the 95% confidence intervals across all repetitions. 3D static hand posture reconstruction represents hand configuration for the best time point of the fold resulting in the highest average R² for both the small sphere and the cylinder objects in Mk-Jo.