Research ArticleNEUROSCIENCE

The somatosensory cortex receives information about motor output

See allHide authors and affiliations

Science Advances  10 Jul 2019:
Vol. 5, no. 7, eaaw5388
DOI: 10.1126/sciadv.aaw5388

Abstract

During voluntary movement, the somatosensory system not only passively receives signals from the external world but also actively processes them via interactions with the motor system. However, it is still unclear how and what information the somatosensory system receives during movement. Using simultaneous recordings of activities of the primary somatosensory cortex (S1), the motor cortex (MCx), and an ensemble of afferent neurons in behaving monkeys combined with a decoding algorithm, we reveal the temporal profiles of signal integration in S1. While S1 activity before movement initiation is accounted for by MCx activity alone, activity during movement is accounted for by both MCx and afferent activities. Furthermore, premovement S1 activity encodes information about imminent activity of forelimb muscles slightly after MCx does. Thus, S1 receives information about motor output before the arrival of sensory feedback signals, suggesting that S1 executes online processing of somatosensory signals via interactions with the anticipatory information.

INTRODUCTION

Tactile sensory information from the world is gained through active exploration with the hand. During voluntary limb movement, the central nervous system is continuously inundated with both somatosensory signals arising from changes in the external world (exafference) and those from our actions (reafference). Theoretically, the forward model proposes that the sensory system receives a copy of motor commands to extract exafference from sensory inputs (1). Human psychophysical experiments showed that voluntary limb movement decreases sensitivity to tactile stimulation on the moving limb (2, 3). Electrophysiological experiments in animals demonstrated that somatosensory-evoked potentials in the primary somatosensory cortex (S1) attenuate during voluntary limb movement more than during passive movement (46) and during the movement preparation period as well (4, 6). These results imply that the movement-generating neural circuitry may affect somatosensory processing. Thus, the central nervous system not only passively receives somatosensory inputs but also actively processes these inputs during voluntary limb movement.

To understand how the central nervous system extracts important information for the agent during voluntary limb movement, it is critical to consider the interactions between the somatosensory and motor systems (7). Neuronal activity in S1 during voluntary movement is quite different in firing pattern from passive movement (810). A functional imaging study has shown that S1 activity observed during voluntary limb movement but not during passive movement correlates with the activity of motor-related areas (11). Artificial excitation of the primary motor cortex (M1) in rodents induces excitatory or inhibitory effects on responses to tactile stimuli (1214) or whisker stimulation (1517) at multiple levels of the somatosensory system. These findings imply that inputs from motor-related areas might modulate the flow of somatosensory signals from the periphery to S1. However, sensorimotor interactions have been examined only by comparing somatosensory processes under different conditions (e.g., voluntary versus passive movements or artificial excitation of motor-related areas versus no excitation) such that inputs from motor-related areas have not been clearly dissociated from somatosensory signals that are ascending from peripheral afferents to reveal online sensorimotor interactions during voluntary limb movement. Thus, how and what information the somatosensory system receives from the motor system during voluntary limb movement have yet to be elucidated.

Here, we simultaneously recorded activities of motor and somatosensory cortices and activities of an ensemble of afferent neurons in the dorsal root ganglion (DRG) of monkeys during reach and grasping movements. Using a decoding algorithm, we clarified the temporal and spatial contributions of activities of the motor cortex (MCx) and peripheral afferents to S1 activity. Before movement initiation, S1 was found to encode imminent muscle activity similar to M1, while during movement, S1 activity reflected an integration of information from MCx and peripheral afferents. These results suggest that S1 receives information about motor output from the motor system to predict the consequence of the action.

RESULTS

Simultaneous recording of central and peripheral neural signals

We trained two monkeys to perform a reaching and grasping task in which they held a button down with the right hand for a predetermined time and subsequently moved it to a target lever (fig. S1A). We simultaneously recorded multiregional physiological signals as described below (Fig. 1 and fig. S1C). Electrocorticography (ECoG) signals in S1 and MCx, including dorsal (PMd) and ventral (PMv) aspects of the premotor cortex and M1, were recorded using a surface electrode array with 32 electrodes (fig. S1B). Electromyography (EMG) signals were obtained with pairs of wire electrodes in 12 and 10 forelimb muscles of Monkeys T and C, respectively. Activities of an ensemble of sensory afferent neurons (25 to 39 units in Monkey T and 11 to 15 units in Monkey C) that innervate the forelimb were recorded with two 48-channel electrode arrays that were implanted into two DRGs in cervical segments C7 and C8 in Monkey T and C6 and C7 in Monkey C. An ensemble of afferent neurons included muscle spindles, tendon organs, and cutaneous receptors. Forelimb joint kinematics (shoulder, elbow, and wrist joints) were recorded with an optical motion capture system.

Fig. 1 Simultaneous recording of cortical and peripheral activities showed movement-related modulation of these activities.

(A) An example of simultaneous recording of three trials by Monkey T. Top row and second row: Power spectrograms of the respective electrodes in M1 and S1. Third row: EMG signal from the shoulder muscle. Fourth row: The instantaneous firing rate of an ensemble of peripheral afferents. Bottom row: Three forelimb joint angles along the extension-flexion axis (shoulder, solid line; elbow, dashed line; wrist, dotted line). (B) Modulations of cortical and peripheral activity in Monkey T aligned to movement onset. Top and second: Modulations of high-γ activity in M1 and S1. Third: EMG signals of 12 forelimb muscles. Fourth: Instantaneous firing rate of a neural ensemble of peripheral afferents. Bottom: Six forelimb joint angles (shoulder, solid lines; elbow, dashed lines; wrist, dotted lines). Thin lines represent the activity in each electrode (M1 and S1), units (afferent), and muscles, and thick lines represent their respective averages. The vertical solid, dotted, and dashed lines represent times of the onset of movement, pulling the lever, and the end of the movement, respectively. Arrowheads represent the peak times of respective activities. AU, arbitrary units.

The simultaneous recordings enabled us to directly compare the temporal profiles of movement-related modulations of multiregional physiological signals (Fig. 1B and fig. S1C). It is generally accepted that cortical high-γ activity represents activity of the neuronal population beneath the electrode (18). High-γ activity (60 to 180 Hz) in M1 and S1 encoded the kinematic variables of movements (fig. S2). When we aligned multiregional physiological signals to the timing of movement onset, cortical high-γ activity, activity of forelimb muscles, and neuronal firing of peripheral afferents each exhibited their own respective temporal patterns (Fig. 1B and fig. S1C). High-γ activity in S1 showed three peaks (Fig. 1B, second row): the first peak around movement onset (arrowhead), the second peak around lever pulling (double arrowhead), and the third peak around the end of the movement (triple arrowhead). Most of the peripheral afferents exhibited a large peak of activity around lever pulling (Fig. 1B, double arrowhead in the fourth row). This peak occurred slightly earlier than the second peak of S1 activity, suggesting that activity in peripheral afferents might be transmitted to S1 around the time of pulling the lever. On the other hand, all three peaks of high-γ M1 activity occurred earlier than the corresponding S1 activity peak (Fig. 1B, first row). M1 activity also increased before the onset of S1 activity, as well as that of PMd and PMv (Fig. 1B and fig. S1, C and D), implying that MCx activity might be transmitted to S1.

Decoding of S1 activity

