Preserved neural dynamics across animals performing similar behaviour (2024)

Subjects and behavioural tasks

Monkeys

We trained four monkeys (monkeys C, M and J: male, Macaca mulatta; monkey T: male, M. fascicularis; aged 6–10 years) to sit in a primate chair and make reaching movements using a customized planar manipulandum. The movement of a cursor on a computer screen was mapped to the motion of the handle of the manipulandum and the behavioural task was run through custom software in Matlab (The Mathworks). Monkeys C, M and J were trained to perform a two-dimensional centre-out reaching task for at least several months before the neural recordings, ensuring they had reached expert performance. Monkeys C, M and T were trained on a more complex random target sequential reaching task. In the centre-out task, the monkey moved its hand to the centre of the workspace to begin each trial. After a variable waiting period, the monkey was presented with one of eight outer targets. The targets were equally spaced in a circle and selected randomly with uniform probability. Then, an auditory go cue signalled the animals to reach to the target. Monkeys were required to reach the target within 1 s after the go cue and hold for 0.5 s to receive a liquid reward, except for monkey J, who was trained without the instructed-delay period or the 0.5 s target hold time and therefore made larger movements (Extended Data Fig. 1a, right). For the centre-out task, there were 12 sessions for monkey C, 6 sessions for monkey M and 3 sessions for monkey J.

In the random target task, the monkeys made four consecutive reaches to random targets within a 10 × 10 cm² workspace in each trial. Each target was presented sequentially in a random location within an annulus with 5 cm inner radius and 15 cm outer radius of the previous target to enforce minimum and maximum reach lengths. Monkeys received a liquid reward during a short break after each successful sequence of four random target acquisitions. There was no explicit auditory go cue and only a brief hold period within the target (100 ms) and then a brief delay period (100 ms) before the next target was presented. These short constraints helped to enforce that the monkeys made separate, directed movements but did not require that the monkey necessarily stop between movements. For the random target task, there was one ‘reference’ session for monkey C, six sessions for monkey M and four sessions for monkey T. As the monkeys performed these tasks, we recorded the position of the endpoint at a sampling frequency of 1 kHz using encoders in the joints and digitally logged the timing of task events, such as the go cue. Portions of the centre-out reaching data have been previously published and analysed in refs. ^26,28,46,65. Portions of the random target data have been previously published and analysed in refs. ^31,32.

Mice

