Introduction

The motor cortex, a brain region generating motor commands, is widely known for its relation to movement kinetics^1,2 and kinematics^3–5. Moreover, it has been reported to also carry substantial sensory information^6–13. Activity in the motor cortex can be strongly influenced by reference frame¹⁴, hand trajectory¹⁵, and stimuli for target selection¹⁶; there even exists evoked activity during visual replay related to the tracking movements¹⁷. However, it is unclear how such sensory input contributes to pre-movement activity in the motor cortex.

From the dynamical systems perspective, preparatory population activity develops toward a movement-specific optimal subspace to set initial states seeding the motor generation^18–20. However, sensory information may change the distribution of neural states during movement in an orthogonal fashion without affecting the motor output.

To investigate the neural dynamics shaped by concurrent sensorimotor signals, we recorded population activity in the primary motor cortex (M1) from monkeys performing a flexible manual interception task. Unlike previous studies constraining interception at a fixed location^9,21, our task demands predictive spatiotemporal mappings to displace a body effector to a trial-varying location. We found that the activity of most neurons was jointly tuned to both reach direction and target motion via directional selectivity shifts, gain modulations, offset adjustments, or their combinations. Strikingly, such mixed sensorimotor selectivity exists throughout the entire trial, in contrast to the gradient of sensory-to-motor tuning from cue to movement epochs in posterior parietal cortex (PPC) that we recently reported ^22,23. Principal component analysis (PCA) of the neural population revealed a clear orbital neural geometry in low-dimensional space at movement onset: The neural states were distributed in reach-direction order, and formed ringlike structures whose slopes were determined by target-motion conditions. Neuronal simulation indicates that these characteristics of neural population dynamics could be derived from the mixed sensorimotor selectivity of single neurons. A recurrent neural network (RNN) trained with proper input-output mappings offers insights into the relationship between neuronal modulation and neural geometry in a dynamical system. We propose that this sensory modulation occurs at both single-neuron and population levels as a general element of neural computations for predictive sensorimotor control.

Results

Mixed tuning of M1 single neurons during flexible manual interception task

Three monkeys (Macaca mulatta, C, G, and D, male, weight 7 to 10 kg) were trained to perform a delayed manual interception task (Figure 1A), which was modified from that utilized in our recent studies ^22,23. To initiate a trial, the monkey needed to hold the center dot of a standing touch screen for 600 ms. Then, a peripheral target appeared, either stationary or rotating around the center dot at a constant angular speed. The monkey was required to wait during a random delay (400-800 ms) until the central dot went dark (GO) and then to immediately reach to the target. Once the monkey touched the screen (Touch), the target stopped and another dot showed the touching location, in red for a successful interception or in blue for a failure. The error tolerance between the target and the touch location (reach endpoint) was 3 cm. There were five targetmotion conditions, consisting of clockwise (CW) conditions - 240 °/s and -120 °/s, counterclockwise (CCW) conditions 120 °/s and 240 °/s, along with 0 °/s (static) condition. We define the target-speed magnitude as TS_mag (0 °/s, 120 °/s, and 240 °/s) and the target-speed direction as TS_dir (CCW and CW). These target-motion conditions were interleaved trial by trial, and initial location of the target was random. The reach endpoints of a well-trained monkey distributed uniformly around the circle (Figure 1B, Rayleigh’s test: p=0.36; data from monkey C, 772 correct trials in one session). Because the reach direction was defined as the angle of the reach endpoint, for simplicity, we divided the circular space equally into eight sectors (45° per each), and grouped trials according to the eight reach directions and five target-motion conditions.

Figure 1.

We recorded neural data with Utah arrays from monkeys C, G, and D (implanted sites are shown in Figure 1C, and all datasets listed in Table S1) and hand trajectories from monkeys C and G (Figure S1), during the interception task. The hand trajectory was launched to the final interception position. The temporal profiles of hand velocity were unimodal bell-like curves and similar across target-motion conditions (correlation coefficients is 0.97±0.05 for monkey C and 0.99±0.02 for monkey G).