In the above analyses, only the temporal patterns of averaged movement-related modulations were compared between different physiological signals. By a simple comparison of the temporal patterns among different regions (Fig. 1B and fig. S1, C and D), it is not possible to elucidate the temporal and spatial profiles of signal integration. To investigate how preceding activity in MCx, including M1 and premotor cortices, and that in an ensemble of peripheral afferents account for S1 activity, we decoded S1 activity from combined activities in MCx and peripheral afferents using a sparse linear regression (SLiR) algorithm. We built a linear model using combined activities in MCx and peripheral afferents to account for subsequent S1 activity. The model accurately reconstructed the overall temporal pattern of activity in S1 better than a model built from surrogate shuffled control data [Monkey T (mean, n = 16 signals): correlation coefficient Rdata = 0.61, Rshuffled = 0.02; Monkey C (mean, n = 12 signals): Rdata = 0.53, Rshuffled = 0.00; Fig. 2, A and B]. The reconstruction accuracy of models built from combined activities in MCx and peripheral afferents was superior to that of models built from activity in MCx alone or from activity in peripheral afferents alone (fig. S3). These results indicate that both MCx and peripheral afferent activities were necessary information for the reconstruction of S1 activity.

Fig. 2 Combined activities in MCx and afferents account for S1 activity.

(A) Reconstruction of S1 activity using combined activities in MCx and peripheral afferents of Monkey T. Examples of two successive movements. Black line, the observed activity in S1; red line, the reconstructed activity; R, correlation coefficient between the observed and reconstructed activities; vertical solid lines, movement onset; vertical dashed lines, end of movement. (B) Mean decoding accuracy pooled across ECoG electrodes in S1. The correlation coefficient between the observed and reconstructed traces from the recorded data compared with the correlation coefficient between the observed traces and the traces reconstructed from random shuffling of activity (Monkey T: mean, n = 16 signals; Monkey C: mean, n = 12 signals; *P < 0.001). Superimposed bars, mean. (C) Average modulation of the observed S1 activity in Monkey T (S1, gray) and the reconstruction using combined activities in MCx and peripheral afferents (MCx + Afferent, purple) aligned to movement onset. Vertical solid lines, movement onset; vertical dashed lines, average time of end of movement; solid lines, means; shading, SD. (D) Variance accounted for (VAF) between the observed and reconstructed traces throughout the duration of S1 modulation (−100 to 1400 ms around movement onset for Monkey T, −100 to 1500 ms for Monkey C).

We then aligned reconstructed waveforms to the point of movement onset to examine whether the model from the combined activities reconstructed movement-related modulation of activity in S1 (Fig. 2, C and D). The average high-γ activity recorded in S1 increased before movement onset (premovement activity; Fig. 1B and fig. S1D), which is consistent with those obtained in single-cell recordings of S1 neurons during self-generating movements (8, 9, 19). The model from the combined activities accurately reconstructed premovement activity in S1. The decoding model further reconstructed the peaks of S1 activity during movement. The model from combined activities in MCx and peripheral afferents reconstructed movement-related modulation of S1 activity throughout the duration of S1 modulation better than models from MCx activity alone or peripheral afferent activity alone (−100 to 1400 ms around movement onset for Monkey T, −100 to 1500 ms for Monkey C; fig. S4, A and C). These results suggest that S1 integrates information from both MCx and peripheral afferents.

Decomposition of S1 activity

To elucidate the contributions of MCx and peripheral afferent activities to the reconstructed S1 activity, we calculated each component of the reconstructed activity from either MCx or peripheral afferent activity and respective weight values in a decoding model that was built from the combined activities (Fig. 3, A and B). Like the observed activity in S1, the MCx component rose before movement onset [Monkey T (mean, n = 16 signals): latency of activity TS1 = −75.8 ms, TMCx + Afferent = −79.9 ms, TMCx comp. = −76.3 ms; Monkey C (mean, n = 12 signals): TS1 = −76.7 ms, TMCx + Afferent = −79.8 ms, TMCx comp. = −77.1 ms]. On the other hand, the Afferent component increased only after movement onset (Monkey T: TAfferent comp. = 106.3 ms; Monkey C: TAfferent comp. = 89.5 ms). The Afferent component increased later than the observed activity in S1, the reconstructed activity from combined activities in MCx and peripheral afferents, or the MCx component (Fig. 3C). The area of the Afferent component above the baseline during the premovement period (−100 to 0 ms around movement onset) was much smaller than that of the MCx component (Fig. 3D). Furthermore, when we built a decoding model that reconstructed activity in S1 from activity in peripheral afferents alone, the model also failed to reconstruct premovement activity in S1 (fig. S4, A, B, E, and F). These results indicate that premovement activity in S1 was decoded from activity in MCx but not from those in peripheral afferents.

Fig. 3 MCx and afferent activities account for S1 activity in temporally different ways.

(A) Average modulation of the observed S1 activity in Monkey T (S1; gray), the reconstruction using combined activities in MCx and peripheral afferents (MCx + Afferent; purple), and each component in the reconstruction (MCx component, blue; Afferent component, green) aligned to movement onset (from the same session in Fig. 2C). Vertical solid line, movement onset; vertical dashed line, average time of end of movement; solid lines, mean; shading, SD. (B) Magnification of (A). Horizontal dashed line, threshold for the onset of activity; arrows, onset times. (C) Onset times of the observed S1 activity (gray), of the reconstruction using combined activities in MCx and peripheral afferents (purple), and of each component (MCx component, blue; Afferent component, green) (Monkey T: mean, n = 16 signals; Monkey C: mean, n = 12 signals: *P < 0.01). n.s., not significant. Superimposed bars, mean. (D) The areas above the baseline of the observed S1 activity (S1, gray), the reconstruction using combined activities of MCx and peripheral afferents (MCx + Afferent, purple), and each component (MCx component, blue; Afferent component, green) in 100-ms sliding time windows. Dashed line, average time of end of movement; error bar, SE.

On the other hand, activities in both MCx and peripheral afferents could account for activity in S1 during movement. Both activities contributed to decoding the second (200 to 300 ms) and third (800 to 900 ms for Monkey T, 1000 to 1100 ms for Monkey C) peaks of S1 activity (Fig. 3D). The areas of the MCx and Afferent components were above the baseline until the end of S1 modulation (1400 ms after movement onset in Monkey T, 1500 ms in Monkey C) (Fig. 3D; see Materials and Methods). In some time windows, the areas of the Afferent component were as much as those of the MCx component in Monkey T [PMCx comp. vs Afferent comp. = 0.55 (t15 = 0.61, time = 200 to 300 ms), 0.14 (t15 = −1.54, time = 300 to 400 ms), 0.74 (t15 = −0.34, time = 400 to 500 ms), Student’s paired t test; Fig. 3D]. Furthermore, the models from combined activities in MCx and peripheral afferents reconstructed S1 activity during movement (0 to 1000 ms around movement onset) better than models from activity in MCx alone or from activity in peripheral afferents alone (fig. S4D). Thus, the activities in both MCx and peripheral afferents contributed to the reconstruction of S1 activity during movement.

We next examined how the MCx and Afferent components intermingled in the reconstructed activity in individual S1 electrodes. Figure 4A shows S1 activity and the MCx and Afferent components in eight electrodes over S1 of Monkey T. The first principal component could explain 90% of the total variance of the MCx component, suggesting that the MCx component had a similar temporal profile among all S1 electrodes (Fig. 4B). On the other hand, the temporal profile of the Afferent component among the electrodes was more varied than the MCx component. In addition, the variance of the Afferent component size among different S1 electrodes was much larger than that of the MCx component size (Fig. 4C). Thus, while the MCx component was evenly distributed over the S1 electrodes, the Afferent component was distributed in spatially varied patterns. These results suggest that most of the variability in activity in different S1 electrodes after movement initiation was likely caused by differences in the Afferent component.