Four 8–16-week-old mice were trained to perform a forelimb reaching and pulling task (similar to refs. ^38,66) for approximately one month, following habituation to head-fixation and the recording setup. In each trial, mice had to reach and pull a joystick positioned about 1.5 cm away from the initial hand position. The joystick appeared, without any cue, in one of two positions (left or right, less than 1 cm apart). Mice could then self-initiate a reach to the joystick and pull it inwards to get a liquid reward. The joystick was weighted with either a 3 or a 6 g load (light or heavy), making up a total of four trial types (two joystick positions by two loads). Each trial type was repeated 20 times before task parameters were switched to the next trial type without any cue. Each session consisted of two repetitions of each set of four trial types presented in the same order, making up 2 × 4 × 20 = 160 trials. Trials with incorrect responses (for example, pushing the joystick past a threshold, 5 mm) or timeout (the lack of pull or push for 10 s) were marked as unsuccessful. All joystick operations were programmatically controlled using a custom-written open-source Python package: (https://github.com/janelia-pypi/mouse_joystick_interface_python). Mice were maintained on a 12/12 h (08:00–20:00) light/dark cycle and recordings were made between 09:00 and 15:00. The holding room temperature was maintained at 21 ± 1 °C with a relative humidity of 30–70%.

There were two sessions for mouse 38, one session for mouse 39, two sessions for mouse 40 and one session for mouse 44. Movement kinematics were tracked using markerless video-based pose estimation. Annotation of behaviour was accomplished using Janelia Automatic Animal Behavior Annotator⁶⁷. Briefly, behaviour was recorded using two synchronized high-speed (500 frames s⁻¹), high-resolution monochrome cameras (Point Grey Flea3; 1.3 MP Mono USB3 Vision VITA 1300; Point Grey Research) with 6–15 mm (f/1.4) lenses (C-Mount), placed perpendicularly in front and to the right of the animal. A custom-made near-infrared light-emitting diode light source was mounted on each camera. Video was recorded using custom-made software developed by the Janelia Research Campus Scientific Computing Department and IO Rodeo. This software controlled and synchronized all facets of the experiment. For the main analyses, light and heavy trials were pooled together because we focused on the reaching phase of the task and the location of the joystick does not depend on its weight. Note that in Extended Data Fig. 6a we repeated the main analysis to demonstrate preserved latent dynamics during the pulling phase, considering all four conditions.

Neural recordings

Monkeys

All surgical and experimental procedures were approved by the Institutional Animal Care and Use Committee of Northwestern University under protocol no. IS00000367. We implanted 96-channel Utah electrode arrays in the primary motor cortex (M1) or dorsal premotor cortex (PMd) using standard surgical procedures. Throughout the paper, neural recordings from these two subregions were pooled together and denoted as motor cortex. This allowed us to ensure that we could evaluate overt and covert dynamics within the same population. Implants were done in the opposite hemisphere of the hand animals used in the task. Monkeys M and T received two arrays in M1 and PMd simultaneously. Monkey J received a single array in M1. Monkey C received two sets of implants: one array in the right M1 while performing the task using the left hand and, following removal of this original implant, two arrays simultaneously in the left M1 and PMd while using the right hand (respectively, monkeys C_R and C_L in our previous work²⁶). Note that for all across-individual analyses, C_R and C_L are considered the same animal.

Neural activity was recorded during the behaviour using a Cerebus system (Blackrock Microsystems). The recordings on each channel were band-pass filtered (250 Hz–5 kHz) and then converted to spike times on the basis of threshold crossings. The threshold was set to 5.5× the root-mean-square activity on each channel. We also manually spike sorted the recordings from monkeys C, Mand T to identify putative single neurons. Monkey J had fewer well-isolated single units than the other monkeys, so rather than spike sorting we directly applied the multi-unit threshold crossings acquired on each electrode. However, it has been shown that the latent dynamics estimated from multi-unit and single neuron activity are similar⁶⁸, an observation that holds true for aligning latent dynamics with CCA²⁶ (note that we refer to both single neurons and multi-units simply as units). We included multiple experimental sessions from each monkey: for the centre-out reaching task, eight from monkey C_L, four from monkey C_R, six from monkey M and three from monkey J(example data in Extended Data Fig. 1); for the random target task, one ‘reference session’ from monkey C, six from monkey M and four from monkey T (example data in Extended Data Fig. 8). These experimental sessions were chosen on the basis of the high number of units or trials and blind to the behaviour of the animal. For the centre-out reaching task, the average number of units included for each monkey was: monkey C_L, 277 ± 14 (mean ± s.e.m.; range, 210–345); monkey C_R, 85 ± 4 (range, 73–92); monkey M, 117 ± 4 (range, 106–130); and monkey J, 63 ± 9 (range, 54–81). For the random target task, the average number of units included was: monkey C_L, 280 (one session only); monkey M, 127 ± 9 (range, 101–153); and monkey T, 49 ± 8 (range, 30–66). A more detailed description of the behavioural and neural recording methods is presented in ref. ²⁶.

Mice

All surgical and experimental procedures were approved by the Institutional Animal Care and Use Committee of Janelia Research Campus. A brief (less than 2 h) surgery was first performed to implant a three-dimensional-printed headplate⁶⁹. Following recovery, the water consumption of the mice was restricted to 1.2 ml per day, to train them in the behavioural task. Following training, a small craniotomy for acute recording was made at 0.5 mm anterior and 1.7 mm lateral relative to bregma in the left hemisphere. A neuropixels probe was centred above the craniotomy and lowered with a 10° angle from the axis perpendicular to the skull surface at a speed of 0.2 mm min⁻¹. The tip of the probe was located at 3 mm ventral from the pial surface. After a slow and smooth descent, the probe was allowed to sit still at the target depth for at least 5 min before initiation of recording to allow the electrodes to settle.

Neural activity was filtered (high-pass at 300 Hz), amplified (200× gain), multiplexed and digitized (30 kHz) and recorded using the SpikeGLX 3.0 software (https://github.com/billkarsh/SpikeGLX). Recorded data were preprocessed using an open-source software KiloSort 2.0 (https://github.com/MouseLand/Kilosort) and manually curated using Phy (https://github.com/cortex-lab/phy) to identify putative single units in each of the primary motor cortex and dorsolateral striatum. A total of six experimental sessions (from four mice; Extended Data Fig. 5) with simultaneous motor cortical and striatal recordings were included in this work. The average number of motor cortical units included for each mouse was: mouse 38, 98 ± 4 (range, 95–102); mouse 39, 64; mouse 40, 75 ± 5 (range, 70–80); and mouse 44, 55. The average number of striatal units included for each mouse was: mouse 38, 100 ± 13 (range, 87–112); mouse 39, 108; mouse 40, 74 ± 5 (range, 69–79); and mouse 44, 110.

Data analysis

We used a similar approach for both monkey and mouse data. In all the analyses, we only considered the trials in which the animal successfully completed the task within the specified time and received a reward. We concatenated trials in time for subsequent analyses—that is, no trial-averages were taken. For the monkey centre-out reaching task and the mouse reaching and pulling task, an equal number of trials to each target was randomly selected (eight targets for the monkeys and two targets for mice, except in Extended Data Fig. 6a, for which four targets were considered). Trial order was randomized to eliminate the possible effect of the passage of time. Within each trial, we isolated a window of interest that captured most of the movement, starting 50 ms before movement onset and ending 400 ms after movement onset. To analyse covert behaviour in monkeys, we used a window that spanned the movement planning period, which started 400 ms before movement onset and ended 50 ms after movement onset. Importantly, all of our results held when changing the analysis windows within a reasonable range.

For the monkey random-walk task, each reach could start and end anywhere within the workspace. To define movements (conditions) that could be matched across animals, we first segmented the workspace into 12 circular subsections. Each subsection was then divided into six equal sectors and targets in the same angular sector were grouped together, creating 72 possible target conditions. We separated the sequences of four consecutive reaches and considered each reach as a separate movement. To assign each movement to a target condition, we first assigned each movement to one of the subsections on the basis of the starting position of the given movement, excluding movements that started more than 2 cm from the centre of the subsection. We then recentred the movements so that they started in the centre of each subsection and reached outwards towards their target. The movement was then assigned to a sector and target condition on the basis of the angle of target. To study the preservation of latent dynamics across monkeys performing similar behaviour, we needed to match movements (reach conditions) across sessions for different monkeys. To maximize the number of matched movements, we compared all sessions for Monkey M and Monkey T against a reference session for Monkey C_L that had the most successful trials. We matched movements in each pair of sessions by minimizing the mean squared error (MSE) between pairs of movements, excluding matches that had MSEs above the threshold of 2% of MSEs calculated for all possible pairs of movements. If the matched movements had different corresponding target conditions, we used the target condition label from the reference session. After this process was completed, we excluded target conditions with less than six matched movements, such that paired sessions had up to 29 shared target conditions. Because these movements were more ballistic than in the centre-out task, we examined a window starting at movement onset and ending 350 ms after movement onset.

All the analyses were implemented in Python using open-source packages such as numpy, matplotlib, sci-kit, scipy and pandas^{70,71,72,73,74} and custom code. As we were analysing existing datasets on an individual basis, no explicit planning of sample size, group randomization or blinding was performed.

Behavioural correlation

To assess the behavioural stereotypy of a given animal, we calculated hand trajectory correlations (Pearson’s correlation) of every pair of trials within a session (Extended Data Fig. 1b and Extended Data Fig. 5b). The distributions in Fig. 2k inset illustrate these correlations pooled across all the monkey centre-out and mouse reaching and pulling sessions included in this work. To determine the behavioural similarity across pairs of sessions from different monkeys or mice (Fig. 2k), we similarly calculated correlations to compare all pairs of trials from the two sessions.

Neural population latent dynamics

To estimate the latent dynamics associated with the recorded neural activity in each session for both mice and monkeys, we computed a smoothed firing rate as a function of time for each unit. We obtained these smoothed firing rates by applying a Gaussian kernel (σ = 50 ms) to the binned square-root transformed spike counts (bin size 30 ms) of each unit. We excluded units with a low mean firing rate (less than 1 Hz mean firing rate across all bins) but we did not perform any further exclusions, for example, based on lack of modulation or behavioural tuning. For each session, this produced a neural data matrix X of dimension n by T, where n is the number of recorded units and T the total number of time points from all concatenated trials on a given day; T is thus given by the number of targets per day × number of trials per target × number of time points per trial. We performed this concatenation as described above after randomly subselecting the same number of trials for all targets for each animal (15 trials for monkey centre-out, six for monkey random walk, 22 for mouse reaching and pulling). For each session, the activity of n recorded units was represented as a neural space—an n-dimensional sampling of the space defined by the activity of all neurons in that brain region. In this space, the joint recorded activity at each time bin is represented as a single point, the coordinates of which are determined by the firing rate of the corresponding units. Within this space, we estimated the low-dimensional latent dynamics by applying PCA to X. This yielded n PCs, each a linear combination of the smoothed firing rates of all n recorded units. These PCs are ranked on the basis of the amount of neural variance that they explain. We defined an m-dimensional, session-specific manifold by only keeping the leading m PCs, which we referred to as neural modes. We chose a manifold dimensionality m = 10, based on previous studies examining motor cortical recordings during upper limb tasks^5,26,46. Across all datasets, a ten-dimensional manifold explained about 60% of the neural variance for each of the monkey motor cortex (Extended Data Fig. 1c), mouse motor cortex and mouse striatum (Extended Data Fig. 5e). Note, however, that our results held within a reasonable range of dimensionalities, similar to refs. ^26,33,46 (Extended Data Figs. 2f and 4b). We computed the latent dynamics within the manifold by projecting the time-varying smoothed firing rates of the recorded neurons onto the m = 10 PCs that span the manifold. This produced a data matrix L of dimensions m by T.

Aligning latent dynamics through CCA

We addressed our hypothesis that different animals performing the same behaviour would share preserved latent dynamics by aligning the dynamics using CCA^26,75. CCA was applied to the latent dynamics of each pair of sessions afterconcatenating the same number ofrandomly ordered trials to each target(condition, in the case of the sequential reaching task). For details on using CCA to align latent dynamics, see ref. ²⁶.

We measured the similarity in latent dynamics across animals by computing the across-animal correlations as the canonical correlations (CCs) across all pairs of sessions from any two different monkeys or mice. To establish the strength of the across-animal correlations, we computed an upper bound defined by the within-animal correlations, which we calculated as the 99th percentile of the correlations between two randomly selected subsets of trials within any given session over 1,000 samples. The ‘control’ correlations represent a lower bound for the CCs. We computed these by shuffling the targets across the two sessions and using a randomly selected control window (more details in the ‘control analyses’ section below) in each trial, rather than the movement or preparatory epochs.

Note that to summarize each comparison to a single datapoint (for example, in Fig. 2k and Extended Data Figs. 2hand 6d), we computed the mean of the top four CCs of the latent dynamics²⁶. In Fig. 2k, we used this approach to establish a relationship between the strength of preservation of the latent dynamics and the consistency of behaviour, quantified as the mean trajectory correlation of all possible pairs of trials across two animals. Furthermore, when showing pairs of ‘aligned’ trajectories across animals, such as inFig. 2e and Extended Data Fig. 3, the CCA axes were made orthogonal using singular value decomposition for visualization purposes.

Finally, we showed that preserved latent dynamics could be uncovered across a broad range of manifold dimensionalities. In Extended Data Fig. 2f we repeated the alignment analysis for manifold dimensionalities m = 2–19.

Decoding analysis

To test whether the aligned latent dynamics maintain movement-related information, we built standard decoders to predict hand trajectory during overt behaviour. If the aligned latent dynamics across different animals were behaviourally relevant, they would allow predicting time-varying hand trajectories even if the methods used to identify them (PCA and CCA) are not supervised, that is, they do not attempt to optimize decoding performance. We compared the predictive accuracy of three different types of decoders: (1) a within-animal decoder trained and tested (using ten-fold cross-validation) on two non-overlapping subsets of trials from each session of each animal; (2) an across-animal ‘aligned’ decoder that was trained on the aligned dynamics from one animal and tested on another, a comparison we performed on each pair of sessions from two different animals; (3) an across-animal ‘unaligned’ decoder that was trained on the latent dynamics from one animal and tested on another without aligning the dynamics using CCA. We also performed a similar analysis to predict the upcoming target during covert movement preparation in monkeys (Fig. 4f).

Hand trajectory decoders were LSTM models with two LSTM layers, each with 300 hidden units, followed by a linear output layer. The models were implemented with Pytorch⁷⁶ and trained for ten epochs with the Adam optimizer, with a learning rate of 0.001. Upcoming target classifiers were Gaussian Naïve Bayes models¹² (the GaussianNB class in ref. ⁷²). We included three bins of recent latent dynamics history, for a total of 90 ms, in the input of both the decoders and the classifiers. These extra neural inputs incorporate information about intrinsic neural dynamics and account for transmission delays. The R² value, defined as the squared correlation coefficient between actual and predicted hand trajectories, was used to quantify decoder performance. Moreover, in Extended Data Fig. 4d we verified that our choice of across-animal decoder accuracy metric did not influence the observation that preserved latent dynamics are informative about behaviour by also computing a variance accounted for (VAF) metric, defined as:

where y_i represents the actual value of the predicted variable, ŷ_i its predicted value and \(\bar{y}\) its mean. For this analysis, we normalized hand trajectories by the length of the reaches (determined by the 99th percentile of their hand positions along each axis) because monkeys had workspaces of different sizes.

The hand trajectory was a two-dimensional signal in monkeys and a three-dimensional signal in mice. We built separate decoders to predict hand trajectories along the x, y (and z for mice) axes. We then reported the average performance across all axes. For target classification, we reported the mean accuracy of the classifier (the score() method).

To test how many dimensions of the aligned latent dynamics were needed for accurate across-animal decoding of behaviour, we repeated the decoding analysis in the monkey centre-out dataset for manifold dimensionalities m = 1, 2…,14 (Extended Data Fig. 4b).

Finally, we performed a control analysis to ensure our across-animal decoding results were not biased by sharing similar trials for both alignment and decoder training. We split the full dataset of one animal into three non-overlapping sets: one to align the latent dynamics, one to train the decoder and one to test the performance across animals. Extended Data Fig. 4c shows the result of this analysis for the monkey centre-out data. Despite having aligned the latent dynamics only using half of the data, the impact on decoding performance is negligible.

Control analyses

Alignment of latent dynamics with random behavioural windows

To establish a ‘behaviourally irrelevant’ window as control, we randomly selected windows of similar length to our behavioural windows (450 ms) along the entire duration of the intertrial and trial periods combined. This ensured we had samples of dynamics in the neural population with realistic statistics but that they were not directly coupled to shared behaviour across individuals. We used this window to provide a lower-bound control for the alignment of neural population latent dynamics (‘control’inFigs. 2f,g,j, 3d,e and 4b,e and Extended Data Figs. 2b–d,g, 3 and 8d).

Aligning latent dynamics through Procrustes analysis

We used CCA to align the latent dynamics in all the analyses. However, to ensure that our results hold regardless of the specific method used for alignment, we replicated the main result using Procrustes analysis⁷⁷. Procrustes finds the best transformation that minimizes the sum of the squares of the differences between the two input datasets. Following a procedure identical to the CCA analysis, we aligned the dynamics from two different datasets using Procrustes analysis (the scipy.spatial.procrustes class in ref. ⁷³) and then correlated the aligned dynamics to yield a metric comparable to that of the CCA (Extended Data Fig. 2g,h). Note that the degrees of preservation of latent dynamics obtained with CCA and Procrustes analysis are largely similar.

Neural variance explained by aligned latent dynamics

We measured the percentage of neural variance explained by the preserved latent dynamics using a method we devised in ref. ³³. Briefly, we ‘reconstructed’ the preserved neural activity by projecting the aligned latent dynamics along the CC axes back to the PCA space (the neural manifold)and then to the original neural state space. We then measured the difference between the total neural variance and the variance of these reconstructed signals using an approach similar to that in ref. ⁷⁸. By repeating this procedure iteratively for an increasing number of manifold dimensions m, we measured the neural variance explained by each dimension of the aligned latent dynamics. Using this approach, we found that preserved latent dynamics explain a significant fraction of the neural population variance (Extended Data Fig. 2e).

Surrogate datasets with TME

We established a lower-bound control by aligning the latent dynamics from randomly selected windows sampled across different task conditions and behavioural epochs (see above). In addition to this control, we also used TME to generate surrogate neural data as another lower-bound control²⁹. TME produces surrogate data that preserve the second-order statistics of the actual neural data (that is, covariance across time, across neurons or across experimental conditions) but are otherwise random (Extended Data Fig. 2a). Aligning these surrogate data through the same procedure as the original data shows significantly lower correlations for monkey centre-out task, random-walk task and mouse reaching and pulling task (Extended Data Fig. 2b–d).

Aligning topological structure in neural population activity

To test whether the topological structure in the produced movements is sufficient to produce preserved latent dynamics, we quantified the degree of similarity in latent dynamics across individuals that could be uncovered when aligning the static, topological features of the neural population activity, rather than the dynamics of the movements, using a technique developed in ref. ²⁶. To align the topological structure of neural population activity, we time-averaged the activity for each neuron during the execution epoch of each trial in the monkey centre-out reaching task. We then analysed the time-averaged data with the previous methodology by performing PCA to find a neural manifold and using CCA to align each pair of sessions (Extended Data Fig. 7a). This procedure led to well-aligned ‘topological representations’ (example in Extended Data Fig. 7b). To directly test whether aligning the topological structure of neural population activity is sufficient to uncover preserved latent dynamics, we projected the latent dynamics on the CC axes found through this (static) topological alignment and calculated the pairwise correlations of the resultant projected latent dynamics. These correlations were remarkably lower than those obtained through alignment of the time-varying latent dynamics (Extended Data Fig. 7c,d).

Control analyses on the numbers of conditions and neurons

To establish that the preservation of latent dynamics holds across different degrees of task complexity, we calculated the correlations for increasing numbers of subsampled target conditions for each pair of sessions in the monkey random target task (Fig. 3f and Extended Data Fig. 8b). We randomly subsampled different combinations of target conditions and calculated the degree of preservation of the latent dynamics for up to 10,000 combinations for each number of conditions.

To establish that preserved latent dynamics can be uncovered regardless of the specific measured neurons, we also calculated the correlations for varying numbers of neurons in the random target task (Fig. 3g and Extended Data Fig. 8c). For each pair of sessions, we eitherrandomlykept neurons (Fig. 3d) or randomlydropped neurons(Extended Data Fig. 8c) in increments of ten until we ran out of measured neurons for either session and repeated this process 50 times, calculating the degree of preservation at each step. For both analyses, we calculated the mean correlations for the top four CCs across all subsamples for each pair of sessions.

Comparison of different but related tasks

The central hypothesis of this study is that preserved latent dynamics are the basis for the generation of similar behaviour across individuals from the same species. Here, we sought to further support this hypothesis by showing that the latent dynamics produced by two individuals engaged in the same task are more similar than the latent dynamics produced by the same individual performing two different but related tasks. To this end, we compared our results to our previous study on the relationship of neural population activity underlying different but related wrist manipulation or reach-to-grasp tasks in monkeys³³(Extended Data Fig. 9).

Recurrent neural network models

Model architecture

To show that the preservation of latent dynamics across animals engaged in the same task is not a trivial consequence of similar behaviour, we trained RNNs to perform the same centre-out reaching task as the monkeys. These models were implemented using Pytorch⁷⁶. Similar to previous studies simulating motor cortical dynamics during reaching^27,79,80,81, we implemented the dynamical system \(\dot{{\bf{x}}}=F\left({\bf{x}},{\bf{s}}\right)\) to describe the RNN dynamics:

where x_i is the hidden state of the ith unit and r_i is the corresponding firing rate following tanh activation of x_i. All networks had N = 300 units and I = 3 inputs, a time constant τ = 0.05 s and an integration time step dt = 0.01 s. The noise η was randomly sampled from the Gaussian distribution \({\mathscr{N}}(\mathrm{0,0,2})\) for each time step. Each unit had an offset bias, b_i, which was initially set to zero. The initial states x_t=0 were sampled from the uniform random distribution \({\mathscr{U}}\left(-\mathrm{0.2,0.2}\right)\). All networks were fully recurrently connected, with the recurrent weights J initially sampled from the Gaussian distribution \({\mathscr{N}}\left(0,\frac{g}{\sqrt{N}}\right)\), where g = 1.2. The time-dependent inputs s fed into the network had input weights B initially sampled from the uniform distribution \({\mathscr{U}}\left(-\mathrm{0.1,0.1}\right)\). These inputs consisted of a one-dimensional fixation signal which started at 2 and went to 0 at the go cue and a target signal that remained at 0 until the visual cue was presented. The two-dimensional target signal (2 cos θ^target, 2 sin θ^target) specified the reaching direction θ^target of the target.

The networks were trained to produce two-dimensional outputs p corresponding to x and y positions of the experimentally recorded reach trajectories, which were read-out via the linear mapping:

where the output weights W were sampled from the uniform distribution \({\mathscr{U}}\left(-\mathrm{0.1,0.1}\right)\).

Model training

Networks were trained to generate positions of reach trajectories using the Adam optimizer⁸² with a learning rate l = 0.0005, first moment estimates decay rate β₁ = 0.9, second moment estimates decay rate β₂ = 0.999 and ϵ = 1 × 10^–8. The loss function L was defined as the MSE between the two-dimensional output and the target positions over each time step t, with the total number of time steps T = 400. The first 50 time steps were not included to allow network dynamics to relax:

To examine whether two networks could have different latent dynamics while producing the same motor output, we devised a network with more constraints to perform the behavioural task with distinct latent dynamics (Fig. 5a). We added a loss term that penalized the CC between the latent dynamics of the ‘constrained’ network being trained and those of another previously trained ‘standard’ network during movement execution:

where c_i is the ith CC. To examine different degrees of preserved latent dynamics, we trained the networks at varying values of α = 0, 0.05, 0.25 or 0.50.

Networks were trained until the average loss of the last ten training trials fell below a threshold of 0.2 cm², for at least 50 and up to 500 training trials, with a batch size B = 64. Each batch had equal numbers of trials for each reach direction. We clipped the gradient norm at 0.2 before the optimization step. Both standard and constrained training were performed on ten different networks initialized from different random seeds. The same set of random seeds was used for constrained networks at different values of α.

Connectivity analyses

By increasing the value of α, we were able to decrease the preservation of the latent dynamics while keeping behavioural performance constant. To examine how this changed the underlying connectivity, we calculated the variance and dimensionality of the weight changes in the recurrent weights J following training (Fig. 5f,g).

Statistics and reproducibility

We compared the performance of various within-animal and across-animal movement decoders and classifiers using two-sided Wilcoxon’s rank sum tests. We replicated the core findings across two species (mice and monkeys), four behaviours (a centre-out reaching task, a sequential reaching task and a reach, grasp and pull task, along with during covert movement planning) and two brain regions (motor cortex and dorsolateral striatum). Experiments on each species were performed independently in two different laboratories and by different scientists. The mice experiments were done in a single cohort, whereas the monkey data were collected in two sets of experiments (one for the centre-out task, another for the random reaching task), each spanning 2 years. Overall, our neural recordings and behavioural data are in good agreement with related published studies. All attempts at replication were successful.

Reporting summary

Further information on research design is available in theNature Portfolio Reporting Summary linked to this article.