Notably, we found that neuronal directional tuning during motor execution (MO±100 ms, MO for movement onset) was modulated by target motion, according to our statistical criteria mainly in three ways: preferred direction (PD) shift, gain, and offset addition (Figure 1D and Figure S2-S4, Methods). PD-shift neurons had their PDs shifted in moving-target conditions compared to the static-target condition. As illustrated by an example neuron (Figure 1D, PD shift), whose PDs corresponding to CCW conditions (red and yellow) and CW conditions (blue and purple) exhibited obvious differences, TS_dir rather than TS_mag dominated this modulation. Gain-modulation neurons exhibited reach-direction tuning multiplied by target speed: while their directionality remained invariant, the neuronal responses at PD differed across target-motion conditions. This modulation was dominated by TS_dir as well. The turning curves of the example neuron (Figure 1D, Gain) had higher responses at PD in CW conditions (blue and purple) than in the others (green, yellow, and red), indicating a varying tuning depth for reach direction. Neurons with addition modulation underwent changes of offset activity induced by target speed (both TS_dir and TS_mag). As shown by the example neuron (Figure 1D, Addition), this effect was almost equal for all reach directions.

The activity of these example neurons could be well explained by PD shift, gain, and additive models (Figure 1E, Methods), respectively. Nevertheless, it was difficult to classify all neurons with mixed sensorimotor selectivity into one of these three groups exclusively, because many of them experienced a mixture of two or three of above modulations. Moreover, the proportions of three groups were almost equal (Table 1). We found that the adjusted R₂ of a full model (0.64 ± 0.21, mean ± sd.) was larger than that of the PD shift (0.59 ± 0.23), gain (0.56 ± 0.21), and additive (0.47 ± 0.24) models for monkey C (95 neurons), as well as for all three monkeys (rank-sum test, p<0.01, Figure S5). These target-motion modulations of neuronal directionality suggested the participation of sensory signals in shaping neural dynamics during interception execution.

Table 1.

Single neurons were classified by comparing firing rate patterns between target-speed conditions (more details in Methods). The first column is the mean±std. of recording unit number in one monkey across sessions. The other columns give the ratio of target-speed modulation unit across sessions.

Coding of sensory and motor information in neural populations

To quantify target motion information embodied in neural response in motor cortex, we performed a series of decoding analyses on the neural data from monkey C (n=95, 772 correct trials). To begin with, neural data were utilized to train a support vector machine (SVM) classifier for targetmotion conditions (one in five), and another for reach directions (one in eight; Methods). As Figure 2A shows, the decoding accuracy of the target-motion condition increased quickly and peaked at over 70% around GO, while the decoding accuracy of reach direction climbed and reached a plateau of about 80% before MO. These results are similar to those of the other two monkeys and the pseudo-populations (Figure S6), suggesting that target motion information existed in M1 throughout almost the entire trial. According to the demixed PCA results, while reach direction occupied top components, there was interaction components that should not be neglected and even explained more variance than target speed components (Figure S7)²⁴.

Figure 2.

Then, another decoding analysis was implemented to probe the potential interaction between reach-direction and target speed during execution period. We trained a reach-direction decoder (one in eight) to check if the decoder of one certain target-motion condition could be transferred to other conditions (Figure 2B left). It turned out that the performance of the transferred decoder deteriorated more significantly for TS_dir. (CCW vs. CW, mean ± sd. of accuracy 0.26±0.06), compared with that for TS_mag (120 vs. 240, 0.50±0.06, paired t-test, p<0.01) and for target’s state (static vs. motion 0.55±0.06, paired t-test, p<0.01). This result suggests that the coding of reach direction was rather sensitive to TS_dir, but containing similarities for static and moving targets. As control, we also compared the neural coding rules across different reach-direction conditions. We trained a target-speed decoder (one in five), and similarly checked the transferred decoding accuracy (Figure 2B right). We observed that the target-speed decoder was locked with reach direction, as the transferred decoding accuracy diminished with increasing difference of reach direction. These results qualitatively imply the interaction as that target speeds affected the reach-directional tuning, especially by TS_dir..

To explore how sensory information influences neural dynamics while preserving motor output, we performed PCA on the normalized population activity. We obtained the trial-averaged neural trajectories (5 target speeds x 8 reach directions 40 conditions) after TO or after GO (Figure S8A). The distance between neural trajectories grouped by reach-direction conditions was larger than the distance grouped by target-speed conditions especially after GO (Figure S8B), consistent with the demixed PCA results that M1 primarily encoded reach direction rather than target speed. In addition, the distance of neural trajectories between CCW and CW was much larger than the distance between 120°/s and 240 °/s conditions, indicating that TS_dir dominated the target-speed effect, agreeing with the decoding results (Figure S8C).