Fig. 4 MCx and afferent activities account for S1 activity in spatially different ways.

(A) Average modulation of the observed activity in S1 electrodes of Monkey T (gray) and respective MCx (blue) and Afferent (green) components aligned to movement onset. Vertical lines, movement onset; solid lines, mean. (B) Principal component analysis (PCA) of the MCx and Afferent components in all S1 electrodes. Cumulative variance explained by the PCA components plotted against the number of feature dimensions used. Error bars, SE. (C) Coefficient of variation (CV) of areas above the baseline of the MCx and Afferent components in all S1 electrodes throughout the duration of S1 modulation (Monkey T: n = 21 sessions; Monkey C: n = 6 sessions; *P < 0.05). A pair including an outlier (CV of the Afferent component, 5.56; CV of the MCx component, 0.37 in Monkey C) is not used in the analysis. Superimposed bars, mean.

We then asked to what extent activities in the separate regions of MCx (i.e., M1, PMd, and PMv) contributed to the reconstruction (fig. S5). Weight values given to M1 activity in the model were significantly larger than those given to PMd or PMv activity in our model (fig. S5C). In addition, the area of the M1 component above the baseline was substantially larger than that of the PMd or PMv component throughout the duration of S1 modulation [Monkey T (mean, n = 16 signals): area normalized by S1, AMCx comp. = 53.1%, AM1 comp. = 52.0%, APMd comp. = 0.3%, APMv comp. = 3.8%; Monkey C (mean, n = 12 signals): AMCx comp. = 75.1%, AM1 comp. = 60.1%, APMd comp. = 16.3%, APMv comp. = 6.4%; fig. S5D]. Furthermore, M1 also contributed more to the decoding of S1 activity during the premovement period than did PMd and PMv (Monkey T: AMCx comp. = 77.8%, AM1 comp. = 76.4%, APMd comp. = 0.3%, APMv comp. = 5.1%; Monkey C: AMCx comp. = 68.5%, AM1 comp. = 56.0%, APMd comp. = 18.2%, APMv comp. = 5.0%; figs. S5E and S6, A to C). Thus, the M1 component contributed to most of the MCx component, while the PMd or PMv component contributed little. We further examined M1 electrodes that predominantly contributed to the decoding output by analyzing a weight value given to each M1 electrode. Electrodes with the highest weight value were located just anterior to the S1 electrode whose activity was decoded [one-way analysis of variance (ANOVA), P < 0.001 for all S1 electrodes of both monkeys; fig. S6D]. Thus, activity in S1 before movement onset could not be explained by ongoing activities in the premotor cortex and peripheral afferents but could be at least partially accounted for by activity in a subset of M1 regions. These results suggest that S1 receives inputs from M1 or inputs from a common driving source with M1 before movement onset.

Encoding muscle activity by premovement activity of S1

Our data suggest that before movement onset, S1 potentially encodes the same information as M1 does. In a previous study, we showed that a linear regression model built from high-γ activity in M1 could reconstruct imminent muscle activity (20). In the present study, muscle activity began to rise immediately before movement onset and reached an initial peak at movement onset (Fig. 1B, third row). We built a model to reconstruct muscle activity from preceding activity in M1 during the premovement period (−500 to 0 ms around movement onset). Consistent with the previous report (20), the model succeeded in reconstructing the initial increase in the EMG profile (Fig. 5A, left). When we built a model from S1 activity during the premovement period to reconstruct imminent muscle activity, the model also reconstructed the initial rise in muscle activity (Fig. 5A, right). Furthermore, both models succeeded in the reconstruction of muscle activity during movement (Fig. 5C). The reconstruction accuracies from both M1 and S1 were statistically more significant than that of a model built from time-shuffled control data (Fig. 5, B and D), thus indicating that S1 activity before movement initiation encoded temporal changes in muscle activity as well as M1 activity did.

Fig. 5 Premovement activities in M1 and S1 encode muscle activity.

(A) Reconstruction of EMG profile (aligned to movement onset) of a hand muscle during the premovement period (−500 to 0 ms around movement onset) in Monkey C using the activity in M1 or S1 before movement onset (time = 0). Black and red traces show the observed and reconstructed EMG using the activity in M1 (left) or S1 (right), respectively. R, correlation coefficient between the observed and reconstructed EMG; solid lines, mean; shading, SD. (B) Mean decoding accuracy pooled across muscles (Monkey T: n = 12 muscles; Monkey C: n = 10 muscles). Plotted are correlation coefficients between the observed and reconstructed traces from the data (M1, cyan; S1, red) compared with those between the observed traces and the traces reconstructed from random shuffling of activity (Monkey T: n = 12 muscles; Monkey C: n = 10 muscles; *P < 0.01). Superimposed bars, mean. (C) Reconstruction of EMG profile throughout the premovement and movement periods (−500 to 2000 ms around movement onset). The diagrams use a similar format for (A). (D) Mean decoding accuracy of the reconstruction throughout the premovement and movement periods. The diagrams use a similar format for (B).

We next examined how premovement activity in each M1 and S1 electrode contributed to the decoding of muscle activity. Most of the input signals in M1 and S1 were selected by the SLiR algorithm for building the decoding model (fig. S7A). However, high weight values were assigned to a small number of adjacent M1 and S1 electrodes across the central sulcus in the decoding of both proximal and distal forelimb muscles (one-way ANOVA, P < 0.001 for all the muscles of both monkeys; fig. S7B). The result indicates that premovement activity in a core region encoded imminent muscle activity.

Information flow from M1 to S1

We then asked the question: When does M1 or S1 start to encode information about muscle activity? All recorded muscles exhibited an initial EMG burst at the point of movement initiation. The amplitude of the EMG burst varied from trial to trial. If activity in M1 or S1 in the premovement period predicts trial-by-trial variability of the amplitude of the EMG burst, then we could examine the point when M1 or S1 started to encode it. We first tested whether premovement activity in M1 or S1 encoded the EMG burst amplitude by building an encoding model from activity in M1 or S1 before movement initiation that accounted for the variability of the EMG burst amplitude on a trial-by-trial basis. Figure 6A shows scatterplots in which each dot represents the observed EMG burst amplitude versus that reconstructed in a single trial. We then assessed the model using two indices, the slope of the regression line and the correlation coefficient between the observed and reconstructed EMG burst amplitude. Models derived from either premovement activity in M1 or S1 predicted EMG burst amplitude in each trial [slope: Monkey T (n = 12 muscles): SM1 data = 0.16, SM1 shuffled = 0.00, SS1 data = 0.16, SS1 shuffled = 0.00; Monkey C (n = 9 muscles): SM1 data = 0.08, SM1 shuffled = 0.00, SS1 data = 0.13, SS1 shuffled = 0.00; correlation coefficient: Monkey T: RM1 data = 0.40, RM1 shuffled = 0.07, RS1 data = 0.40, RS1 shuffled = 0.06; Monkey C; RM1 data = 0.28, RM1 shuffled = 0.11, RS1 data = 0.35, RS1 shuffled = 0.09; Fig. 6, B and C, and fig. S8, A and B]. For both M1 and S1 activities, both the slope and correlation coefficient were significantly larger than those derived from trial-shuffled datasets (Fig. 6, B and C). Together, these results indicate that, similar to M1, premovement activity in S1 encodes EMG burst amplitude.

