Action kinematics as an organising principle in the cortical control of human hand movement

Hand movements are controlled by neuronal networks in primary motor cortex (M1). The organising principle encoding hand movements in M1 does not follow an anatomical body map, but rather a distributed representational structure in which motor primitives are combined to produce motor outputs. Electrophysiological recordings in primates suggest that M1 neurons encode kinematic features of movements, such as joint position and velocity. Human imaging data concur: relative differences in movement kinematics are mirrored by differences in the associated patterns of M1 activity. However, M1 exhibits well-documented sensory responses to cutaneous and proprioceptive stimuli, raising questions regarding the origins of kinematic motor representations: are they relevant in top-down motor control, or are they an epiphenomenon of bottom-up sensory feedback during movement? Here we show that the kinematic signature of a wide variety of naturalistic hand movements is encoded in human M1 prior to the point of movement initiation. Using a powerful combination of high-field fMRI and MEG, a spatial and temporal multivariate representational similarity analysis revealed that patterns of M1 activity mirrored kinematic, but not muscle-based features of naturalistic hand movements prior to movement onset. Comparable M1 activity was not observed for an ethological action model based functional mappings proposed in M1. Our spatial and temporal analyses provide firm evidence that the top-down control of dexterous movements activates cortical networks in M1 encoding hand kinematics.