However, these neural trajectories were not yet the ideal description, because they were shaped mostly by time. So, to highlight the proposed target-speed effect on reach direction, we focused on four key time windows and snapshotted the neural trajectory as neural state to extract the coding rule at single trial level and in a geometric view. Here, we define the “neural states” as the projection of single-trial data during a specific time bin on principal components^25–27. At MO, the first two principal components (PCs) of the neural states explained the most variance ([24.8%, 13.8%]) and were most related to reach direction (the goodness of fitting reach direction, , Figure 2C). While reach direction was represented by the first two PCs at GO and during movement execution, target speed influenced the tuning pattern of the first three PCs (Figure 2D), suggesting that target motion affects M1 neural dynamics via a topologically invariant transformation.

Orbital neural geometry in latent space

We visualized the neural states in the low-dimensional space spanned by the above three PCs at MO. The projections of the single-trial neural state onto the PC1-PC2 subspace clustered, or rather distributed, in reach-direction order (Figure 3A top and Figure 3B left). Interestingly, the neural states formed ringlike structures under single target-motion conditions and the fitted ellipses exhibited concentric shapes (the fitting goodness of ellipses, R₂ = 0.92 ± 0.01, mean ± sd., see Methods). Moreover, these ellipses tilted in condition-dependent angles, which is particularly evident in the PC2-PC3 subspace (Figure 3A bottom and Figure 3B right). Note that here we performed an linear transformation on all resulting neural state points to make the ellipse of the static condition orthogonal to the z-axis for better visualization.

Figure 3.

Next, we quantified the spatial features of these ellipses by calculating the tilting angles, which were between the normal vectors of the moving-target and static-target conditions. Strikingly, this tilting angle was linearly correlated with target speed (both TS_mag and TS_dir), and the relationship was robust in nine datasets from three monkeys (Figure 3C and Figure S9).

Given these results, we propose that this orbital neural geometry epitomizes the sensorimotor dynamics of M1 at the population level. The sensory input can regularly modulate neural states in an orthogonal dimension (PC3), without interfering with motor generation (in PC1 and PC2).

Population neural geometry relies on neuronal tuning

We then asked whether a group of single neurons with a certain type of mixed sensorimotor selectivity could exhibit the orbit neural geometry.

For this purpose, three neuronal models were constructed for the three typical modulations described above (PD shift, gain, and addition, see Figure 1D and Figure S2-S4). We ran a simulation with these representational neuronal models (Figure 4A). Here, each group consisted of 300 model neurons with their PDs uniformly distributed, and was solely modulated as PD-shift, gain, or addition (Methods). The resulting neural geometry of the three simulation groups showed distinct features (Figure 4B): The single-condition ellipses were inclined with target-motion-dependent angles in the PD-shift and gain groups, similar to the real neural data, but the ellipses in the addition group were layered in parallel. The reach-direction clusters in the first two PCs were conservative in the gain and addition group, but not in the PD-shift group. These results indicate that the neural states of the real data mainly exhibit the geometry feature of the gain modulation group. Nonetheless, we found that a population with a uniform mixture of all three modulations was able to simulate the neural geometry as well (Figure 4C).

Figure 4.

Comparing with the static-target condition, we calculated relative rotation angle and tilting angle between ellipses along with the vertical shift of neural states, in order to quantify the simulated structure (Figure 4D). The results show that the real data yielded a smaller rotation angle than the PD-shift group, a smaller vertical shift than the additive group, but larger tilting angles than all models. The mixed group had the most similar tilting angle, though with moderate performance in rotation angle and state shift.

These simulations suggest that the existence of PD-shift and additive modulation would not disrupt the neural geometry that is primarily driven by gain modulation; rather it is possible that these three modulations support each other in a mixed population.

Recurrent neural network provides dynamic insights