Fig. 6 Premovement activity in S1 encodes the initial burst of muscle activity slightly after M1 does.

(A) Scatterplots of observed peak EMG amplitude of the hand muscle versus reconstructed EMG amplitude using the premovement activity in M1 (cyan) or S1 (red) of Monkey T. Dot, a single trial. Equation of the fitting is shown in the lower right corner. R, correlation coefficient between the observed and reconstructed EMG. (B and C) Average slopes of regression lines (B) and correlation coefficients (C) pooled across muscles (Monkey T: n = 12 muscles; Monkey C: n = 10 muscles; *P < 0.05). Superimposed bars, mean. (D and F) Average slopes of regression lines (D) and correlation coefficients (F) were plotted against the end of 50-ms sliding windows. Solid lines, means; shading, SD; dashed black line, threshold for onset of activity; arrows, onset. (E and G) Onset times of M1 (cyan) and S1 (red) activities that encode peak EMG amplitude pooled across electrodes were less than the onset time of EMG (E, slopes; G, correlation coefficients; Monkey T: n = 12 muscles; Monkey C: n = 9 muscles; *P < 0.01). Onset time of the activity in M1 was earlier than that in S1 (Monkey T: n = 12 muscles; Monkey C: n = 9 muscles; **P < 0.01). Superimposed bar graphs, mean.

We then examined the point at which M1 or S1 started to encode EMG burst amplitude. We built models that reconstructed trial-by-trial variability of EMG burst amplitude from activity in M1 or S1 within overlapping, sliding time windows of 50 ms and calculated the slope of the regression line and the correlation coefficient as above. We plotted two indices against the end of the 50-ms sliding window (Fig. 6, D and F). For encoding models of activity in M1, both indices increased before movement initiation. Onset times of the slope of the regression line and the correlation coefficient were −175.0 and −170.1 ms, respectively, for Monkey T (n = 12 muscles) and −171.7 and −172.0 ms, respectively, for Monkey C (n = 9 muscles) (Fig. 6, E and G). For encoding models of activity in S1, both indices also increased before movement initiation (Fig. 6, D and F). Onset times of the slope of the regression line and the correlation coefficient were −146.7 and −139.2 ms, respectively, for Monkey T (n = 12 muscles) and −131.7 and −137 ms, respectively, for Monkey C (n = 9 muscles) (Fig. 6, E and G). The onset times for both M1 and S1 were much earlier than the onset time of EMG signals (Monkey T: −47.9 ms; Monkey C: −33.6 ms; fig. S1D). Thus, both M1 and S1 encoded information about subsequent EMG signals. In addition, the onset times for both M1 and S1 were earlier than those for peripheral afferents (fig. S8, C to F). Furthermore, the onset times for M1 were earlier than those for S1 based on both the slope of the regression line (paired t test, P < 0.005 for both monkeys; Fig. 6E) and the correlation coefficient (P < 0.005 for both monkeys; Fig. 6G), suggesting that M1 had the information about muscle activity earlier than S1. That activity in M1 reconstructed premovement activity in S1 suggests that the information that M1 encodes might propagate to S1 before movement initiation.

DISCUSSION

Intensive studies have posited the idea that inputs from the motor system to the somatosensory system affect the flow of somatosensory signals from the periphery to S1 during voluntary limb movement (1217, 21). However, how and what information the somatosensory system receives from the motor system have not been clearly documented. Here, we used simultaneous recordings of multiregional neural signals and decoding models to provide insights into the flow of information that S1 receives during multijoint goal-directed forelimb movements (fig. S9). Before movement initiation, S1 activity was decoded from preceding activity in MCx, not from peripheral afferent activity. Throughout the movement, S1 activity could be decoded from both MCx and afferent activities. Furthermore, S1 encoded information about the future state of the limb, which presumably has interaction with upcoming sensory feedback signal about the actual state of the muscle activity. These sensory-motor interactions in S1 form the online processing of sensory-motor integration.

In the previous studies using analyses that simply compared onset latencies of muscle activity and neuronal activity (9, 22), evidence for the temporal and spatial profiles of signal integration remained elusive. A variety of inputs from different regions might individually contribute to the generation of activity in S1. In the present study, we built a decoding model that used a linear regression algorithm to reconstruct the movement-related activity in S1 from the weighted sum of preceding activities in MCx and peripheral afferents (Fig. 2). This model reconstructed activity in S1 from both MCx and peripheral afferent activities with an accuracy (correlation coefficient) of 0.61 in Monkey T and 0.51 in Monkey C (Fig. 2), suggesting that S1 integrates signals from both inputs. Since the accuracy of the reconstruction was not 100%, other information sources are also likely to be involved in generating S1 activity. Furthermore, our analyses indicated that premovement activity in S1 was accounted for by MCx activity and not by peripheral afferent activity, while the modulation during movement was explained by both activities (Fig. 3 and figs. S4 to S6). Decomposition of reconstructed S1 activity also showed that, while the MCx component was evenly distributed over all S1 electrodes, the Afferent component was distributed to S1 electrodes to various extents (Fig. 4). By fitting multiregional neural signals to a linear model, we showed that motor- and sensory-related signals propagated separately to S1 in the temporally and spatially different patterns.

Single neuronal activity in S1 during active movement of the forelimb has a different firing pattern from that during passive movement (8, 9, 19). These neurons are characterized by tuning properties related to the direction of voluntary movement of the hand and arm (23). Unique activities that are detected only during active movement have been considered to be a putative efference copy signal. However, these studies analyzed neuronal activities during movement so that these activities contained information from both motor-related areas and peripheral afferents. To clarify what information S1 receives from MCx, it is crucial to dissociate the effect of MCx from that of peripheral afferents. The present study showed that premovement activity in S1 could be explained, at least in part, by activity in M1 and was independent of activity in peripheral afferents. This result suggests that the premovement activity is isolated from the effect from peripheral afferents. In addition, premovement activity in S1 encoded imminent muscle activity slightly after M1 did. These results suggest that S1 receives information about motor output before receiving sensory feedback signals.

One functional imaging study of blood oxygen level dependent (BOLD) signals in human subjects has examined activity in S1 during voluntary movement without proprioceptive signals with ischemic nerve block (24). Their correlation analysis showed that activity in S1 was linked to that in the premotor cortex. However, the low temporal resolution of human imaging studies cannot dissociate motor commands from sensory signals in M1. Furthermore, anatomical evidence documented that long-range cortico-cortical connections from the premotor cortex to S1 are much sparse (25). Our results using electrophysiological recording with millisecond accuracy show that activity in M1 rather than that in premotor cortices accounted for activity in S1 not only before movement onset but also during movement (figs. S5 and S6). In addition, a linear regression analysis that accommodates time direction (see Materials and Methods) is able to reveal the signal transfer from M1 to S1. Thus, it is reliable to consider that S1 receives information regarding descending motor commands from M1 rather than from the premotor cortex. This idea was also proposed by Witham et al. (26), using directed coherence analysis, in which a β-band oscillation propagated from M1 to S1 in monkeys during a finger flexion task. The directed coherence phase from M1 to S1 documented a delay of 40 ms, which was very similar to the delay (35 ms) we observed between M1 and S1 for encoding the initial EMG burst in this study (Fig. 6, E and G).