Main text
Mounting evidence supports the encoding of movements in M1 based on kinematics and synergistic muscle activation, rather than the anatomy of the peripheral musculature 4,5 . Measurements from individual M1 neurons in non-human primates reveal the encoding of multiple kinematic features, such as speed, direction, and position in the same cells in a time-varying manner 6 .
The same neuronal populations have been shown to encode instantaneous features during motor execution, as well as the target kinematic end point and upcoming movement trajectory 7,8,9,10 .
In the human brain, evidence of neuronal tuning to multiple kinematic features has been reported during the production of intended movements from M1 microelectrode recordings made in tetraplegic patients 14 . The encoding of kinematic features of hand movements in M1 has also been supported by human imaging studies. Patterns of fMRI activity in sensorimotor cortex have been shown to mirror the relative differences in the final joint configuration across a range of prehensile movements 11 . Similarly, the representational structure of fMRI activity in M1 during finger flexion is consistent with patterns of finger co-use during naturalistic hand movements 3 .
However, the functional relevance of kinematic encoding in M1 to human motor control remains a fundamental unknown. As well as their role in motor output, M1 neurons exhibit rapid and integrative responses to somatosensory signals 12,15 . Kinematic information is inextricably linked to proprioceptive and tactile signals: specific patterns of movement are associated with specific patterns of sensory feedback. Are kinematic motor representations reported in human M1 functionally relevant in the process of top-down motor control, or an epiphenomenon generated by bottom-up sensory feedback during human movement production?
We addressed this question using a spatiotemporal multivariate representational similarity analysis to ask where in the human brain and when during movement production are the kinematics of human hand movements encoded? This approach combined high-field fMRI and MEG data with kine-  matic data glove recordings made during a broad repertoire of prehensile and non-prehensile hand movements. Probing recordings of human brain activity with high spatial resolution from fMRI and high temporal resolution from MEG offers a powerful means to identify the location and timing of kinematic information encoding. Together this information was used to dissociate the relevance of kinematic information in M1 to top-down or bottom-up processes in motor control, as well as a relevance of alternative muscle-based or ethological action based models.
Ten right-handed participants performed a range of 26 prehensile and nonprehensile hand movements 16,17 (Table 1, Video S1) in two fMRI sessions (1.5 hours total fMRI data per participant), two MEG sessions (1.5 hours total MEG data per participant), and a behavioural testing session (35 minutes kinematic data recording). In each session participants wore a right-handed 14-channel fibre optic data glove; kinematic data were recorded throughout all sessions. Electromyography (EMG) data were acquired during MEG sessions to validate the movement onset measures calculated from the data Figure 1: Spatial and temporal evidence for the encoding of hand kinematics in primary motor cortex prior to movement production. Kinematic, muscle, and ethological action models of hand movement were used in a spatiotemporal representational similarity analysis. Top row: fMRI data show that kinematic information was encoded consistently in primary motor cortex across all ten participants; complementary MEG data revealed temporal encoding of kinematic information (blue box) prior to movement onset (green line). The muscle model (middle row) and ethological action model (bottom row) showed very limited evidence of encoding, outside of M1 in the post-central gyrus and offered no evidence of significant temporal encoding during movement production. MEG temporal searchlight plots: data presented are from beta band analysis; full analysis presented in Figure 4, green line -movement onset defined by the data glove; blue regions -significant peaks in representational similarity between MEG data and the motor model; dashed line -correlation noise ceiling. EMG onset violin plots based on data presented in Figure S10. Model matrices reproduced in a larger format in Figure S2. show consistent encoding of kinematic information in the left motor cortex, contralateral to movement. Heatmaps were constructed from individually thresholded cortical searchlights for each participant, derived using their own kinematic model (C) (Omnibus threshold, α = 0.01, maximum accuracy distribution calculated from peak correlation value across 10,000 searchlight permutations with label-switching). Supra-threshold range of Spearman's ρ for each participant presented in Table S2. glove.
To probe the spatial and temporal correspondence between patterns of brain activity and hand kinematics, data glove recordings were used to construct a kinematic model quantifying the similarity of the kinematic signals measured during each of the 26 movements ( Figure 1: Top row, Figures S2 and S3).
The kinematic model quantified the distance between the displacement measures for each movement pair across the 14 channels (Pearson correlation), subject to a Fisher Z-transformation and averaged across the 14 recording channels. The resulting kinematic model exhibits strong split-half and intersession consistency within participant ( Figure S1). A grand average of the kinematic model across sessions and participants was subject to non-classical multidimensional scaling for visualisation of the relative dissimilarity of each movement across two dimensions (Video S2). In both the spatial and temporal representational similarity analysis, the kinematic model was investigated 6 . CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. : MEG temporal representational similarity analysis searchlight in motor cortex reveals encoding of kinematic information prior to movement onset. Temporal MEG searchlight analysis reveals evidence of the kinematic encoding of hand movements prior to movement onset (green bar). Distinct significant peaks in the correspondence between the kinematic model and MEG data (blue) were observed in the beta band (-210 ms to -85 ms) and the alpha band (-175 ms to -115 ms). An additional significant peak in the beta band analysis was observed after movement onset (1260 ms -1350 ms). No such significant peaks were observed for the muscle model or ethological action model. Green line -movement onset defined by the data glove; blue regions -significant peaks in representational similarity between MEG data and the model; dashed line -correlation noise ceiling.
alongside two other models. A muscle based model was constructed from high-density EMG recordings (15 channels) made in an independent cohort of 10 participants performing the same range of hand movements (Figure 1: Middle row). An additional ethological action model classified movements into precision prehensile, power prehensile, and non-prehensile, based on the notion of ethological maps in primate M1 13 (Figure 1: Bottom row).
We first used high-resolution fMRI data to perform a cross-validated cortical surface-based searchlight representational similarity analysis to find evidence for the spatial encoding of kinematic information during movement.
In each participant and each cortical searchlight, the unsmoothed pattern of fMRI activity during movement was used to construct a representational dissimilarity matrix (RDM) 18  This analysis assessed where the relative dissimilarities in the kinematic recordings across the different hand movements were mirrored by the relative differences in the pattern of fMRI activity elicited by performing the same movements. The searchlight revealed a strong and consistent representational similarity in the contralateral pre-central region of the anatomical hand-knob 19 across participants ( Figure 1). Specifically, the fMRI searchlight results revealed the consistent encoding of the kinematic information in Brodmann Area 4 during the production of hand movements across participants (Table 2) 20 .
Inspection of the single-subject cortical searchlight results for the kinematic model highlights the consistent and spatially limited correspondence of the kinematic model and fMRI data at the level of individual participants and models in contralateral M1 (Figure 2A). A highly comparable result was also observed using the kinematic model constructed from the data glove recordings made in the behavioural testing session ( Figure S4), highlighting the applicability of this result to real-world hand use in an upright sitting position. No such consistent representational similarity was observed in the corresponding searchlight of movement-related activity in the ipsilateral hemisphere at the group level ( Figure 2B and Figure S4B). Ultra high field fMRI data analysed at the level of individual subjects offered detailed spatial resolution, revealing the encoding of kinematic information in the hand knob region of M1. However, fMRI offers relatively poor temporal resolution to understand the point in time at which the kinematic model matches the pattern of brain activity. The boundary between motor and somatosensory cortex is increasingly blurred by evidence of sensory processing in M1 12 and motor modulation of sensory afferents 21 . The observed representational similarity of fMRI activity and hand kinematics may result from top-down control of motor function or bottom-up proprioceptive information passed back to M1 and S1. In order to dissociate the driving force behind the kinematic model fit observed in the fMRI data, a temporal representational similarity analysis of MEG data was used to identify the point during movement preparation or execution at which kinematic information is encoded in the M1.
A cross-validated fixed-effects representational similarity analysis was applied, comparing a group average kinematic model derived from data glove recordings made in the MEG scanner to the pattern of alpha (7-14 Hz), beta (15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30), and gamma (30-100 Hz) band MEG brain activity in M1 ( Figure S9) in 20 ms sliding windows during movement preparation and execution. The muscle model and ethological action model were assessed in equivalent analyses. In light of the interest in contrasting the kinematic and muscle models, the kinematic and muscle models were each assessed in a partial correlation to discount the contribution of the other.
In both the alpha and beta band analysis, there was significant correspondence between the MEG data and the kinematic model both preceding movement onset, and in the case of the beta band, after movement onset ( Figure 4). In the beta band, the kinematic model mirrored the pattern of brain activity in a significant peak from -210 ms to -85 ms relative to movement onset (peak Spearman's ρ: 0.32). There was also a significant peak in the correspondence between the kinematic model and MEG data in the beta band after movement onset in the beta band (1260 -1350 ms; peak Spearman's ρ: 0.35).
In the alpha band a correspondence between the kinematic model and MEG data was observed prior to movement onset from -175 ms to -115 ms (peak Spearman's ρ: 0.37). No significant peaks were observed in the correspondence between the M1 MEG signal in the alpha, beta, or gamma band and either the muscle model or the ethological action model (Figure 4).

10
. CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under   . CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019.
kinematic information occurs prior to the onset of a movement. In other words, the relative differences in the kinematic structure of a range of different hand movements is encoded in M1 up to 210 ms before the onset of movement can be detected in the hand.
Information contained in the kinematic model showed temporally distinct correspondence to information contained in the alpha and beta bands of the MEG data. From 210 ms to 90 ms before movement is detected, the representational structure in the M1 beta band corresponds significantly to the representational similarity of the kinematics of the upcoming movement.
The correspondence between the kinematic model and the information con- where ERD in M1 corresponds to increased activation in the region 22 , with post-motion event related synchronisation in M1 26 . Here we demonstrate that there is a link between information contained in the alpha and beta frequencies in M1 before movement onset and the subsequent kinematics of hand movements (Figures 1 and 4), suggesting that important information about the upcoming motor command may be encoded within these oscillations 25,27 . A significant peak was also observed post-movement beta-band analysis; it be could be speculated that this peak represents information encoding of kinematic information around the time PMBR is known to occur after movement, or could reflect afferent inputs to M1 28 ; these two possibilities cannnot be dissociated by the current experiment.
In contrast to alpha and beta frequencies, we observed no concurrence between the information contained in the gamma-frequency and the kinematic model. An increase in the amplitude of gamma oscillations has previously been reported during motor execution: movement-related gamma synchronisation (MRGS) 29,30 . In contrast to alpha and beta frequencies, evidence from studies of gamma oscillations report changes only after movement onset, and therefore would not be implicated in the encoding of information in M1 prior to movement, consistent with the data herein 31,32,33 .
Hand kinematics have previously been investigated in the context of human fMRI. Relative differences in target joint position at the end of a hand movement have been shown previously to mirror the relative differences in the fMRI signal in a broad region sensorimotor cortex 11 . Additional work considering unidigit and multidigit flexion has demonstrated that patterns of M1 fMRI activity associated with such movements are better explained by kinematic models of digit co-use than by competing muscle-based models 3 . In the present study we have used MEG to fundamentally extend on these findings, demonstrating a top-down role for kinematic encoding prior to movement onset. Furthermore by using a kinematic model that compared the displacement trajectory of each movement rather than a single joint position, it was possible to contrast the kinematic features of both prehensile and non-prehensile movements to explore cortical encoding relevant to a full range of naturalistic hand use. Moreover, by capitalising on the gains in spa-

13
. CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint tial and temporal resolution afforded by using 7T fMRI and a limited field of view, we have been able to firmly pinpoint the spatial location at which kinematic information is encoded to the motor region of the anatomical hand knob, corresponding principally to Brodmann Area 4 19 . The fMRI spatial searchlight analysis did not reveal evidence of consistent encoding of kinematic information in ipsilateral M1 across participants ( Figure 2). Previous fMRI studies provide evidence for the activation of ipsilateral M1 during the production of individual uni-digit movements 34,35 but not multi-digit sequences of uni-digit movements 36 . This study considered a broad array of naturalistic hand movements, engaging a wide variety of hand kinematics, involving simultaneous and/or sequential movement of different digits. It is possible that unlike sequences of uni-digit movement, these more complex movements do not drive the circuits of ipsilateral M1 as uni-digit movements do 34,35 .
Previous studies have made direct comparisons between muscle-based models and kinematic models, arguing for the latter as an organising principle in the encoding of hand movements 3,11 . Here, a model constructed from an independently acquired set of high-density EMG recordings did not reveal any evidence for the spatial or temporal encoding of information on the basis of differences in muscle activity across the range of 26 hand movements under study. In addition, the kinematic model showed a superior representational similarity to the fMRI dataset than the muscle model in primary motor cortex ( Figure 3). As with previous studies, these findings do not rule out the existence of muscle representations in M1, but rather support the existence of highly organised muscle representations structured around movement kinematics rather than anatomy. The assertion perhaps explains the fractures and repetitions observed in muscle representations during the search for an M1 body map 2 .
14 . CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint frequency band 42 , potentially consistent with the notion that during observation of biological motion, M1 may encode kinematic information.
The data presented in this study rely on complementary information acquired from BOLD fMRI and MEG, though the remit of this work does not extend to fusion of the two modalities. BOLD fMRI provides only an indirect measure of neuronal activity based on haemodynamic changes associated with the execution of a given task 43 , which can be resolved with a relatively high degree of spatial specificity with 7T imaging. In contrast, MEG reflects a more direct, temporally-rich, measure of neuronal activity.
While the origins of the measured signals differ, compelling recent evidence provides non-coincidental data to support the notion of shared information across MEG and fMRI measures of brain activity across a wide range of frequency bands 44 ; similar correspondences have been reported from invasive electrocorticography data 45 . However, the spatial component of MEG data must be inferred from mathematical modelling. Despite advances in the context of MEG source localisation, this feature of MEG analysis limits the spatial specificity of the measured signals, which integrate information across relatively large tissue volumes in comparison with fMRI 46 . We have harnessed the spatial and temporal strengths of fMRI and MEG, which in combination provide greater insight regarding the encoding of movements in M1 than the sum of their individual parts.
Here we apply a rich multi-modal design with multivariate analysis to demonstrate that the encoding of kinematic information in human M1 occurs prior to the onset of a wide range of naturalistic hand movements, in contrast to competing muscle and ethological action models. Mounting evidence for the encoding of complex kinematic information in M1 from this and other work continues to blur the boundary between primary somatosensory and primary motor cortex: even M1 neurons have been shown to rapidly con-16 . CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint solidate sensory torque information across multiple joints 15 . The notion of kinematic representation in M1 is compatible with recent evidence of the tight integration of information across the central sulcus 47 , whereby S1 encodes the current body state, while M1 encodes the kinematics necessary to achieve the intended body state. Such a system of motor control would see kinematic information encoded prior to movement onset as a prediction for the future sensory inputs expected by S1 when a movement has been achieved 48 . . CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint

18
. CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint . CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019.

22
. CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint . CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint . CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint . CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint . CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint . CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019.

Motor task and kinematic data acquisition
During all sessions participants were engaged in a motor task involving the production of a range of 26 hand movements (Table 1)  CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint the MRI scanner, and in the MEG scanner, as detailed below.
None of the grasping tasks in this study engaged participants with real objects; previous work has differentiated motor activity with or without real objects in anterior intraparietal sulcus, but not primary motor cortex: as such an object-free study design seemed appropriate for a study focusing on M1 54 .

Kinematic recording session
During the behavioural testing session participants performed five runs of the motor task. Participants were seated at a desk with their right forearm supported on a memory foam mount, while wearing the data glove. Participants viewed instructions presented on a 14 inch laptop display. Each movement block comprised a 3 second video of the movement to be produced, a 1 second preparation period and 8 subsequent movement trials; each comprising 1.6 seconds of movement (green expanding/contracting bar), followed by a 0.8 second rest period (red static bar). The transition of the bar from red to green was defined as the go signal. A break period of up to 15 seconds was permitted between each movement block; participants advanced the task with a key-press using their left hand. Excluding break periods each task run was 10 minutes and 3.2 seconds in duration. The four task runs yielded 33 minutes and 16.8 seconds of kinematic data recording per participant. The kinematic dissimilarity matrices were averaged across task runs withinparticipant to yield an fMRI, MEG, and behavioural kinematic model for each participant. The split-half consistency and inter-session consistency of these models is outlined in Figure S1. A grand average across participant and session kinematic model (Figures 1 and S2) was computed and subject to hierarchical clustering; this resulting clustering was applied to all visualisations of the kinematic model.

Muscle model
An independent EMG dataset was acquired in order to construct a model of hand movement dissimilarity on the basis of muscle activity in the hand.
An independent cohort of ten participants (Age range: 20-30; mean age: 31 . CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint 25.1; age SD: 3.57; 5 female) undertook a more detailed EMG recording than was feasible during the MEG session, while performing the same 26 hand movement. EMG data were acquired using a Biosemi Active 2 system with a 32 channel headbox (Biosemi B.V. Amsterdam). Muscle activity was recorded using touchproof flat active electrodes. Electrodes 1-15 were placed as labelled in Figure S16, electrode 16 was used to rereference the EMG data in subsequent analysis and was placed on the bony protrusion of the elbow. There were also CMS and DRL electrodes, which served as a ground/reference during recording in the Biosemi software; they were placed on the palmar side of the wrist. The EMG data were recorded at 2048Hz.
The EMG recording sessions mirrored the design and setup of the kinematic recording session outlined above and were informed by previous fMRI MVPA studies of digit flexion 3 . Five runs were recorded in total, each containing 26 trials (one for each of the movements). The EMG data were processed using Fieldtrip 55 . EMG data were rereferenced to electrode 16, rectified and low-pass filtered (fourth order Butterworth filter: 40 Hz), and epoched relative to earliest measured muscle onset in any EMG channel using an adaptive threshold (activity duration threshold: 200ms) (Hooman Sedghamiz: Matlab File Exchange: Automatic Activity Detection in Noisy Signals using Hilbert Transform.) This results in individual trials of 2.0s in duration. These trials were baselined using the fixation cross window at the start of each trial. EMG trial data were then subject to multivariate noise normalisation by weighting channels in trial by the error covariance across the different channels in order to more accurate quantify the true differences between the muscle activity across different movements. The normalised trials were averaged into 5 folds. A Mahalanobis distance comparing each of the muscle activity of each of 26 different movement types was calculated using a cross-validated leave-one-out approach. In each iteration, the muscle 32 . CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint activity patterns from one fold were assigned to fold A and the muscle activity data from the remaining four folds were assigned to fold B; distances were calculated between all possible pairs of the 26 movement muscle activity recordings across these two folds (Equation (3)). Distance measures were calculated across all possible pairs of cross-validation folds. An average muscle model across all ten participants' data was generated and used to probe the spatial and temporal encoding of muscle based dissimilarities in the brain using fMRI and MEG (Figures 1 and S2).

Ethological action movement model
An alternative ethological action based model was constructed based on more recent evidence of ethological maps in primate M1 13 , and therefore categorises movements on the basis of their specific action, namely prehensile movements, sub-categorised into precision grip and power grip, and non-prehensile movements 17 (Figure 1). The ethological action model was subject to hierarchical clustering for visualisation.

MRI data acquisition
MR data were acquired using a Siemens 7T Magnetom system (Siemens Healthcare, Erlangen, Germany) with a 32-channel head coil. Blood oxygenation level dependent (BOLD) fMRI was acquired with a T2*-weighted multi-slice gradient echo planar imaging (EPI). True axial slices were po- structural MRI data were acquired to facilitate BOLD EPI slice placement 33 . CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint and for cortical surface reconstruction (TR/TE: 2200/2.82 ms, isotropic resolution: 1.0 mm, GRAPPA factor = 2). An additional gradient echo BOLD EPI acquisition of 4 volumes was acquired using posterior-anterior phaseencoding direction for distortion correction.

fMRI behavioural task
During the fMRI acquisitions participants performed a total of ten runs of the motor task (5 runs per MRI session). Participants were led supine with their right forearm supported against their right hip and their elbow supported by a foam pad, while wearing the data glove. Participants viewed instructions via a mirror mounted on the transmit coil and a projector screen mounted at the end of the bore. Each movement block comprised of a 3 second instruction screen ("Prepare to Move"), a 3 second video of the movement to be produced, and a 1 second further instruction screen ("Move"), followed by 5 movement trials, each comprising 1.6 seconds of movement (green expanding/contracting bar), followed by a 0.4 second rest period (red static bar). Each movement block was 17 seconds. In addition to the movement blocks, 8 rest blocks were included in each task run; rest blocks were of equivalent duration to movement blocks and comprised of a 3 second instruction screen ("Rest"), a 3 second video of a static resting hand, and a 1 second further instruction screen ("Rest"), followed by the same period of expanding and contracting bar visual stimuli as the fMRI movement blocks.
Rest blocks were positioned randomly in each run, excluding self-adjacency.

Structural MRI data preprocessing
MPRAGE data were subject to reorientation, bias-field correction and brain extraction using the FMRIB Software Library (FSL) fsl anat tool 56,57,58 prior to cortical surface reconstruction using FreeSurfer Version 5.3.0 59,60 .

34
. CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint 2.10. fMRI data analysis 2.10.1. fMRI preprocessing and general linear modelling fMRI data were subject to standard preprocessing, including motion correction with MCFLIRT 61 , brain extraction using BET 57 , and high pass temporal filtering (100 second threshold). fMRI data were not subject to spatial smoothing. All fMRI data were subject to manual independent components analysis denoising 62 . Distortion correction was undertaken using FSL Topup to estimate a fieldmap image for use in FSL FUGUE 63 . Undistorted BOLD EPI data were co-registered with structural MPRAGE data using Boundary-Based-Registration from FMRIB's Linear Registration Tool (FLIRT) implemented in epi reg 64,61,65 . Example fMRI timeseries from a single voxel located in the anatomical hand knob is presented for four participants on a single session in Figure S15.
For each participant and each fMRI run, fMRI data were analysed using a first-level general linear modelling (GLM) approach implemented in FSL FEAT 58 using FMRIBs Improved Linear Model (FILM) to estimate time series autocorrelation and pre-whiten each voxel. Each of the 26 movements was modelled with a separate boxcar regressor with gamma-HRF convolution and its temporal derivative, giving a total of 52 regressors. Parameter estimates were calculated, contrasting each movement type against the rest condition; these voxel-wise maps and an estimate of the residuals from the GLM were resampled into the respective participants' structural space and used in subsequent representational similarity analysis (RSA).

fMRI multivariate noise normalisation
In order to account for the spatial structure of the noise inherent to fMRI data, spatial prewhitening of the parameter estimates from each participant 35 . CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint and each fMRI task run was conducted. The residuals (R) from the first-level GLM analysis provided an estimate of data not fit by the model regressors across voxels (V) and time (T), from which a V x V covariance matrix ( Σ) can estimate the noise structure across voxels (Equation (1)) 66 .
The noise covariance structure was combined with the voxel-wise parameter estimates (P) for a given movement type (k) to generate a spatially prewhitened parameter estimate (P * k : Equation (2)): . CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint to improve spatial specificity. Spatially pre-whitened parameter estimates were extracted from the resulting volumetric region corresponding to each searchlight.

fMRI cross-validated distance measures
Within each searchlight the similarity between each of the spatially prewhitened voxel-wise parameter estimates corresponding to each of the 26 different movement types was calculated using a cross-validated leave-oneout approach to avoid the possibility of over-fitting the data 68,69 . In each iteration, the parameter estimate maps from one fMRI task run was assigned to fold A and the parameter estimate maps from the remaining nine task fMRI runs were assigned to fold B; squared Euclidean distances were calculated between all possible pairs of the 26 movement parameter estimate maps across these two folds (Equation (3) The correspondence between the fMRI-derived RDM in each searchlight and the candidate kinematic and theoretical models was assessed using a Spearman's rank correlation, with the resulting ρ (rho) value was plotted in each searchlight's central vertex on the cortical surface. For statistical 37 . CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint inference a fixed effects randomisation test 18 was applied on the individual participant level: correlations using 10,000 condition-label randomisations were undertaken in each searchlight. From each of the permutations, the spatial peak ρ value (rho) was extracted from across the cortical surface, forming a maximum accuracy distribution from which an omnibus threshold (α = 0.01) was extracted. The resulting thresholded ρ-value surface maps for each participant were resampled onto the Human Connectome Project 32k surface (S1200.L.pial.MSMAll.32k fs LR.surf.gii), binarised and used to form a heatmap corresponding to the spatial distribution of the each model fit across participants. In light of the interest in contrasting the kinematic and muscle models, a comparison of the corresponding unthresholded Spearman's ρ cortical surface maps was undertaken using a Wilcoxon signed-rank test (one-sided), subject to FDR correction (α = 0.05).

fMRI motion considerations
Variability in the magnitude of fMRI motion across different movement conditions has the potential to influence the observed pattern of results. The potential for noise induced by participant motion has been mitigated in a number of ways. First, all data has been subject to ICA denoising to remove any characteristic motion artefacts 62 . Second, the multivariate analysis of fMRI data emoloyed herein used spatial prewhitening of the parameter estimates to account for voxel-wise variability in noise to not downweight voxels with high error variance and to account for noise covariance between voxels 66 . Finally DVAR values were calculated for each fMRI timeseries (D: temporal derivative of time courses, VARS: root mean squares variance over voxels). These values quantify for each frame of an fMRI acquisition the magnitude of signal intensity change in comparison in volume N compared with volume N-1, as per the following formula: 38 . CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint Where I i is image intensity at locus − → x on frame i; angle brackets denote the spatial average over the whole brain 70 . DVARs are able to quantify corruption of fMRI acqusition due to head motion. DVAR values were extracted for volumes corresponding to each of the 26 hand movements for all participants; the resulting distribution of DVAR values is presented in Figure S13. The profiles of very limited motion across participants during each session of around 10 minutes in duration also demonstrate high quality data acquisition ( Figure S14).

MEG data acquisition
MEG signals were measured continuously at 1200Hz during the motor task using a whole-head 275-channel axial gradiometer CTF MEG system (CTF, Vancouver, Canada) located inside a magnetically shielded room. An additional 29 reference channels were recorded for noise cancellation purposes and the primary sensors were analysed as synthetic third-order gradiometers 71 . Three electromagnetic coils were placed on three fiduciary locations (nasion, left and right pre-auricular) and their position relative to the MEG sensors were recorded continuously during each experimental block. The head surface and fiducial locations were digitized using an ANT Xensor digitizer (ANT Neuro, Enschede, Netherlands) prior to the MEG recording.

MEG behavioural task
During the MEG data acquisitions participants performed a total of ten runs of the motor task (5 runs per MEG session). Participants were sitting upright with their right forearm and elbow supported on a foam armrest, while 39 . CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint wearing the data glove. Participants viewed instruction on a back-projected screen in front of them from a projector mounted outside the shielded room. Each movement block comprised of a 2 second period with a central fixation cross, a 3 second video of the movement to be produced, and a 1 second instruction screen ("Prepare to Move") followed by five movement trials, each comprising 1.6 seconds of movement (green expanding/contracting bar), followed by a 0.8 second rest period (red static bar). Each movement block was 18 seconds. The order of movement blocks was randomised within each task run; each movement was presented once per task run.

Data glove movement onset detection: MEG sessions
The 14 channels of data glove recordings collected during the MEG sessions were aligned with the MEG acquisitions using the onset of the green bar and were epoched alongside the MEG data. Epoched data glove recordings were subject to onset segmentation using an adaptive threshold (activity duration  CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint (Donders Institute for Brain, Cognition and Behaviour, Radboud University Nijmegen, the Netherlands. See http://www.ru.nl/neuroimaging/fieldtrip).
Co-registration was performed in a two stage process: first the fiducial locations were marked on the T1 structural for that participant; the head digitization data was then used to align the data with the MRI, subject to manual adjustment. Alignment was undertaken independently for data from the two MEG sessions.
Data from each movement type were epoched from the 10 task runs and concatenated into a new dataset containing 10 blocks, each containing 5 movement trials. The fixation cross and movement trials were epoched from the overall block. The movement trials were defined relative to the data glove defined movement onset time (movement trial time: 2 s; pre-onset time: 0.5s, post-onset time: 1.5s). The fixation cross period was used as a baseline for the 5 movement trials within each movement block. A high pass filter of 1 Hz and a low pass filter of 100 Hz were applied. MEG analyses were conducted across three frequency bands: alpha (7-14 Hz), beta (15-30 Hz) and gamma (30-100 Hz). All of the movements trials for a given movement type were concatenated across the 10 task runs, creating a dataset comprising 50 repeats of a movement. At this point the data was visually inspected and those trials containing artefacts were removed from further analysis up to a maximum of 10 trials, such that the minimum number of movements trials per movement included in further analysis was 40.

MEG source reconstruction
In order to reconstruct oscillatory activity at brain locations directly comparable across participants, the individual anatomical MRI was non-linearly warped to the MNI MRI template. The MNI template was divided into a 41 . CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint warp was used to warp the template grid into the anatomical space of each participant. Sensor leadfields were calculated using a semi-realistic volume conduction model based on the individual anatomy 72 . The temporal evolution of source activation at each location in the brain was estimated using a linearly constrained minimum variance (LCMV) beamformer algorithm 73 with the optimal dipole orientation at each voxel estimated using singular value decomposition (SVD). Virtual sensors were then reconstructed from all 3294 voxels by multiplying the sensor level data by the corresponding set of optimised weights. At this stage data were subject to multivariate noise normalization 74,75 : we calculated the error covariance matrix at sensor level and then used this combined with the filters from the LCMV to create the virtual sensor data. This means that sensors with more noise would be down-weighted compared to those with less noise. At this stage the data was also down-sampled to 600 Hz to reduce computational cost.

MEG temporal representational similarity analysis
The MEG data were split to produce 10 partitions and then averaged within each partition to perform a cross-validated representational similarity analysis to avoid the possibility of over-fitting the data 68 Motion parameters for all MEG acquisitions were extracted and analysed to rule out the possibility of excessive head motion as a potential driving force behind any observed patterns of brain activity. Rotational and translational displacement for each participant and each experimental session are presented in Figure S11. In addition, the motion parameters during each movement block were extracted and the resulting distribution is presented across the 26 different movement types ( Figure S12). The profiles of motion across participants demonstrate a high quality data acquisition.

Electromyography with MEG
Electromyography (EMG) data were acquired simultaneously with MEG data. Three surface EMG electrodes were attached to the right hand underneath the data glove, positioned on abductor pollicis brevis (APB), first 44 . CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint dorsal interosseus (FDI) and abductor digiti minimi (ADM). The area under the electrodes was exfoliated and cleaned with alcohol prior to data acquisition. EMG signals were recorded at 1200Hz. EMG data were initially subject to a bandpass filter (20-500Hz) and a notch filter (50 Hz). EMG data were epoched and baselined alongside the MEG data. Epoched EMG data were subject to manual artefact rejection. Signals from the three electrodes during each epoch were independently subject to a Hilbert transform and smoothing (5 ms window) prior to activity onset segmentation using an adaptive threshold (activity duration threshold: Video S1: Compilation of instructional videos used at the beginning of movement blocks in all testing sessions. Movement labels are provided for reference only; labels were not included during the task (VideoS1.mov).
Video S2: Visualisation of multidimensional scaling of grand average kinematic model constructed across participants and sessions. (VideoS2.mov).

45
. CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under  . Figure S1: Data-driven kinematic models constructed for each participant and each session type exhibit strong split-half and inter-session consistency. Muscle model reproducibility data also presented.

46
. CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under

47
. CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint Figure S2: The average kinematic model across participants for each session type, the muscle model derived from an independent cohort of participants, and the categorical ethological action model.

48
. CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint

49
. CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint Figure S4: fMRI searchlight analysis conducted using the kinematic model constructed from data glove recordings made during the behavioural testing session. Evidence of the encoding of kinematic data in contralateral primary motor cortex persists using independent data glove recordings while participants were sitting upright in a more naturalistic position. Comparison with data presented in Figures 1 and 2. Supra-threshold range of Spearman's ρ for each participant presented in Table S2.

50
. CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint

51
. CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint Figure S6: Individual participant cortical searchlight results using an muscle model of movement encoding. Cortical heatmaps of the left (A) and right (B) hemisphere, show limited but consistent encoding of an action model in the left post-central gyrus, contralateral to movement. Heatmaps were constructed from individually thresholded cortical searchlights for each participant using a single categorical action model (C) (Omnibus threshold, α = 0.01, maximum accuracy distribution calculated from peak correlation value across 10,000 searchlight permutations with label-switching.) Supra-threshold range of Spearman's ρ for each participant presented in Table S2.

52
. CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint Figure S7: Individual participant cortical searchlight results using an ethological action model of movement encoding. Cortical heatmaps of the left (A) and right (B) hemisphere, show limited but consistent encoding of an action model in the left post-central gyrus, contralateral to movement. Heatmaps were constructed from individually thresholded cortical searchlights for each participant using a single categorical action model (C) (Omnibus threshold, α = 0.01, maximum accuracy distribution calculated from peak correlation value across 10,000 searchlight permutations with label-switching.) Supra-threshold range of Spearman's ρ for each participant presented in Table S2.

53
. CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint

54
. CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint

55
. CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint

56
. CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint Figure S11: MEG motion timeseries across participants and sessions. Blue: X-axis, orange: Y-axis, yellow: Z-axis 57 . CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint

58
. CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint

59
. CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint Figure S14: fMRI motion displacement plots. Plots of absolute (blue) and relative (orange) motion calculated using FSL MCFLIRT for each participant and each fMRI task run; motion correction was undertaken prior to ICA denoising.

60
. CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint Figure S15: Example fMRI axial slices and single voxel timeseries from 4 participants. fMRI timeseries data extracted from a single voxel (red crosshairs) for four participants. Data presented were subject to high-pass filter (100 seconds). fMRI data were not subject to spatial or temporal smoothing.

61
. CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint Figure S16: EMG recordings from an independent cohort used to generate a muscle model of hand movement. A) Schematic demonstrating electrode placement for EMG recordings on both the palmar and dorsal surface of the hand and forearm covering these muscles: 1 first dorsal interosseus (FDI), 2-3 dorsal interosseus muscles, 4 abductor digiti minimi, 5 abductor pollicis brevis (APB), 6 adductor pollicis, 7-9 lumbrical muscles, 10 flexor carpi ulnaris, 11 flexor carpi radialis, 12-14 flexor digitorum superficialis and flexor digitorum profundus, 15 flexor pollicis longus. B) Example EMG trace from one participant for three different movements (squeeze, grip sphere and finger abduct) showing electrodes 1-15 across the 2s movement time window.

62
. CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint Figure S17: Noise ceiling calculation for spatial searchlight using fMRI data. To assess the spatial consistency in RDMs calculated from fMRI data at each vertex, each participant's RDM was correlated with the average cross-subject RDM; the correlations were then averaged to obtain a vertex-wise upper bound of the noise ceiling.

63
. CC-BY 4.0 International license a certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was not this version posted April 24, 2019. ; https://doi.org/10.1101/613323 doi: bioRxiv preprint