To infer how such modulated subpopulations would interact with each other in a dynamical system, we trained 100 RNN models with random weight initialization (Figure 5A; see Methods). The inputs included motor intention, target location, and a Go signal. Motor intention was defined as an abstract motor command predicted to compensate for sensorimotor delays²⁸, and could be provided by the PPC^29,30, here simplified as the interception location. The network was to generate hand velocity after MO. For a fixed validation set of 500 trials, these trained network models performed well (mean square error of the reach endpoint and the target was 0.0046 ± 0.0027, a. u., mean ± sd., while the radius of target motion was 0.15).

Figure 5.

In these network models, we found comparable features of nodes. First, most activated nodes (for example, Figure 5B) could be classified into the above-mentioned three modulations (Table 1) by the same statistical standard as for the real neural data. Second, the states reduced from node population activity were arranged in a way resembling the actual neural geometry at MO (the fitting goodness of ellipses, R₂ = 0.98 ± 0.05, mean ± sd.; Figure 5C). The tilting angles followed the same pattern as suggested by the actual results (Figure 5D).

With these RNN models, we perturbed the connection between different modulation groups. The connections within and between gain modulation and additive modulation were stronger than the others (Figure 5E), however, further strengthening the connection from gain to additive nodes or within additive nodes would disorder the orbital structure more significantly (Figure S10A and Table S2). We speculate that the influence, and perhaps the formation, of modulations both relate to the robust input and output weight pattern across models (Figure S10B, C). These results support strong and similar roles of gain and additive nodes, but what is even more important is that the three modulations interact each other, so the PD-shift nodes should not be neglected.

Discussion

To reveal how sensory context influence neural dynamics during movement, we recorded M1 population activity from monkeys performing a flexible manual interception that is highly dependent on predictive sensorimotor transformations. Single-neuron activity showed that the movement tuning of M1 neurons varied with target-motion conditions in complicated ways, including PD shift, gain modulation, addition, or their mixtures. Unsupervised dimensionality reduction (PCA) on population activity revealed an orbital neural geometry, where neural states of single trials spontaneously formed ellipses and were tilted according to target-motion conditions. Such a geometry, which could be simulated with a group of representational neurons for gain modulation alone, also emerged in the RNN with appropriate input-output mappings. As suggested by a dynamic and connected population as the RNN, the interactions between modulations were sophisticated in preserving the encoding of motor output. These results reveal sensory modulation of pre-movement neural dynamics at the single-neuron and population levels, and bridge the neuronal mixed sensorimotor selectivity and a low-dimensional neural geometry, compactly representing the sensorimotor co-function.

Previous studies have showed that M1 activity is seldom evoked by purely visual stimuli, with only weak tuning^8,31,32, and has transient responses to behaviorally relevant visual cues^7,10,33–35. In this study, the modulated neural dynamics, though distinguished according to target-motion conditions, is ultimately shaped for motor generation rather than sensory representation. It is a pity that we did not record sufficient muscle data during interception to rule out the representation of the motor cortex for kinetic variables such as force^1,27,36 , but we believe that this would not invalidate our main conclusions.

Previous studies have revealed interactions between motor kinematic variables such as reach direction and hand velocity^15,37,38. Although the effect of arm-dependent speed (hand velocity) on reach-direction tuning may exist, it cannot fully explain neural activity during interception. First, we found hand velocity to vary with TS_mag, but the largest gap in neuronal activity was related to TS_dir. Second, for a new dataset with trials resampled to ensure similar distributions of hand velocities in five target-motion conditions, the orbital neural geometry remained similar and the target-motion gain model provided a better explanation compared to the hand-velocity gain model (Figure S11).

In the present study, we focused on the orbital neural geometry by taking ‘snapshots’ of neural activities with PCA. In addition to unsupervised PCA, we tried supervised dimensionality reduction methods like dPCA²⁴. However, the dPCA results (Figure S7) showed that the condition-independent temporal component explained the largest variance, followed by the component of reach direction, whereas the interactive and target-speed components only covered very small proportions. The dPCA did not distinguish a target-speed dimension, probably because reach directions were continually distributed around a circle and target speed was primarily mixed with reach directionality.