Although most of the peripheral afferents exhibited an increase in activity after movement onset, we also found a small number of peripheral afferents that exhibited an increase in neuronal firing before movement initiation (Fig. 1B). The source of these impulses might be muscle spindles elicited by γ motor commands (27). Alternatively, preparatory postural movements might evoke discharges of peripheral afferents before limb movement. However, our present decoding results indicated that these discharges have only a very small impact on the decoding of premovement S1 activity (Fig. 3).

Premovement activity in S1 encodes imminent muscle activity as effectively as premovement activity in M1 (Figs. 5 and 6 and figs. S7 and S8). On the analogy of the role of M1 in the control of movement, it could be considered that S1 might have the ability to control the limb movement directly. The barrel cortex, S1 in rodents, controls whisker movements and innervates whisker motor neurons via the spinal trigeminal nuclei (28). However, we consider that this notion is not the case with S1 in the monkey according to the following evidence. First, the thresholds of intracortical microstimulation of S1 to evoke the limb movement are much greater than those of M1 (29). Second, the corticospinal fibers from S1 project to the dorsal horn of the spinal cord but are absent in the ventral horn where motoneurons are located (30). Third, spike- or stimulus-triggered averaging of EMG activity from S1 neurons showed that the excitatory output effects on muscle activity from S1 are almost absent, and the effects are predominantly inhibitory even if there are (31). Thus, it is more natural to consider that information about motor output is used for sensory processing in S1 rather than direct control of limb movements.

Multiple studies have shown that both somatosensation and somatosensory-evoked potentials in S1 are attenuated during voluntary movement (26). This means that the effect of M1 inputs on S1 activity could be suppressive. However, our decoding analysis showed that most of the M1 component has a positive value (fig. S5, A and B), suggesting that S1 receives excitatory effects from M1. Explanations of this incongruence are that suppressive inputs from M1 would have minimal effects on S1 neurons with low spontaneous activity and/or that inputs from other cortical areas such as the prefrontal cortex to subcortical structures might be involved in the gating of somatosensory information. On the other hand, optogenetic activation of vibrissal M1 has been reported to enhance neuronal responses in the barrel cortex to mechanical stimulation (15, 16). These neuronal responses were enhanced only when M1 neurons were activated before vibrissae stimulation, not after the stimulation (15). These findings support our demonstration that, during voluntary movement, S1 activity was decoded from M1 before the arrival of inputs from peripheral afferents. Thus, M1 inputs might have an excitatory influence on neuronal responsiveness to upcoming somatosensory inputs. Together, these results suggest that the central nervous system implements online processing of somatosensory signals under the influence of motor contexts during voluntary movement.

How is information about motor output transmitted to S1? Several pathways might transmit information about motor output from M1 to S1. In the present study, M1 encoded information related to future muscle activity 35 ms earlier than S1 (Fig. 6, E and G). This time lag is longer than a delay attributed to a monosynaptic connection; hence, a direct connection from M1 is not the primary pathway conveying the information to S1. On the basis of the conduction time of the signal, the polysynaptic cortico-cortical connection is more reasonable. Pyramidal neurons in vibrissal M1 provide strong excitatory effects to pyramidal neurons in the barrel cortex via an intracortical disinhibition circuit (32). Although it has not been shown that M1-S1 connectivity in the forelimb area in primates is the same as the rodent vibrissal system, the cortico-cortical pathway is a strong candidate. Another possible polysynaptic cortico-cortical pathway from M1 to S1 is via the posterior parietal cortex. M1 sends densely cortico-cortical projections to area 5, which interconnects with areas 1 and 2 (33).

Other possible pathways that might convey information about motor output to S1 are through subcortical structures. M1 neurons project many axonal fibers to subcortical structures, including the thalamic nuclei and dorsal column nucleus (34, 35). Electrical stimulation of motor cortex evoked unit activity in short (less than 7 ms) and long latency ranges (about 20 ms) in the thalamic somatosensory relay nuclei in the rat (14). The conduction time of the longer latency response is in the same ballpark as the difference in times for M1 and S1 to encode the information about EMG activity in our results. Thus, it is quite possible that the information is relayed through the thalamic somatosensory relay nuclei. However, as the thalamic somatosensory relay nuclei are devoid of corticothalamic fibers from M1 (35), another multisynaptic route might be responsible for activating neurons in the thalamic somatosensory relay nuclei.

M1 generally inhibits neurons that project to the medial lemniscus in the main cuneate nucleus of cats and rats (13, 36). The inhibitory effect is presumably via cuneate inhibitory interneurons that receive excitatory inputs from M1 and send inhibitory outputs to projection neurons (37). However, electrical stimulation of M1 activates cuneate neurons only when the joint controlled by a cortical site in M1 topographically corresponds to the receptive field of the cuneate neurons (38). The result implies that information about the motor map of M1 might transfer to the cuneate nucleus. Thus, the pathway through the cuneate nucleus is also likely to be a pathway conveying the excitatory effects of M1 on S1 in the present study.

While the timing of activity between M1 and S1 suggests that information is conveyed through polysynaptic connections between M1 and S1, the monosynaptic connection is likely to remain disputable. If monosynaptic inputs from M1 to S1 are not strong enough to immediately activate S1 neurons, then the delay between M1 and S1 activities (Fig. 6, E and G) might be explained by the time required for M1 subthreshold inputs to change the state of S1 neurons. After this lag, S1 activity begins to reflect monosynaptic inputs from M1.

Another possible mechanism underlying the decoding of S1 from M1 is that S1 receives inputs from a common source after M1 receives them. Previous literature has shown that area 5 and the secondary somatosensory cortex send projections to both M1 and S1 (25). Since S1 activity was accounted for by the adjacent M1 region, it seems to be reasonable to consider that both M1 and S1 receive common inputs via these cortico-cortical pathways at a similar timing, which, however, is not consistent with the delay between M1 and S1 for encoding muscle activity (Fig. 6, E and G). Alternatively, a co-modulator might send information directly to M1 and the same information indirectly to S1 with a delay. Whether S1 receives information about motor output from M1 or a common driving source should be examined by the circuit manipulation experiment, which is left for future research.

In conclusion, multiregional recording of a sensory-motor closed-loop system revealed that S1 receives information about motor output before the arrival of sensory feedback signals. This result provides insight into the online processing of somatosensory information under voluntary movement.

MATERIALS AND METHODS

Experimental design

We hypothesized that S1 integrates information from both MCx and peripheral afferents because there are anatomical pathways between them. To test this hypothesis, we conducted simultaneous recordings of the activities in S1, MCx, and an ensemble of afferent neurons in two behaving monkeys. We stopped the recording when we could not detect any spiking activity from peripheral afferent recording. We used data for the analysis when the number of unit activity of peripheral afferents was more than 20 and 10 in Monkeys T and C, respectively. All quantified data were included, and no outliers were excluded other than those mentioned.

Animals

We used one adult male monkey (Monkey T: weight, 6 to 7 kg) (Macaca fuscata) and one adult female monkey (Monkey C: weight, 5 to 6 kg) (Macaca mulatta). The animals were housed individually in temperature-controlled environments on a 12-hour light/12-hour dark cycle. The experiments were approved by the experimental animal committee of the National Institute of Natural Sciences and animals were cared for and treated humanely in accordance with National Institutes of Health guidelines.

Behavioral task