Our simulations with neuronal models and RNN demonstrate that differences in modulations and the interaction between them can be essential. In fact, modulations in the form of PD shift and gain have been widely found in the motor cortex and PPC regions contributing to sensorimotor transformation. For example, motor cortex neurons experience gain and PD shift modulation by arm posture^36,39,40. Neurons in PPC have eye and hand gain fields for a visually-guided reach plan^41–43, integrating target position related to gaze and hand position to form motor intention^29,44,45. Therefore, we conclude that three types of modulations should all be involved in sensorimotor computation in predictive motor control. In addition, recent studies show that interactions of PPC-motor cortex circuity are involved in motor planning, spatial transformation and motor selection^46–48. The PPC-M1 circuit, as a key part of cortico-subcortical networks for the predictive sensorimotor control⁴⁹, is a topic for future studies.

STAR Methods

Key resources table

Resource availability Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Dr. He Cui (hecui@cibr.ac.cn).

Materials availability

This study did not generate new unique reagents

Data and code availability

The example datasets and code have been deposited in Mendeley Preview, and will be made public as soon as the article is officially published.

Experimental model and study participant details

Monkey

Three adult male rhesus macaques (monkey C, D, and G, Macaca mulatta, 7-10 kg) were used in this study. The monkeys sat in a primate chair to perform the task. The stimuli were back projected onto a vertical touch screen (Elo Touch system; sampling at 100 Hz, spatial resolution <0.1 mm) about 30 cm in front of the monkeys. Hand trajectory was tracked by an optical camera (VICON Inc.) with an infrared marker on the fingertip. All animal maintenance and procedures were in accordance with NIH guidelines and approved by the Institutional Animal Care and Use Committee (IACUC) of Institute of Neuroscience, CAS.

Method details

Task and behavior

The monkeys were trained to perform a flexible manual interception task in a dark room. The task paradigm was modified based on the visually-guided reaching interception task in a previous study ²². At the beginning, the monkey held the green center dot of a touch screen for 600 ms to initiate a trial (Figure 1A). Then, a green target dot appeared at a random location, moving in a circular path centered on the center dot. The center dot turned dark as a GO cue after a random delay (400-800 ms); then the monkey could intercept the target at any moment within 150-800 ms after the GO cue. Once any peripheral location was touched, the target stopped. The tolerance range of the touch endpoint for correct trials was within 3 cm of the target. The monkey would be rewarded with juice after each correct trial. Targets moving clockwise (CW; -240 °/s, -120 °/s) and counterclockwise (CCW; 120 °/s, 240 °/s), as well as a stationary target (0°/s), were pseudo-randomly interleaved trial by trial. Additional target speeds (-360 °/s, -180 °/s, 180 °/s, 360 °/s added) were introduced in subsequent sessions (Table S1).

Data collection

After the monkeys were adequately trained for the task (successful rate > 90%), head-posts were implanted stereotaxically under anesthesia (induced by 10 mg/kg ketamine, then sustained by 2% Isoflurane). After few weeks of recovery and adaptation, the monkeys were implanted with Utah microelectrode arrays (Blackrock Microsystems, Salt Lake City, UT) in the motor cortex of the hemisphere contralateral to the handedness (Figure 1C, 128-channel array for monkey C, 96-channel array for monkey G and D). The location of recording areas were identified by Magnetic Resonance Imaging (MRI) and cortical sulcus features. Neuronal activity was recorded via a Blackrock Microsystem 256-channel recording system, sampled at 30 kHz. We collected 85 ± 16, 45 ± 8, 98 ± 18 well-isolated units from monkey C, G, and D across sessions (mean ± sd. more details in Table S1), respectively.

Peri-stimulus time histograms (PSTHs)

The PSTHs and spike rasters of single neurons are shown in Figure 1D and Figure S2-S4. All trials were classified into forty conditions, eight reach-endpoint sectors by five target-motion conditions. Condition-averaged firing rates were calculated with 50-ms bins and smoothed with a Gaussian kernel (standard deviation = 20 ms). The standard error of firing rates was estimated from the 10 bootstrap samples in the trials of corresponding condition.

Quantification and statistical analysis

Classification of neuronal tuning properties

To classify target-motion modulation for single neuron reach-direction tuning, we used the direction-averaged firing rate around movement onset (MO ± 100ms) to calculate three indexes for each target-motion conditions: the preferred direction (PD), the tuning depth, and the offset activity, and classified PD shift, gain, and addition groups in Table 1. The PD of each neuron for different target-motion conditions was calculated by the weighted sum of averaged neuronal firing rates in eight reach-direction sectors. The tuning depth of each neuron was determined by the range (max - min) of firing rates in corresponding target-motion conditions. The offset activity of each neuron was calculated by the mean firing rate (MO ± 100 ms) of all target-motion conditions.