Our basic methods for the behavioral task, surgery, and recording of neuronal and kinematics signals have been described previously (20, 39). Two monkeys were operantly conditioned to perform a reach-to-grasp task with the right hand (fig. S1A). Each monkey started a trial by putting its hand on a home button for a predetermined time (2 to 2.5 s). After receiving the cue, they reached for a joystick lever and pulled it to get a reward [time from movement onset: Monkey T, 342 ± 39 ms (mean ± SD, n = 3317 trials); Monkey C, 375 ± 47 ms (n = 934 trials)]. The targets were placed 25 and 20 cm apart from the starting point for Monkeys T and C, respectively. We recorded times for releasing the home button, pulling the lever, and pushing the home button. Monkey T performed the task for 24 sessions of 10 min each, in which they conducted 142.4 ± 7.4 (mean ± SD) trials per session. Among the 24 sessions, we recorded the neuronal activity in peripheral afferents, ECoG, EMG, and kinematics in 21 sessions and ECoG, EMG, and kinematics in 3 sessions. Monkey C performed the task for 40 sessions of 10 min each, in which 127.1 ± 15.9 (mean ± SD) trials were conducted per session. Among the 40 sessions, we recorded the neuronal activity in peripheral afferents, ECoG, EMG, and kinematics in 7 sessions. We recorded ECoG, EMG, and kinematics in 33 sessions but did not record times for pulling the lever in these sessions.

Monkey T made reaching movements to one of two target locations 18 cm apart from one another along the right-left axis of the body and grasped one of the objects, with three different shapes, in each session. Monkey C made reaching movements to one target location and grasped one object in each session during which we recorded the neuronal activity in peripheral afferents. In sessions where we recorded ECoG, EMG, and kinematics only, Monkey C made reaching movements to one target location (22 sessions), one of three locations (6 sessions), one of four locations (1 session), one of six locations (1 session), and variable locations (3 sessions) on the coronal plane in each trial.

Surgery

We used a mixture of xylazine (0.4 mg/kg) and ketamine (5 mg/kg) to induce anesthesia, and then isoflurane (exhaled level, 1 to 2%) and nitrous oxide gas (1 to 2%) to maintain anesthesia. During the implantation of electrode arrays into DRGs, the monkeys were paralyzed using pancuronium bromide (0.2 mg/kg per hour; Mioblock). Expiratory CO2 levels were continuously maintained within the physiological range (3.3 to 4.2%), and the depth of anesthesia was checked by monitoring expiratory CO2 levels and heart rate. Dexamethasone, ketoprofen, and ampicillin were postoperatively administered.

For EMG recording, pairs of Teflon-insulated wire electrodes (AS 631, Cooner Wire) were secured into the forelimb muscles on the right side using silk suture. The wire electrodes were implanted in the deltoideus posterior, triceps, biceps, brachioradialis, extensor carpi radialis, extensor digitorum communis, palmaris longus, flexor digitorum profundus, flexor digitorum superficialis, flexor carpi ulnaris, extensor digitorum 2 and 3, and adductor pollicis of Monkey T and the pectoralis major, deltoideus posterior, triceps brachii longus, triceps lateralis, biceps, brachioradialis, extensor carpi radialis, extensor digitorum communis, flexor digitorum profundus, flexor carpi ulnaris, abductor pollicis longus, and adductor pollicis of Monkey C. Since we observed the electrocardiogram mixed with activities in the pectoralis major and biceps in Monkey C, we did not use these data in analyses.

To implant a grid electrode array on the cortical surface, we made a craniotomy to expose the premotor cortex, M1, and S1 on the left side. We implanted a 32-channel grid electrode array, in which the diameter of each electrode was 1 mm and the interelectrode distance was 3 mm (Unique Medical), beneath the dura mater (fig. S1B). We placed the ground and reference electrodes over the ECoG electrode so that they contacted the dura. After implanting the array, we fixed a connector on the skull via dental acrylic.

To implant electrode arrays into DRGs, we bilaterally exposed the C3 through the Th2 vertebrae and inserted stainless screws into the lateral mass of each vertebra on both sides. After we dissected a lateral mass of C5-Th1 segments on the right side, we implanted two multielectrode arrays, consisting of 48 platinized-tip silicon electrodes with 0.1 to 1 megohm at 1 kHz, 1 mm in length, 400 μm in interelectrode distance, and in a 5 × 10 configuration (Blackrock Microsystems), into two cervical DRGs (Monkey T: C7 and C8; Monkey C: C6 and C7) on the right side using a high-velocity inserter. We placed reference wires over the dura of the spinal cord. After implanting the arrays, we attached a connector to the spine using dental acrylic.

Movement recordings

Forelimb movements were recorded using an optical motion capture system that used 12 infrared cameras (Eagle-4 Digital RealTime System, Motion Analysis). The spatial positions of the reflective markers (4- or 6-mm-diameter spheroids) were sampled at 200 Hz. Ten markers were attached to the surface of the forelimb using mild adhesive. Positions of the 10 markers were as follows: the left shoulder (marker 1; m1), the center of the chest (m2), the right shoulder (m3), the biceps (m4), the triceps (m5), the lateral epicondyle (m6), medial to m6 (m7), the radial styloid process (m8), the ulnar styloid process (m9), and the metacarpophalangeal joint of digit 2 (m10). We calculated flexion/extension (FE) of the shoulder, adduction/abduction (AA) of the shoulder, FE of the elbow, pronation/supination (PS), FE of the wrist, and radial/ulnar (RU) of the wrist (table S1). To reduce noise from various sources, we applied a low-pass filter with a cutoff frequency of 5 Hz to temporal changes in the joint angles.

Neural recordings and spike detection

EMG signals were amplified using amplifiers (AB-611J, Nihon Kohden) with a gain of ×1000 to 2000 and were sampled at 2000 Hz in Monkey T and 1000 Hz in Monkey C. Temporal filtering of the signals was carried out with a second-order Butterworth band-pass filter (1.5 to 60 Hz). The signals were rectified and computed in 5-ms bins corresponding to the sampling rates of the motion capture system. A smoothed curve of the signals was then calculated using a moving window process with a window length of 11 bins.

ECoG signals were amplified using a multichannel amplifier (Plexon MAP system, Plexon) with a gain of ×1000 and sampled from each electrode at 2000 Hz in Monkey T and 1000 Hz in Monkey C. Temporal filtering of the signals was carried out with a second-Butterworth band-pass filter (1.5 to 240 Hz). We computed short-time fast Fourier transform on moving 100-ms windows of the preprocessed signals. We used a 200-Hz frequency step size to match the sampling rate to that of the motion capture system. We computed power normalized to the averaged power in each session and calculated an averaged power in high-γ bands (high-γ 1, 60 to 120 Hz; high-γ 2, 120 to 180 Hz). SLiR analysis showed that high-γ power of ECoG signals in M1 encoded the kinematic variables (fig. S2). Moreover, the same analysis indicated that high-γ power of ECoG signals in S1 encoded the kinematic variables immediately before cortical activity (fig. S2). As a neuronal ensemble activity of M1 neurons encoded the kinematic variables (40), we used the high-γ power of ECoG signals as representative of neural activity in cortical areas in the analyses. Data from one electrode of Monkey C were not used for the analysis because of high noise (black circle in fig. S1B).