Based on the PD, tuning depth and offset activity of single neurons, a neuron was classified as ‘PD shift’, if its PDs were significantly different between interception and static-target condition (Watson-Williams test in circular data, CircStat by⁵⁰); a neuron was classified as ‘gain’ if its tuning depths were significantly different between interception and static-target condition (two-tailed Wilcoxon signed-rank test, p<0.05); a neuron was classified as ‘offset’ if its offset activities were significantly different between interception and static-target condition (two-tailed Wilcoxon signed-rank test, p<0.05).

Population decoding

The population activity of the motor cortex was used to decode target motion and reach direction by support vector machine (SVM). Neuronal firing rate was soft-normalized as

where raw firing rates was divided by the range of firing rates plus five ¹⁸. We trained two SVM classifiers (MATLAB function ‘fitcecoc’, 10-fold cross-validation) to decode reach direction (one in eight) and target motion (one in five) of single trials in 100-ms window sliding with 50-ms step (Figure 2A). The temporal decoding was repeated ten times to obtain the mean and standard deviation of decoding accuracy.

We tested the generalization of reach-direction and target-motion decoders (SVM, MATLAB function ‘fitcecoc’) in different conditions during execution period (MO ± 100 ms, Figure 2B). The decoder predicted single-trial reach direction (45°, 1/8 chance level) or target motion (1/5 chance level) across conditions based on normalized population activity. The reach-direction decoder (Figure 2B left), which was trained by trials in a set of target-motion conditions, was tested with trials from other target-motion conditions (CCW vs. CW, or 120 vs. 240, or static vs. motion, randomly selected without replacement 100 trials from corresponding training-test datasets), and the training-test decoding was repeated 1000 times to compare the distribution of accuracy with paired t-test. The target-speed decoder (Figure 2B right), which was trained by trials in a given reach-direction condition within 45°, was tested with trials from other reach-direction conditions.

Neural state

The population activity was stored in NKT datasets, where N, K, and T denote the number of neurons, trials, and time bins, respectively. We averaged T dimension of neural activity in a 100-ms bin (e.g., the two 50-ms bins around MO), and normalized neuronal dimension by Z-score (MATLAB function ‘zscore’) to get a K*N dataset. After preprocessing, we used PCA to reduce the dimension from K*N to K*C (C is the number of principal components, first two and first three PCs), and fitted the PCs with reach direction (PC = a_1*cos(θ)+a_2*sin(θ)+c) and target speed in Figure 2. Neural states of single trials were colored for target speed or reach direction and fitting ellipses with first two PCs (Figure 2D) and first three PCs (Figure 3) by MATLAB package ‘MatGeom’⁵¹. To show the relative position of ellipses in the best viewing angle, we used an isometric affine transformation to globally map all datapoints on new axes while preserving the proportional relation of distances between points. After this linear transformation, the azimuth and elevation of ellipses changed slightly, but the tilting angle between ellipses remained constant (Figure 3). The tilting angles, rotation angles, and state shift were calculated between ellipses of target-motion conditions and static-target condition in each session, defining CW as negative angles and CCW as positive angles. A set of tilting angles were obtained from corresponding conditions in one dataset, and the linear regression model was used to fit all ellipses angles (θ) and target speeds (sp.) in nine sessions.

Fitting and simulating single-neuron activity

We used PD shift, gain, offset, and full models to fit neuronal activity, basing on a cosine tuning model³. Neurons are fitted by single-trial data. We introduced a special sigmoid function to fit the nonlinear target-speed effects because target-speed direction (CCW vs. CW) has a stronger effect than target-speed magnitude (120 vs. 240)^52,53. The gain and additive models refer to hand-velocity gain studies^37,38,54.

In the gain model, the target-speed effects on the amplitude of cosine tuning are denoted as:

where FR is the firing rate at movement onset (MO ± 100 ms). θ and sp. are the reach direction and target speeds, respectively; θ_pd is the fitted preferred direction of the neuron; a₁, a₂, c₁ are constants to be fitted.

In the additive model, the target speed adjusts the offset activity, as:

with similar symbols to the gain models, plus the new constant a₃.

In the PD shift model, the target-speed effects on PDs are represented as:

with the similar symbols to the above model.

The full model integrates all the three of the above effects:

with constants a₁, a₂, a₃, a₄, a₅.

We fitted neuronal activity with these four models (MATLAB ‘fit’ function), and compared the fitting goodness with adjusted R-squares , where n is the trial number, and p is the degree of the polynomial).

Simulations with model neurons were based on three models to investigate the relationship between neuronal tuning (Figure 4A) and population neural geometry (Figure 4B). We first built three model neuron groups (each n=300) and performed PCA to get the neural state. Then, we repeated this in a mixed group, with 100 each of PD shift, gain and additive model neurons (Figure 4C). Model neurons had three components: cosine-tuning for reach direction (θ_pd,n, PD uniform distribution around circle for each neuron n), Gaussian temporal profiles (t=1:200, σ=30, peak time t_µ,n ∈ N(100, 10), the 100-th bin is MO, random distribution for t_µ,n of neurons), and distinct target-motion modulation for each group. We designed five target speed values (sp. = [1,2,3,4,5]) and 64 reach directions , for 320 trials in total⁵⁵.

The gain-model neuron has a nonlinear function with target-speed effects on the amplitude:

where g_n is a random target-motion gain within [0:1] for different neurons, and modulates targetmotion gain in a sigmoidal function.

The PD-shift model neuron has target-speed change in PD:

where s_n is a random value within [0:1.5] for different neurons, and contributes to a sigmoidal function with target-motion shift on neuronal PD (θ_pd,n).

The additive model neuron has target-speed effects in offset:

where g_n is also a random value within [0:1] for different neurons, and adjusts the offset activity in a sigmoidal function with target-motion.

The firing rates of three groups model neuron (K*N*T, 320*300*200) were averaged at MO (mean T = [50:150]) to get the K*N (320*200) dataset for PCA. As with the real neural data, we selected the first three PCs (K*C, 320*3) to derive the simulated neural state shown in Figure 4.

RNN training

For the inputs of RNN, motor intention appeared from MO-50 ms to MO, and was represented as constant variables in the form of two-dimensional Cartesian coordinates; target location was designed as time varying two-dimensional Cartesian coordinates of the target throughout the entire trial; the Go-signal was a step function jumping from 0 to 1 at GO. The RNN consisted of 200 hidden units, and output hand velocity for accurate interception after MO.

The RNN nodes evolved according to a standard dynamic differential equation:

where τ is a time constant (here 50 ms), x is the activity, r is the firing rates, and u means the combined inputs. The connection matrix J of hidden layer is initiated as random in a normal distribution (mean = 0, , g = 1.5, N = 200), and the matrix B denotes the connection between inputs and hidden units, initiated as zero matrix. The output z is obtained by

where W is the read-out weight, and is expected to reproduce the desired hand velocity generated by bell-shaped physical equation⁵⁶. W was also initiated as zero matrix. During training, the loss function was:

where e is the mean squared error of z and training target. r₁ is regularity ⁵⁷, denoting the magnitude of the nodes’ activity and is calculated as activity squared summed across time. For constants, α = 1e − 7. The optimization was realized with optim.Adam() based on PyTorch, the learning rate was 0.001.

The classification of nodes’ modulation type was finished with the same procedure for real data as above. When perturbating the connection between certain types of modulations, the selected connection was multiplied with 1.5. This perturbation was retained for the entire trial. The perturbed network was tested with the same test set of 500 trials as for the intact network.

Acknowledgements

We thank C. Li, J. Malpeli, C. Zheng, and R. Zheng for helpful comments and discussions; C. Guan for veterinary assistance; and P. Ding, L. Du, and Z. Xiao for administrative support. This work was supported by the National Key R&D Program (Grants 2020YFB1313400 and 2017YFA0701102), National Science Foundation of China (Grants 31871047 and 31671075), and CAS Strategic Priority Research Program (Grant XDB32040103).

Author contributions

Y. Zhang and H. Cui designed the experiment, Y. Zhang and T. Wang collected the data, Y. Zhang and Y. Chen analyzed the data, Y. Chen built RNNs, Y. Zhang performed the simulation, Y. Zhang, Y. Chen, T Wang, and H. Cui prepared the manuscript.

Declaration of interests

The authors declare no conflict of interests.