Neuronal signals of peripheral afferents were initially amplified using the same multichannel amplifier with a gain of ×20,000 and sampled from each electrode at 40 kHz. We extracted filtered waves (150 to 8000 Hz) above an amplitude threshold that was determined by the “auto-threshold algorithm” of the software. We sorted the thresholded waves using semiautomatic sorting methods (Offline Sorter, Plexon), followed by manual verification and correction of these clusters if needed. When the interval of the consecutive spikes was less than 1 ms, the second spikes were removed. To obtain the instantaneous firing rate, we convolved the inversion of the interspike interval with an exponential decay function whose time constant was 50 ms. We computed the firing rate in 5-ms bins, corresponding to the sampling rates of the motion capture system. When we examined the modality of recorded units, we identified some units as muscle spindles, tendon organs, and cutaneous receptors by moving the forelimb, tapping over the muscle belly, and brushing the skin (39).

Sparse linear regression

We applied a Bayesian SLiR algorithm that introduces sparse conditions for the unit/channel dimension only and not for the temporal dimension of the model. High-γ power recorded in S1 electrodes was modeled as a weighted linear combination of the neuronal activity of peripheral afferents and high-γ activity in MCx using multidimensional linear regression as followsyj,T(t)=k,lwj,k,l×xk,T(t+lδ)+bj(1)where yj,T(t) is a vector of activity of an S1 electrode j (two frequency bands of eight and six electrodes in Monkeys T and C, respectively) at time index t in a trial T, xk,T(t + lδ) is an input vector of a peripheral afferent or a cortical electrode k at time index t and time lag lδ (δ = 5 ms) in a trial T, wj,k,l is a vector of weights on a peripheral afferent or a cortical electrode k at time lag lδ, and bj is a vector of bias terms to yj,T. Because we examined how combined activities in MCx and peripheral afferents influenced activity in S1, time lag lδ (Eq. 1) was set to negative values. We used a time window of 100 ms because the prediction accuracy reached a plateau at 100 ms.

To compute the contribution of each cortical area or peripheral afferents to the reconstruction of S1 activity, we calculated each component of reconstructed activity using MCx, premotor cortices, M1, or peripheral afferent activity and their respective weight values in a decoding model that was built from combined activities in MCx and peripheral afferents. For example, the MCx component was calculated as followsy_MCxj,T(t)=k,lwj,k,l×x_MCxk,T(t+lδ)+bj(2)where y_MCxj,T(t) is a vector of the MCx component at an S1 electrode j at time index t in a trial T, x_MCxk,T(t + lδ) is an input vector of a cortical electrode k at time index t and time lag lδ in a trial T, and wj,k,l is derived from a vector of weights in Eq. 1, but with weights assigned to peripheral afferents removed.

The temporal activity of muscles was modeled as a weighted linear combination of high-γ activity in M1 or S1 using the above Eq. 1. In the analysis, yj,T(t) is a vector of EMG of a muscle j (12 and 10 muscles of Monkey T and C, respectively) at time index t in a trial T. xk,T(t + lδ) is an input vector of a channel k at time index t and time lag lδ (δ = 5 ms) in a trial T. wj,k,l is a vector of weights on a channel k at time lag lδ for a muscle j, and bj is a vector of bias terms to yj,T. As we examined how activity in M1 or S1 before movement initiation influenced the initial increase in muscle activity, time lag lδ was set to negative values. We built a model from the activity in M1 or S1 to reconstruct the subsequent muscle activity during the premovement period (−500 to 0 ms around movement onset). Then, we reconstructed EMG by applying the obtained model to the M1 or S1 activity throughout the premovement and movement periods (−500 to 2000 ms around movement onset). We used a time window of 50 ms because the prediction accuracy reached a plateau at 50 ms.

The initial peak EMG amplitude was modeled as a weighted linear combination of the high-γ activity in M1 or in S1 within an overlapping, sliding time window of 50 ms as followsyj,T=k,lwj,k,l×xk,T(t+lδ)+bj(3)where yj,T is a vector of EMG of a muscle j (12 and 10 muscles of Monkey T and C, respectively) in a trial T, xk,T(t + lδ) is an input vector of a channel k at time index t and time lag lδ in a trial T, wj,k,l is a vector of weights on a channel k at time lag lδ for a muscle j, and bj is a vector of bias terms to yj,T. As we examined how the activity in M1 or S1 influenced the initial peak amplitude of muscle activity, time lag lδ was set to negative values. We used a time window of 50 ms so that l was set to −10. To examine the point at which M1 or S1 started to encode EMG burst amplitude, we changed time index t.

Joint angles were modeled as a weighted linear combination of neuronal activities in peripheral afferents or high-γ power in M1 or S1 using multidimensional linear regression as followsyj,T(t)=k,lwj,k,l×xk,T(t+lδ)+bj(4)where yj,T(t) is a vector of kinematic variables j (joint angle) at time index t in a trial T, xk,T(t + lδ) is an input vector of unit k at time index t and time lag lδ (δ = 5 ms) in a trial T, wj,k,l is a vector of weights on a peripheral afferent or a cortical electrode k at time lag lδ, and bj is a vector of bias terms to yj,T. We considered that, in the sensory-motor closed-loop pathway, neuronal activity in M1 evokes activity in muscles, which, in turn, generate the movement of the limb. In a model of high-γ activity in M1, we set the time lag lδ (Eq. 4) to negative values. In contrast, we considered that self-generated movements evoke the neuronal activity of peripheral afferents and high-γ activity in S1; therefore, in models of peripheral afferents or high-γ activity in S1, we set the time lag lδ to positive values. By changing the length of the time window, we attained a time window (400 ms) within which the accuracy of reconstructing joint kinematics reached a plateau. We then used this time window in encoding forelimb kinematics from neural activities.

Data analysis

We built models to reconstruct activity in S1, temporal changes in EMG, or joint kinematics using a partial dataset (training dataset) and tested them using the remainder of the same dataset (test dataset). For reconstruction of S1 activity and joint kinematics, we partitioned continuously recorded data of each session into 24 blocks (one block for 25 s of data). Among the 24 blocks, 20 randomly selected blocks were used for the training dataset and the remaining 4 blocks were used for the test dataset. For reconstruction of EMG signals, five-sixths of the full trials of each session were randomly selected as a training dataset and the remainder were selected as the test dataset. To assess the model, we calculated the correlation coefficient between observed data and their reconstruction in the test dataset. We also calculated variance accounted for (VAF) as followsVAF=1(y(t)f(t))2(y(t)y(t)¯)2(5)where y(t) is a vector of the actual activity in S1 at time index t, y(t)¯ is the mean of y(t), and f(t) is the reconstructed activity at time index t. We performed sixfold cross-validation in the analysis of each session and used averaged values for the analysis. Then, we calculated averaged values of each electrode and each kinematic variable from data taken from 21 (Monkey T) and 7 (Monkey C) sessions in decoding S1 activity and joint kinematics and averaged values of each muscle from data taken from 24 (Monkey T) and 40 (Monkey C) sessions in decoding EMG activity. In control analyses of the model reconstruction, we created surrogate training datasets in which we shuffled temporal profiles of inputs independently across different blocks to generate a model and subsequently tested the model.

We built a model to predict the initial peak amplitude of EMG using a training dataset and tested it using a test dataset. Five-sixths of the full trials were randomly selected as a training dataset and the remainder were selected as the test dataset (Monkey T: 2847 trials for training, 570 trials for test; Monkey C: 4234 trials for training, 847 trials for test). To assess the model, we calculated the correlation coefficient between the observed and reconstructed EMG amplitude in the test dataset. We also calculated the largest possible variance (the first principal component) between the observed and reconstructed EMG amplitude to obtain the slope of the regression line. We drew regression lines that passed the centroid of the data in the plot (Fig. 6A). Among models of high-γ activity in M1 or S1 in sliding time windows in the premovement period (−500 to 0 ms around movement onset), we selected a model that reconstructed the initial peak amplitude of each muscle at the highest accuracy (Fig. 6, B and C). We performed sixfold cross-validation in the analysis. In control analyses of the model reconstruction, we created surrogate training datasets in which we randomized the trials of input array (trial shuffling) to generate a model and subsequently tested the model.

We built a model to predict the initial peak amplitude of EMG from M1, S1, or peripheral afferents activity using a training dataset and tested it using a test dataset. We first determined putatively the same units among different sessions according to the shape of waveforms and the distribution of interspike intervals. Thirty-one units were identified as putatively the same units among nine different sessions in Monkey T. Five-sixths of the full trials were randomly selected as a training dataset and the remainder were selected as the test dataset (1091 trials for the training, 219 trials for the test). To assess the model, we calculated the correlation coefficient and the slope of the regression line between the observed and reconstructed EMG amplitude in the test dataset. We performed sixfold cross-validation in the analysis. In control analyses of the model reconstruction, we created surrogate training datasets in which we randomized the trials of input arrays (trial shuffling) to generate a model and subsequently tested the model.

To obtain the onset time of the activity or the reconstruction, we first calculated the average of the aligned waveform in a test dataset [Monkey T: 21.4 ± 0.7 trials (mean ± SD, n = 21 sessions); Monkey C: 19.8 ± 1.0 trials (n = 7 sessions)]. Then, we calculated a threshold of the averaged aligned waveform by adding an average of activity during the baseline period (−1250 to −250 ms around movement onset) to one-fifth of the amplitude of the activity from 250 ms before to 100 ms after movement onset. If the activity was over the threshold in five consecutive bins, then the first of these bins was set as the onset of the activity. The calculated onset corresponded well with that based on visual inspection. We calculated the average values of the onset from six test datasets in one session and finally obtained their average from the whole sessions. We also used the average activity during the baseline period to calculate the area above the baseline. To obtain the end point of the activity in S1, we first calculated a threshold of the aligned waveform by adding an average of activity during the baseline period (−1250 to −250 ms around movement onset) to 3/10 of the amplitude of the activity from 250 ms before to 100 ms after movement onset. If the activity 1000 ms after movement onset was below the threshold in five consecutive bins, the first of these bins was set as the end point of the activity in S1. The end point of the activity in S1 was 1307 ± 27 ms (mean ± SD, n = 16 signals; eight electrodes, two bands) after movement onset for Monkey T and 1450 ± 51 ms (mean, n = 12 signals; six electrodes, two bands) for Monkey C.

To calculate the point when M1 or S1 started to encode EMG burst amplitude (Fig. 6), we used a common threshold in the encoding of M1 and S1. For both the correlation coefficient and the slope of the regression line, we calculated the encoding magnitude by subtracting the baseline activity during the baseline period (−845 to −255 ms around movement onset) from the magnitude from 245 ms before to 255 ms after movement onset. Then, we calculated the 1/5 (correlation coefficient) or 1/20 (slope) of the encoding magnitude and used a larger one of the values for M1 and S1 as a common threshold. Last, we added the common threshold to the baseline activity for M1 and S1 and used this value as respective thresholds. If the indices were over the threshold in three consecutive bins, then the first of these bins was set as the point when M1 or S1 started to encode EMG burst amplitude. Reconstruction accuracy (slope) of the flexor digitorum profundus muscle of Monkey C was only 0.02 in the decoding from M1 and 0.04 from S1, so we did not calculate the point when M1 or S1 started to encode burst amplitude for this muscle.

To calculate the point when M1, S1, or peripheral afferents started to encode EMG burst amplitude (fig. S8), we used a common threshold for their encoding. For both the correlation coefficient and the slope of the regression line, we calculated the encoding magnitude by subtracting the baseline activity during the baseline period (−845 to −255 ms around movement onset) from the magnitude from 245 ms before to 255 ms after the movement onset. Then, we calculated the 1/5 (correlation coefficient) or 1/20 (slope) of the encoding magnitude and used a larger one of the values for M1 and S1 as a common threshold. Last, we added the common threshold to the baseline activity for M1, S1, and peripheral afferents and used this value as respective thresholds. If the indices were over the threshold in three consecutive bins, then the first of these bins was set as the point when M1, S1, or afferents started to encode EMG burst amplitude.

To calculate the variability of the averaged profile of the MCx and Afferent components, a principal component analysis (PCA) was performed on the MCx and Afferent components throughout the premovement and movement periods (−500 to 2000 ms around movement onset) of eight (Monkey T) or six (Monkey C) electrodes. Values for high-γ 1 and 2 were averaged before PCA was conducted. More than 90% of the total variance was captured by the first principal component for the MCx component and by the first two principal components for the Afferent component.

To obtain a weight value given to each electrode, weight vector wj,k,l in Eq. 1 or 3 in the manuscript was averaged across time points. Values for high-γ 1 and 2 were averaged. The one-way ANOVA was used to determine whether there are any statistically significant differences between the means of weight values of different electrodes.

Statistical analysis

We used the nondirectional paired Student’s t test. An α level of significance was set at 0.05 for all statistical tests. Data are expressed as means ± SE or means ± SD. We used MATLAB R2015b (MathWorks) for the statistical analysis. Data distribution was assumed to be normal, but this was not formally tested. No statistical methods were used to predetermine sample sizes. However, sample sizes were estimated by methodologically comparable previous experiments in the laboratory and are similar to those employed in the field.

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/7/eaaw5388/DC1

Fig. S1. Closed-loop sensory-motor circuits were simultaneously recorded from monkeys.

Fig. S2. Peripheral afferents, M1, and S1 activities encode forelimb joint kinematics.

Fig. S3. Both MCx and peripheral afferent activities contribute to the decoding of S1 activity.

Fig. S4. MCx, not peripheral afferent, activity contributes to the decoding of premovement S1 activity.

Fig. S5. M1 activity is a better predictor of S1 activity than premotor cortex.

Fig. S6. M1 activity contributes to the decoding of premovement S1 activity.

Fig. S7. Premovement activity in a core area encodes muscle activity.

Fig. S8. Premovement activities in M1 and S1 encode EMG burst.

Fig. S9. Proposed temporal dynamics in which S1 receives information about motor output and somatosensory feedback signals.

Table S1. Calculation of the joint angles.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: We thank Y. Yamanishi for animal care and Y. Nishihara, M. Suzuki, and M. Togawa for technical help. Funding: The monkeys used in this study were provided by the National Institute of Natural Sciences through the National BioResource Project of the Ministry of Education, Culture, Sports, Science, and Technology of Japan (MEXT). This work was supported by grants from the “Brain Machine Interface Development” and performed under the Strategic Research Program for Brain Sciences from MEXT and the Japan Agency for Medical Research and Development; the Japan Science and Technology Agency; the Ministry of Education, Culture, Sports, Science, and Technology; Grant-in-Aid for Scientific Research (KAKENHI 23680061 and 25135733 to Y.N.); and Medtronic Japan External Research Institute (to T.U.). Author contributions: T.U. and Y.N. conceived and designed the experiments. T.U. and Y.N. performed the experiments. T.U. analyzed the data and prepared figures. All authors wrote the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View Abstract

Navigate This Article