Interplay between external inputs and recurrent dynamics during movement preparation and execution in a network model of motor cortex

The primary motor cortex has been shown to coordinate movement preparation and execution through computations in approximately orthogonal subspaces. The underlying network mechanisms, and the roles played by external and recurrent connectivity, are central open questions that need to be answered to understand the neural substrates of motor control. We develop a recurrent neural network model that recapitulates the temporal evolution of neuronal activity recorded from the primary motor cortex of a macaque monkey during an instructed delayed-reach task. In particular, it reproduces the observed dynamic patterns of covariation between neural activity and the direction of motion. We explore the hypothesis that the observed dynamics emerges from a synaptic connectivity structure that depends on the preferred directions of neurons in both preparatory and movement-related epochs, and we constrain the strength of both synaptic connectivity and external input parameters from data. While the model can reproduce neural activity for multiple combinations of the feedforward and recurrent connections, the solution that requires minimum external inputs is one where the observed patterns of covariance are shaped by external inputs during movement preparation, while they are dominated by strong direction-specific recurrent connectivity during movement execution. Our model also demonstrates that the way in which single-neuron tuning properties change over time can explain the level of orthogonality of preparatory and movement-related subspaces.


Introduction
The activity of the primary motor cortex (M1) during movement preparation and execution plays a key role in the control of voluntary limb movement [1,2,3,4,5,6,7]. Classic studies of motor preparation were performed in a delayed-reaching task setting, showing that firing rates correlate with task-relevant parameters during the delay period, despite no movement occurring [8,9,10,11,12,13,14,14,15,16,17]. More recent works have shown that preparatory activity is also displayed before non-delayed movements [18], that it is involved in reach correction [19], and that when multiple reaches are executed rapidly and continuously, each upcoming reach is prepared by the motor cortical activity while the current reach is in action [20]. Preparation and execution of different reaches are thought to be processed simultaneously without interference in the motor cortex through computation along orthogonal dimensions [20]. Indeed, the preparatory and movement-related subspaces identified by linear dimensionality reduction methods are almost orthogonal [21] so that simple linear readouts that transform motor cortical activity into movement commands will not produce premature movement during the planning stage [22]. However, response patterns in these two epochs are nevertheless linked, as demonstrated by the fact that a linear transformation can explain the flow of activity from the preparatory subspace to the movement subspace [21]. How this population-level strategy is implemented at the circuit level is still under investigation [23]. A related open question [24] is whether inputs from areas upstream to the primary motor cortex (such as from the thalamus and other cortical regions, here referred to as external inputs) that have been shown to be necessary to sustain movement generation [25] are specific to the type of movement being generated throughout the whole course of the motor action, or if they serve to set the initial conditions for the dynamics of the motor cortical network to evolve as shaped by recurrent connections [26,27,28,29,30,31].
In this work, we use a network modeling approach to explain the relationship between network connectivity, external inputs, and computations in orthogonal dimensions. Our analysis is based on electrophysiological recordings from M1 of a macaque monkey performing a delayed center-out reaching task. The dynamics of motor cortical neurons during reaching limb movements has been shown to be low-dimensional [32]. Here, we develop a low-dimensional description of the dynamics using order parameters that quantify the covariation between neural activity and the direction of motion (see [33,34,35,36] but also [37,38]). Recorded neurons are tuned to the direction of motion both during movement preparation and execution, but their preferred direction and amplitude of the tuning function change over time ( [39,40,41]). Interestingly, major changes happen when the activity flows from the preparatory to the movement-related subspaces. We describe neuronal selectivity during the task by four parameters: two angular variables corresponding to the preferred direction during movement preparation and execution, respectively; and two parameters that represent the strength of tuning in the two epochs. We characterized the empirical distribution of these parameters, and investigated potential network mechanisms that can generate the observed tuning properties by building a recurrent neural network model, whose synaptic weights depend on tuning properties of pre and post-synaptic neurons, and external inputs can contain information about movement direction. First, we analytically derived a low-dimensional description of the dynamics in terms of a few observables denoted as order parameters, which recapitulate the temporal evolution of the population-level patterns of tuning to the direction of motion. Then, we inferred the strength of recurrent connections and external inputs from data, by imposing that the model reproduce the observed dynamics of the order parameters. There are multiple combinations of feedforward and recurrent connections that allow the model to generate neural activity that strongly resembles the one from recordings both at the single-neuron and population level, and that can be transformed into realistic patterns of muscle activity by a linear readout. To break the model degeneracy, we imposed an extra cost associated with large external inputs -that likely require more metabolic energy consumption compared to local recurrent inputs. The resulting solution suggests that different network mechanisms operate during movement preparation and execution. During the delay period, the population activity is shaped by external inputs that are tuned to the preferred directions of the neurons. During movement execution, the localized pattern of activity is maintained via strong direction-specific recurrent connections. Finally, we show how the specific way in which neurons tuning properties rearrange over time produces the observed level of orthogonality between the preparatoryand movement-related subspaces.

Subjects and Task
We analyzed multi-electrode recordings from the primary motor cortex (M1) of two macaque monkeys performing a previously reported instructed-delay, center-out reaching task [42]. The monkey's arm was on a two-link exoskeletal robotic arm, so that the position of the Figure 1: (a) Schematic of the center-out delayed reaching task and definition of the preparatory (blue) and of the movement-related (red) epochs. Black circles represent the time of: target onset; go cue; start of movement; end of movement, averaged across trials. (b) Example of the trial averaged firing rate of two neurons during movement preparation as a function of the target location. Solid lines represent the corresponding cosine fit. For one example neuron, we show its preferred direction (θ A ), defined as the location of the peak of the cosine function; and its degree of participation (η A ), which is proportional to the amplitude of the cosine function. (c) Same as in (a), but for the activity of the same neurons during movement execution. The preferred direction is denoted as θ B and the degree of participation is η B . (d) Scatter plot of the preferred direction during execution (θ B ) vs preparation (θ A ) for all neurons; preferred directions are defined as the location of the peak of the cosine tuning functions (circular correlation coefficient r = 0.4). The heat map represents the joint distribution over (θ A , θ B ) inferred from the data (Eq. 4). (e) Scatter plot of the normalized amplitude of the cosine tuning curve during execution (η B ) vs preparation (η A ) for all neurons (correlation coefficient r = 0.3). Blue line: empirical distribution of η A from kernel density estimation; red line: empirical distribution of η B . (f-g) Trial averaged firing rate for all neurons during movement preparation as a function of θ A (f) and θ B (g), for the condition with target location = 5π/4, indicated by the dashed vertical line. The activity is normalized across condition for each neuron. The activity of one neuron chosen at random is indicated by the orange dots in panels f and g. monkey's hand controlled the location of a cursor projected onto a horizontal screen. The task consisted of three periods (Fig. 1a): a hold period, during which the monkey was trained to hold the cursor on a center target and wait 500 ms for the instruction cue; an instruction period, during which the monkey was presented with one of eight evenly spaced peripheral targets and continued to hold at the center for an additional 1,000-1,500 ms; a movement period, signalled by a go cue, when the monkey initiated the reach to the peripheral target. Successful trials for which the monkeys reached the target were rewarded with a juice or water reward. The peripheral target was present on the screen throughout the whole instruction and movement periods. In line with previous studies (e.g., [21]), the preparatory and movement-related epochs were defined as two 300ms time intervals beginning, respectively, 100ms after target onset and 50ms before the start of the movement.

Patterns of correlation between neural activity and movement direction
Studies of motor cortex have shown that tuning to movement direction is not a timeinvariant property of motor cortical neurons, but rather varies in time throughout the course of the motor action; single-neuron encoding of entire movement trajectories has been also reported (e.g., [39,40]). We measured temporal variations in neurons preferred direction by binning time into 160ms time bins and fitting the binned trial-averaged spike counts as a function of movement direction with a cosine function. As the variability in preferred direction during the delay period alone and during movement execution alone was significantly smaller than the variability during the entire duration of the task ( Fig. 1-figure supplement 1), we characterized neurons tuning properties only in terms of these two epochs. This simplifying assumption is discussed further in Methods and Discussion. We will show later that such a simplification is enough for the model to recapitulate neuronal activity both at the level of single-units and at the population level.
Figs. 1b-c show two examples of tuning curves, where the trial-averaged and timeaveraged activity during movement preparation and execution is plotted as a function of the location of the target on the screen, together with a cosine fitting curve. In the two epochs, neurons change their preferred direction -denoted, respectively, as θ A and θ B -and their degree of participation, which is proportional to the amplitude of the cosine tuning function -denoted as η A and η B . A scatter plot of the preferred directions in preparatory (θ A ) and execution (θ B ) epochs for all neurons is shown in Fig. 1d, while a similar scatter plot for the degrees of tuning in preparatory (η A ) and execution (η B ) epochs is shown in Fig. 1e. Neurons preferred direction in the two epochs are moderately correlated (circular correlation coefficient r = 0.4, p = 3 10 −5 ), and so are their degree of participation, even if to a lesser degree (correlation coefficient r = 0.3, p = 10 −4 ). Neuronal tuning to movement direction is reflected in a population activity profile that is spatially localized, as shown in Figs.1.f-g, were we plot the normalized firing rate for all neurons, as a function of their preferred direction, separately for the two epochs. An example of the activity of a single neuron in the two epochs is highlighted in red. It illustrates how the activity of single neurons is drastically rearranged across time, while the population activity remains localized around the same angular location, which corresponds to the location of the target on the screen.

Recurrent neural network model
The experimental observations described above motivated us to build a network model in which neuronal selectivity properties match the ones observed in the data. In particular, we built a network model in which neurons are characterized by the four same parameters we use to fit the preferred directions θ A , θ B and degrees of tuning, η A , η B , leading to a 4dimensional selectivity space. The subspace defined by the coordinates (θ A , η A ) is denoted as map A, while the subspace defined by (θ B , η B ) is denoted as map B. We will denote the Figure 2: (a) Schematic depiction of the network model. Each shaded disk corresponds to one unit in the network. The direction of the blue and red arrows within each disk represents the neurons preferred direction during preparatory (θ A ) and execution (θ B ) epochs, respectively, while the thickness of the arrows represents the degree of participation (η A , η B ) in these epochs. The orange arrows show the strength of synaptic connectivity (defined by equation 5) between units. It contains both a symmetric component, that depends on the distance between preferred directions of pre and post-synaptic neurons separately for the two epochs (see strong connections between the two top left neurons), and an asymmetric component that depends on the distance between the preferred preparatory direction of the pre-synaptic neuron and the preferred execution direction of the post-synaptic one (see connections between bottom right neurons). (b) Spatially localized activity of the network at a given time t, shown in three different 2-dimensional subspaces of the 4-dimensional selectivity space. Left: activity plotted on map A (θ A , η A ) at fixed θ B , η B . Center: activity plotted on map B (θ B , η B ) at fixed θ A , η A . Right: activity plotted as a function of θ A , θ B , at fixed η A , η B . (c) Network activity at a given time t, plotted as a function of θ A , at fixed θ B , η B and for different values of η A . (d) Network activity at a given time t, plotted as a function of θ B , at fixed θ A , η A and for different values of η B . coordinate vector by and refer to the units in the network as neurons, even though each unit in the model could instead represents a group of M1 neurons with similar functional properties, and likewise the connection between two units in our model could represent the effective connection between the two functionally similar groups of neurons in M1.
In this model, neurons with coordinates (selectivity parameters) x are described by their firing rate r(x, t) whose temporal evolution is given by: where τ is the time constant of firing rate dynamics. [ ] + is the threshold-linear (a.k.a. relu) transfer function that converts synaptic inputs in firing rates, and I tot (x; t) is the total synaptic input to neurons with coordinates x. We set the time constant to τ = 25ms, which is of the same order of magnitude of the membrane time constant, and we checked that for values of τ in the range 10ms − 100ms our results did not quantitatively change. The total input to a neuron is The first term in the r.h.s. of (2) is the recurrent input, which depends on the firing rates of presynaptic neurons r(x , t), and on J(x, x ), the strength of recurrent connections from neurons with coordinates x to neurons with coordinates x. I ext (x; t) is the external input, and ρ(x) is the probability density of x. The recurrent term in I tot is the sum of single neuron contributions: here, we assumed that the cortical network is sufficiently large that instead of summing over the contributions of single neurons, each one at a given coordinate x = {θ A , θ B , η A , η B }, we can integrate over a continuous distribution of x. The probability density of the coordinates ρ(x) is set to match the empirical distribution of the preferred direction and degree of participation. The preferred directions are not significantly correlated with the degrees of tuning ( Fig. 1-figure  supplement 2), and we therefore took them to be independent: The distribution of preferred directions was well fitted by: as shown in Fig. 1-figure supplement 2, while for the sake of simplicity we assumed ρ p (η A , η B ) = ρ pA (η A ) ρ pB (η B ) and estimated ρ pA , ρ pB non-parametrically using kernel density estimation. The last two ingredients that we need to specify to define the network dynamics are the recurrent couplings and external inputs. The strength of synaptic connections from a pre-synaptic neuron with preferred directions (θ A , θ B ) and participation strengths (η A , η B ) to a post-synaptic neuron with preferred directions (θ A , θ B ) and participation strengths (η A , η B ) is given by where j 0 represents a uniform inhibitory term; j A s and j B s measure the amplitude of the symmetric direction-specific connections in map A and map B, respectively; and j a measures the amplitude of asymmetric connections from map A to map B. In the Methods, we explain how the dynamics of a network with such recurrent connectivity can reproduce the spatiallylocalized population activity that we observed in the data, and we provide an intuition for the role of the different coupling parameters. A schematic depiction of the network is shown in Fig. 2a. We parameterized the external input in analogy with the recurrent input: where C 0 represents an untuned homogeneous external input, C A and C B represent untuned map-specific inputs to maps A and B, respectively, A and B represent directionnally tuned inputs to maps A and B, respectively, and Φ ext is the direction encoded by external inputs. An example of the spatially-localized population activity at a given time t visualized in different 2-D subspaces of the 4-D space with coordinates x = {θ A , θ B , η A , η B } is shown in Fig. 2.b, blue and an example of tuning curves from the network's activity is shown in Fig. 2.c.

Order parameters
The simplicity of the model described by Eqs. (1-6) allowed us to derive a low dimensional description of the model dynamics in terms of a few population-level signals denoted as order parameters. The order parameters quantify the average population activity r 0 (t), and the degree to which the population activity is localized in map A (r A (t)) and B (r B (t)). These parameters are defined by the following equations: Note that r A and r B are the first Fourier coefficient of the population rate over the domain θ A and θ B , respectively.
These order parameters can be both computed in the model and from data. Therefore, they can provide us with a simple tool to compare model and data, and to infer model network parameters from data. Fig. 3a shows the dynamics of the order parameters computed from the data during the delayed reaching task. As expected, we see that during the preparatory period the activity activity is localized in map A (r A (t) > 0) but not in map B (r B (t) ∼ 0). Around the time of the go cue, network activity dynamically reorganizes, such as the network becomes now strongly localized in map B (r B (t) > 0), while the degree of modulation in map A slowly decreases towards zero. We can think of map A and map B as two different coordinate systems that the network can use to produce distinct patterns of population activity, that are both localized around the location of the target on the screen, but in different ways. A spatially localized activity profile associated with either r A (t) > 0 or r B (t) > 0 is denoted as a bump. A definition of this term provided in the Methods. At different times, the bump of activity has different shapes -either localised in map A, or in map B, or in both maps -but its location remains constant. Next, we studied possible mechanisms underlying these observations.

Tuned states in the model: External inputs vs recurrent connectivity
We next investigated the network model to understand the mechanisms underlying the observed dynamics. We first derived equations governing the temporal evolution of the order parameters (see Methods). We then analyzed the stationary solutions of the equations in the absence of tuned inputs (C A = C B = A = B = 0), in the space of the three parameters defining the couplings strength: j 0 , j A s , j B s . We found that, similarly to the ring model [43,44], there exists two qualitatively distinct regions in the space of parameters Phase diagram of the model, shown as a function of the parameters j A s , j B s , that describe how strongly maps A and B are embedded in recurrent synaptic connectivity. Different curves correspond to bifurcation lines for different values of the parameter j a , modulating the asymmetric term in synaptic connectivity. The area underneath each line (gray) corresponds to the homogeneous phase; the area beyond each line (white) corresponds to the phase where the network exhibit a localized activity pattern even in absence of tuned external inputs. (c) Weak coupling scenario: Here, the observed dynamics of the order parameters results from external inputs that are tuned to θ A during movement preparation and to θ B during movement execution. (d) Strong coupling scenario: Here, strong recurrent connections sustain the localized pattern of activity; a change in untuned external inputs makes the activity to be localized along map A during preparation, and along map B during preparation. characterizing recurrent connectivity (Fig. 3b). For weak recurrent connectivity (gray area in Fig. 3b), the activity in the network is uniform in the absence of tuned inputs. Thus, in this region, tuned network activity must rely on external inputs. For strong recurrent connectivity (white area in Fig. 3b), network activity is tuned, even in the absence of tuned inputs. In this case, the location of the bump in network activity is determined by initial conditions. These two regimes are separated by a bifurcation line, where the uniform solution becomes unstable due to a Turing instability [43,44].
This analysis shows that a network activity profile that is localised in a given map, say map A (r A (t) > 0), can be sustained either thanks to the strong recurrent connectivity (i.e. a strong enough parameter j A s in Eq. 5), or thanks to the external input term proportional to A (t) (Eq. 6). Hence, our model can potentially generate the observed dynamics of the order parameters (Fig. 3a) for different choices of recurrent connections and external inputs parameters, ranging in between the following two opposite scenarios: • Recurrent connections are absent: j A s = j B s = j a = 0. In this case, the system is simply driven by feedforward inputs. The localized activity is the result of external inputs that selectively excites or inhibits specific neurons during motor preparation and execution, thanks to sufficiently large values of the parameter A (t) > 0 during preparation and B (t) > 0 during execution. This scenario is depicted schematically in Fig. 3c.
• The network is strongly recurrent, with strongly direction-specific connections, and external inputs are untuned ( A = 0, B = 0). Analytical study of the dynamics (Methods) shows that, in order for such a system to exhibit tuning, the strength of the couplings has to exceed a critical value shown in Fig. 3b. Moreover, when the synaptic strength exceeds this value, the activity can be localized in map A during preparation and in map B during execution simply as the result of homogeneous external inputs changing their strength, without the need for tuned external inputs to selectively excite/inhibit specific neurons. This scenario is depicted in Fig. 3d.
We next turn to the question of which of these two scenarios best describes the data.
Fitting the model to the data Our next goal is to infer the strength of recurrent connections (j 0 , j A s , j B s ) and feedforward inputs (C 0 (t), C A (t), C B (t), A (t), B (t)) that best describes the data. To do so, we imposed that our network reproduce the dynamics of the population-level order parameters (r 0 , r A , r B ). Specifically, for a given set of recurrent connections and feedforward inputs, we can compute analytically the dynamics of the order parameters (Methods). Based on these calculations, we build an iterative procedure that minimizes the reconstruction error E rec , i.e. the squared difference between the predicted and observed dynamics of the order parameters. The results show that there is a large region of the (j 0 , j A s , j B s )-space where the reconstruction error has essentially the same value ( Fig. 4-figure supplement 2). Solutions range from a a network with zero recurrent connections (i.e. a purely feedforward scenario, shown in Fig. 4.a-c) to a network with strong direction-specific connections.
To break the degeneracy of solutions, we added an energetic cost to the reconstruction error. This reflects the idea that a biological network where the computation is only driven by long-range connections from upstream areas is likely to consume more metabolic energy than a network where the computation is processed through shorter-range recurrent connections. Our inference procedure minimizes the cost function: where E ext is the total magnitude of external inputs (see Methods, eq. 31). The hyperparameter α describes the relative strength of these two terms. For very small values of α, the algorithm mainly minimizes external inputs, and the resulting reconstruction of the dynamics is poor. For very large values of α, multiple solutions are found, that yield similar dynamics (see Fig. 4c and Fig. 4i). Fig. 4-figure supplement 1 shows the results we obtained by using intermediate values of α, that are small enough to yield a unique solution, but large enough to guarantee a good reconstruction of the dynamics. Interestingly, all solutions cluster in a small region of the parameter space, that is close to the bifurcation surface in such a way that the parameter j B s is close to (either larger or smaller) its critical value; j A s is much smaller than its critical value, and j a is zero. We show the results obtained with two distinct values of α in Fig. 4d and Fig. 4g; the corresponding dynamics of the inferred external inputs in Fig. 4e and Fig. 4f; and the resulting dynamics of the order parameters from mean field equations in Fig. 4f and Fig. 4i. In particular, Fig. 4d corresponds to a smaller value of α, i.e. we are penalizing more strongly for large external inputs: the resulting coupling parameters are stronger than their critical value. Fig. 4g corresponds to a larger value of α: the strength of the couplings is below the critical line, yielding larger external inputs (Fig. 4h) and a smaller reconstruction error ( Fig. 4.i). The purely feedforward case is added in Fig. 4.a-c for comparison. While a purely feedforward network requires a strong input tuned to map B right before the movement onset ( Fig. 4b), in our solution such input is either absent (Fig. 4e) or much weaker than the corresponding untuned input (Fig. 4h). Our analysis therefore suggests that recurrent connectivity plays a major role in maintaining the degree of spatial modulation observed in the data.
Importantly, the same analysis applied to a second dataset recorded from a different macaque monkey performing the same task yielded qualitatively similar results (see

The model generates realistic neural and muscle activity
We have shown that the model is able to reproduce the dynamics of the order parameters computed from data. Here, we ask how the activity of single neurons in our model compares to data, and whether a readout of neural activity can generate realistic patterns of muscle activity. We simulated the dynamics of a network of 16000 neurons. To each neuron i, we assigned the coordinates θ B so as to match the empirical distribution of coordinates (eq 3 and Fig. 5-figure supplement 3.a). We added a noise term to the total current that each neuron is subject to, modelled as an Ornstein-Uhlenbeck process (see Methods, eq. 33).
First, we checked that the dynamics of the order parameters -previously computed by numerically integrating the analytical equations -is also correctly reconstructed by simulations ( Then, we computed tuning curves during the preparatory and execution epochs from simulations, and estimated the values θ Next, we showed that the network activity closely resembles the one from recordings. At the level of the population, we used canonical correlation analysis to identify common  patterns of activity between simulations and recordings, and to show that they strongly correlate (see Fig. 5.a and Methods). At the level of single units, for each neuron in the data we selected a neuron in the model network with the closest value of θ A , θ B , η A , η B variables. A side-by-side comparison of the time course of responses shows a good qualitative agreement (Fig. 5.c). We also noted that both the activity from recordings and from simulations present sequential activity and rotational dynamics (Fig. 5-figure supplement  4). Finally, a linear readout of the activity from simulations can generate realistic patterns of muscle activity, which closely match electromyographic (EMG) signals recorded from Monkey Bx during a center-out reaching movement ( Fig. 5.b).
So far in this section, we discussed results from a network with strong couplings; Fig. 5-figure supplement 5 shows that the dynamics of a purely feedforward network is almost indistinguishable from from the strongly recurrent case -although canonical correlations between the activity of a purely feedforward network and the data are slightly lower than the ones for networks with strong recurrent couplings, during movement execution (average canonical correlation of 0.74 for the purely feedforward network, and 0.82 for the strongly recurrent one). We also noticed that networks with coupling parameters below the bifurcation surface are robust to noise in the θ variables. Conversely, the solutions above the bifurcation surface that do not require tuned external inputs are very sensitive to such noise, and the location of the bump of the activity drifts towards a few discrete attractors. When the noise is weak, the drift happens on a time scale that is much larger that the time of movement execution; the larger the level of the noise, the faster the drift.

PCA subspaces dedicated to movement preparation and execution
In the previous sections, we considered the low-dimensional description of the dynamics given by the order parameters. Here, we use the commonly used principal component analysis (PCA) to compare the dimensionality of data and network model. In particular, we ask whether our formalism can explain the orthogonality of the preparatory and movementrelated linear subspaces. As previously reported in the literature [22,21], the top panel in Fig. 6.a shows that the top principal components of the preparatory activity explain most of the variance of the preparatory activity, but very little variance of the movement-related activity; vice versa, the top movement-related principal components explain very little variance of the preparatory activity ( Fig. 6.a, bottom panel). The activity of our network model also displays this property ( Fig. 6.b). We also quantified the level of orthogonality of the two subspaces using the alignment index as defined in [21], that measures the percentage of variance of each epoch's activity explained by both sets of PCs. Fig. 6.c shows that the alignment index is much smaller than the one computed between two subspaces drawn at random (random alignment index, explained in the Methods), in both data and network simulations. Note that model with uniform degrees of participation (η A = η B = 1) displays a higher overlap between the preparatory and movement-related subspaces, as shown in Fig. 6-figure supplement 2.c. In the Methods section, we explain how the dimensionality of the activity depends on the connectivity of the network.

Discussion
Orthogonal spaces dedicated to movement preparation and execution Studies on the dynamics of motor cortical activity during delayed reaching tasks have shown that the primary motor cortex employs an 'orthogonal but linked strategy' [22, 21] to coordinate planning and execution of movements.
In this work, we explored the hypothesis that this strategy emerges as result of a specific recurrent functional architecture, in which synaptic connections store information about two distinct patterns of activity that underlie movement preparation and movement execution.
We built a network model in which neurons are characterized by their selectivity properties in both preparatory and movement execution epochs (preferred direction and degree of participation), and in which synaptic connectivity is shaped by these selectivity properties in a Hebbian-like fashion.
A strong correlation between the selectivity properties of the preparatory and movementrelated epochs will produce strongly correlated patterns of activity in these two intervals and a strong overlap between the respective PCA subspaces. We inferred the distribution of these tuning features from data and showed that the correlation between the preparatory and movement-related patterns of activity is small enough to allow for almost orthogonal subspaces, which is thought to be important for the preparatory activity not to cause premature movement. At the same time, the correlation is non-zero, and that allows the activity to flow from the preparatory to the movement-related subspaces with minimal external inputs, as summarized in the next section.

Interplay between external and recurrent currents
We analytically described the temporal evolution of the population activity in the lowdimensional space defined by maps A and B in terms of a few order parameters, which can be easily computed from data. Different combinations of the strength of direction-specific recurrent connections and of tuned external inputs allow the model to accurately reproduce the dynamics of the order parameters from data. We argue that solutions that require less inputs from areas upstream of the motor cortex are favorable in terms of metabolic energy consumption. With the addition of a cost proportional to the magnitude of external inputs, we find solutions where recurrent connections are strong and direction specific. In the resulting scenario, during movement preparation, an external input tuned to map A sustains a population-level activity localized in map A, and pins the location of the peak of activity. During movement execution, the localized activity is sustained mostly by recurrent connectivity; the correlation between preferred directions θ A and θ B allows the activity in map B during movement execution to be localized around the same location as it was in map A during movement preparation. The inferred strength of recurrent connectivity is close to the critical value above which the recurrent network can generate localized patterns of activity in the absence of tuned inputs. Solutions well beyond the bifurcation line require an implausible fine tuning of the recurrent connections, as heterogeneity in the connectivity causes a systematic drift of the encoded direction of motion on a typical time scales of seconds -the larger the structural noise in the couplings, the faster the drift, as has been extensively studied in continuous attractor models [45,46,47,48,49,50]. It has been shown that homeostatic mechanisms could compensate for the heterogeneity in cellular excitability and synaptic inputs to reduce systematic drifts of the activity [47] and that short-term synaptic facilitation in recurrent connections could also significantly improve the robustness of the model [48] -even when combined with short-term depression [49]. While a full characterization of our model in the presence of structural heterogeneity is beyond the scope of this work, we note that solutions close but below the bifurcation line are stable with respect to perturbations in the couplings. In this case, tuned inputs are present also during movement execution, but their magnitude is much weaker than the untuned ones: direction-specific couplings are strong and amplify the weak external inputs tuned to map B, therefore playing a major role into shaping the observed dynamics during movement execution.
Our prediction that external inputs are direction-specific during movement preparation but mostly non-specific during movement execution needs to be tested experimentally. Interestingly, it agrees with several previous studies on the activity of the primary motor cortex during limb movement. In particular, [30] showed that changes in neural activity that characterize the transition from movement preparation to execution reflect when movement is made but are invariant to movement direction and type. [51] measured the correlations in firing between motor thalamic (m-Th) and motor cortical cells, in monkeys performing the same task considered in the present work. In order to investigate if the observed cofiring patterns resulted from neurons tuning properties, the authors looked at the effect of the difference in preferred directions on cofiring patterns between pairs of M1-mTh cells. They reported significant correlation between cofiring patterns and difference in preferred direction solely before movement, but not during movement execution (see Appendix, Fig. S8 of ref. [51]). While this result is in line with our prediction on the role of recurrent vs feedforward inputs, further simultaneous recordings of M1 and upstream regions, as well as measures of synaptic strength between motor cortical neurons, will be necessary to rule out alternative explanations.

Cosine tuning
While cosine tuning functions have been extensively used to describe the firing properties of motor cortical neurons (e.g., [52]), and they were also hypothesized to be the optimal tuning profile to minimize the expected errors in force production [53], more recent work has emphasized that tuning functions in the motor cortex present heterogeneous shapes. Specifically, the authors of [54] showed that tuning functions are well fitted by a sum of a cosine-modulated component, and an unstructured component, often including terms with higher spatial frequency. They also argued that the unstructured component is key for a readout of motor cortical activity to reproduce EMG activity. In this work, we showed that a noisy recurrent network model that is based on cosine tuning reproduces well the R 2 coefficient of the cosine fit of tuning curves from data (see a comparison between Fig.  1-figure supplement 2.e-f and Fig. 5-figure supplement 2.c). Moreover, tuning curves from simulations present heterogeneous shapes (Fig. 5-figure supplement 3), including bimodal profiles, especially for neurons with low values of η A , η B . Indeed, we showed that the R 2 of the cosine fit strongly correlates with the variables η, both in the data and in the model. Neurons with a low degree of participation have low R 2 coefficient of the cosine fit, and present heterogeneous tuning curves. We also showed that a linear readout of the activity of our model can very accurately reconstruct EMG activity. Finally, our model could easily be extended to incorporate other tuning functions. While we leave the analysis of a model with an arbitrary shape of the tuning function to future work, we don't expect our results to strongly depend on the specific shape of the tuning function.

Comparison with the ring model
The idea that the tuning properties of motor cortical neurons could emerge from directionspecific synaptic connections goes back to the work of Lukashin and Georgopulos [55]. However, it was with the theoretical analysis of the so called ring model [56,43,57,44] that localized patterns of activity were formalized as attractor states of the dynamics in networks with strongly specific recurrent connections. Related models were later used to describe maintenance of internal representations of continuous variables in various brain regions [58,59,60,61,62,63,50,64,65,66] and were extended to allow for storage of multiple continuous manifolds [67,68,69] to model the firing patterns of place cells the hippocampus of rodents exploring multiple environments. While our formalism builds on the same theoretical framework of these previous works, we would like to stress two main differences between our model and the ones previously considered in the literature. First, we studied the dynamic interplay between fluctuating external inputs and recurrent currents, that causes the activity to flow from the preparatory map to the movement-related one and, consequently, neurons tuning curves and order parameters to change over time, while at the population level the pattern of activity remains localized around the same location. Then, we introduced an extra dimension representing the degree of participation of single neurons to the population pattern of activity; this is an effective way to introduce neuronto-neuron variability in the responses, which decreases the level of orthogonality between the preparatory and movement-related subspaces, and yields tuning curves whose shape resembles the one computed from data -in contrast with the classic ring model, where all tuning curves have the same shape.

Representational vs dynamical system approaches
There has been a debate whether neuronal activity in motor cortex is better described by so-called 'representational models' or dynamical systems models [40,70]. Representational models [70,71] are models in which neuronal firing rates are related to movement parameters [1,52,72,73,74,75,76,77,78,79,80]. Dynamical systems models [29] are instead recurrent neural networks whose synaptic connectivity is trained in order to produce a given pattern of muscle activity. Such models have been argued to reproduce better the dynamical patterns of population activity in motor cortex [70].
In our model, firing rates are described by a system of coupled ODEs, and the synaptic connectivity is built from the neuronal tuning properties, in the spirit of the classical ring model discussed above, and of the model introduced by Lukashin and Georgopoulos [55]. Our model is thus a dynamical system where kinematic parameters can be decoded from the population activity. Importantly, our model is constrained to reproduce the dynamics of a few order parameters that are a low-dimensional representation of the activity of recorded neurons. In contrast to kinematic-encoding models, our model can recapitulate the heterogeneity of single-unit responses. Moreover, as in trained recurrent network models, a linear readout of the network activity can reproduce realistic muscle signals.
The advantage of our model with respect to a trained RNN is that it yields a lowrank connectivity matrix which is simple enough to allow for analytical tractability of the dynamics. The model can be used to test specific hypotheses on the relationship between network connectivity, external inputs and neural dynamics, and on the learning mechanisms that may lead to the emergence of a given connectivity structure. The model is also helpful to illustrate the problem of degeneracy of network models. An interesting future direction would be to compare the connectivity matrices of trained RNNs and of our model.

Extension to modeling activity underlying more complex tasks
Neurons directional tuning properties have been shown to be influenced by many contextual factors that we neglected in our analysis [81,82], to depend on the acquisition of new motor skills [83,84,85] and other features of movement such as the shoulder abduction/adduction angle even for similar hand kinematic profiles [37,38], or the speed of movement for similar trajectories [40]. A large body of work [38,86,87,81,88,89,90] has shown that the activity in the primary motor cortex covaries with many parameters of movement other than the hand kinematics -for a review, see [91]. More recent studies [70,92,90,93] have also suggested that the largest signals in motor cortex may not correlate with task-relevant variables at all.
Our model can be extended to more realistic scenarios in several ways. A simplifying assumption we made is that the task can be clearly separated into a preparatory phase and one movement-related phase. A possible extension is one where the motor action is composed of a sequence of epochs, corresponding to a sequence of maps in our model. It will be interesting to study the role of asymmetric connections for storing a sequence of maps. Such a network model could be used to study the storage of motor motifs in the motor cortex [94]; external inputs could then combine these building blocks to compose complex actions. Moreover, incorporating variability in the degree of symmetry of the connections could allow to model features that are not currently included, such as the speed of movement.
In summary, we proposed a simple model that can explain recordings during a reaching task. It provides a scaffold upon which more sophisticated models could be built, to explain neural activity underlying more complex tasks.
[7] Andrew R Brown and G Campbell Teskey. Motor cortex is functionally organized as a set of spatially distinct representations for complex movements.

Circular statistics
To define the statistical measures of angular variables [95, 96], let us consider a set of angles {θ i } i , and visualize them as vectors on the plane: The mean vector is then defined as from which we can get the mean angular direction θ using the four quadrant inverse tangent function. The length of the mean resultant vector, is a measure of how concentrated the data sample is around the mean direction: the closer R is to 1, the more concentrated the data is. We quantify the spread in a data set through the circular variance given by which ranges from 0 to 1. To compute the correlation between two sets of angular variables, {θ i } i and {φ} i , we used the correlation coefficient [97] defined as: where θ and φ are the mean angular directions of the two sets.

Preparatory and movement-related epochs
In our analysis, we characterized neurons' tuning properties during movement preparation and execution. The choice of considering only two epochs is an approximation based on the observation that neurons preferred directions are more conserved within the delay period alone and within movement execution alone than across periods ( Fig. 1-figure supplement  1). In Fig. 1-figure supplement 1, we considered three temporal epochs: the delay period (blue lines), movement execution (red) and the whole duration of the task (grey). We binned each epoch into N W time bins of length 160 ms: the preparatory and execution epochs consist of N W = 3 bins each, while the whole task consists of N W = 7 bins. For each neuron, we first computed its preferred direction in each time bin, by fitting the trialaveraged and time-averaged firing rate with a cosine function. Then, for each epoch, and for each neuron, we computed the circular variance of preferred directions across the N W bins. The cumulative distribution of circular variances is shown in Fig. 1-figure supplement 1.a.: the variability of preferred direction within the delay epochs and within movement execution is non-zero, although their distribution is skewed towards zero. For each epoch, we computed the median variability in preferred direction and used a bootstrap analysis to check its statistical significance. Let us consider the matrix of single trial activity of a given neuron, for a specific condition, and during a given epoch; the size of the matrix is N W × n, where n is the number of trials per condition (in the data, n ranges between 41 and 47). We repeat the same procedure for all of the 8 conditions, and z-scored the 8 activity matrices across conditions (this is done because the average activity changes across time bins, and in the following we will shuffle time bins). To check if tuning curves are exactly conserved across time-bins, we generated B = 1000 bootstrap samples of the activity matrix for a given condition, where each entry is chosen at random (with repetitions) between the possible N W × n entries of the original matrix. For each bootstrap sample matrix, we computed the variance of preferred direction as we did for the original matrix. This is repeated for all neurons, all epochs, all conditions. The cumulative distribution of all variances of preferred direction is shown in Fig. 1-figure supplement  1.b. Then, for each bootstrap sample, we computed the corresponding median variance of preferred direction; the histogram is shown in Fig. 1-figure supplement 1.c. From  Fig. 1-figure supplement 1.c, it is obvious that the median variance of preferred directions from the data is significantly larger than the one from the bootstrap samples. This suggests that neurons do change their preferred direction within epochs, even though major changes happen between epochs.

Analysis of the model
We studied the rate model with couplings defined by (5) with two complementary approaches. First, an analytic approximation based on mean-field arguments and valid in the limit of large network size allowed us to derive a low-dimensional description of the network dynamics in terms of a few latent variables, and to fit the model parameters to the data. Next, we checked that simulations of the dynamics of a network of N = 10 4 neurons reproduce the results that we derived in the limit N → ∞.
The mean-field equations are derived following the methods introduced in [43,68]. In the limit where the number of neurons is large, the average activity of a neuron with coordinates θ A , θ B , η A , η B is described by equations (1) and (2), where we have defined In order to derive a lower dimensional description of the dynamics, we rewrite (1) in terms of the average activity rate (14), and of the second Fourier components of the activity rate modulated by η A/B : The phase ψ A(/B) is defined so that the parameter r A(/B) is a real nonnegative number, yielding eq 14 for r A(/B) , and the following equation for the phase ψ A(/B) , representing the position of the peak of the activity profile in map A(/B): From (1), we see that the order parameters evolve in time according to the following set of equations, where the total input (2) is rewritten as: .

Stationary states for homogeneous external inputs
We first characterize the model by studying the fixed-points of the network dynamics when subject to a constant and homogeneous external input, and focusing on the scenario where the joint distribution (3) of the θ A , θ B is of the form which fits the empirical distribution of the data well for x = 2/3 ( Fig. 1-figure supplement  2). For now, we leave the distributions ρ pA (η A ) and ρ pB (η B ) unspecified. We first consider the case where the external input to the network is a constant that is independent of θ A , θ B : The stationary solutions of (1,16) are of the form: where we have defined Here, {r 0 , r A , r B } are solutions of the system (15) with the left hand side set to zero. The second term on the r.h.s in the last equation of (20) is obtained from (16) by observing that in the stationary state ψ A = ψ B if either j a > 0 or if θ A and θ B are correlated (as we are assuming here). As in [43,44], we can distinguish broad from narrow activity profiles. The term broad activity profile refers to the scenario where the activity of all of the neurons is above threshold, the dynamics is linear and the stationary state reduces to: By inserting the above equation in (14), we find that the only solution is homogeneous over the maps θ A and θ B : where the notation . represents an average over the measure dµ, i.e.
First, we notice that a nonzero homogeneous state is present if Then, the stability of this state with respect to a small perturbation {δr 0 , δr A , δr B } can be studied by linearizing (15) around the stationary solution, at fixed ψ A = ψ B = 0. The resulting Jacobian matrix is The homogeneous solution (22) is stable if j 0 < 1 and if the couplings parameters {j A s , j B s , j a } satisfy the following system of inequalities: If j 0 ≥ 1, the system undergoes an amplitude instability. For values of {j A s , j B s , j a } that exceed the threshold implicitly defined by (23), there is a region of the four-dimensional space (θ A , θ B , η A , η B ) where the total input I tot (19) is negative. The dynamics is no longer linear and the activity profile at the fixed point is localized around a particular direction, and characterized by positive stationary values of the order parameters r A , r B . This 'bump' of activity can be localized either more strongly in map A (r A > r B > 0), in map B (r B > r A > 0) or at the same level in both maps (r A = r B > 0). Since maps A and B are correlated according to 3, the location of the stationary bump of activity is constrained to be the same in the two maps: ψ A = ψ B ≡ ψ, with arbitrary ψ. The asymmetric term in the connectivity, modulated by the parameter j a , further contributes to aligning the location of the bump in map B to the one in map A. The system can relax to a continuous manifold of fixed points parameterized by ψ that are marginally stable. In this continuous attractor regime, the system can store in memory any direction of motion ψ, as the activity is localized in absence of tuned inputs.
Examples of two-dimensional cross-sections of the four-dimensional activity profile in the marginal phase are shown in Fig. 2.b. Fig. 2.c shows one-dimensional cross-sections of the activity as a function of θ A or θ B , for different values of η A and η B . These correspond to the tuning curves of the model for stationary external inputs. The tuning profile can have either a full-cosine modulation or a rectified cosine modulation, depending on the values of η A , η B . The plots of Fig. 1.f-g correspond to the activity as a function of θ A , θ B in the case of a discrete number of neurons -each neuron with a different value of η A , η B .

Time dependent external input
The activity of motor cortical neurons shows no signature of being in a stationary state but instead displays complex transients. To model the data, we assumed that the neurons are responding both to recurrent inputs and to fluctuating external inputs that can be either homogeneous or tuned to θ A , θ B , with peak at constant location Φ A = Φ B ≡ Φ (eq. 6). The total input 2 that neurons are subject to at a given time t is: where In order to draw a correspondence between the model and the data, we note that the activity of each neuron in the data corresponds to the activity rate at a specific coordinate θ A , θ B , η A , η B in the model. From equations (1) and (24) we see that, if we assume that changes in the external inputs happen on a time scale larger than τ , the modulation of the activity profile in map A at each time t is cos(θ A − Φ) and has amplitude proportional to η A , while the modulation of the activity profile in map B is cos(θ B − Φ) and has amplitude proportional to η B . Accordingly, for each recorded neuron θ A and θ B represent the location of the peak of the neuron's tuning curve computed during the preparatory and movementrelated epoch, respectively, while η A and η B are proportional to the amplitude of the tuning curves.
Fitting the model to the data Neurons activity rates were computed by smoothing the spike trains with a Gaussian kernel with s.d. of 25ms and averaging them across all trials with the same condition; r (i) k (t) denotes the rate of neuron i at time t for condition k, each condition corresponding to one of the 8 angular locations of the target on the screen. Since trials had highly variable length, we normalized the responses along the temporal dimension before averaging them over trials, as follows. We divided the activity into three temporal intervals: from the target onset to the go cue; from the go cue to the start of the movement; from the start of the movement to the end of the movement. For each interval, we normalized the response times to the average length of the interval across trials. We then aggregated the three intervals together. We defined the preparatory and execution epochs -denoted by A and B -as two 300ms time intervals beginning, respectively, 100ms after target onset and 50ms before the start of the movement, in line with [21]. Fig. 3-figure supplement 1 shows that our results do not change qualitatively when the lengths of the preparatory and execution intervals are increased. For each neuron i, we fitted the activity rate averaged across time within each epoch as a function of the angular position Φ of the target with a cosine function: where the parameters θ (i) B represent the neuron's preferred direction during the preparatory and execution epochs. In our the model, the direction modulation of the rates, see (25), is proportional to η A\B , that measures how strongly the neuron participates in the two epochs of movement; hence, we defined η (i) A\B to be proportional to the amplitude of the tuning curve: The scatter plot of η A , η B (Fig. 1.e) shows an outlier with η B ∼ 1. We checked that our results hold true if we discard that point. Fig. 1-figure supplement 2 shows that both η A and η B strongly correlate with the R 2 -coefficient of the cosine fit. In other words, neurons with a higher value of η A , η B are the ones whose tuning curves more strongly resemble a cosine (this holds true in our simulations too, see Fig. 5-figure supplement 2.c, and Discussion). The order parameters (14) are computed at each time t by approximating the integrals with the sums: where N is the number of neurons and n c = 8 is the number of conditions. The angular location of the localized activity ψ (k) A/B (t) can be computed from (14) as However, this estimate is strongly affected by the heterogeneity in the distribution of θ A , θ B -deviating from the rotational symmetry of the model. That is why in Fig. 2 we approximated ψ (k) A/B (t) by the k − th angular location of the target on the screen. We checked that computing ψ (k) A/B (t) by using either method does not affect the dynamics of the order parametersr data A/B (t). We assumed that the observed dynamics of the order parameters r 0 (r), r A (t), r B (t) obeys equations (15,16) with time-dependent external inputs of the form (6). We inferred the value of the parameters {C 0 (t), C A (t), C B (t), A (t), B (t)} t of the external currents and j 0 , j A s , j B s , j a of the coupling matrix that allow us to reconstruct the dynamics of the order parameters r 0 , r A , r B computed from data; since this inference problem is undetermined, we required as further constraint that the model reconstruct the dynamics of the following two additional order parameters: In this way, for given coupling parameters j 0 , j A s , j B s , j a , we can uniquely identify the external currents parameters that produced the observed dynamics. Still, an equally good reconstruction of r 0 , r A , r B , r 0A , r 0B can be obtained for different choices of coupling parameters (2). Hence, we inferred the model parameters by minimizing a cost function composed of two terms: one that is proportional to the reconstruction error of the temporal evolution of the order parameters and the other that represents an energetic cost penalizing large external inputs.
Fitting the model to the data: details.
The fitting procedure was divided in the following steps: 1. The time interval T going from the target onset till the end of the movement was binned into ∆t = 5ms time bins: T = {∆t 1 , ∆t 2 , . . . ∆t T }.
3. For each time step ∆t i , i = 2, . . . , T : • We started from the reconstructed order parameters at the previous time step: and we let the dynamical system (15, 29) with external currents parameters C 0 , C A , C B , A , B evolve for ∆t = 5ms to estimate the order parameters at the current time step: • We inferred the value of the external currents parameters by minimizing the reconstruction error: that quantifies the difference between the order parameters estimated from the data and the reconstructed ones; note that the dependence of the cost function E on C 0 , C A , C B , A , B is implicitly contained in the reconstructed order parameters. We minimized the cost function (30) by using an interior point method algorithm [Byrd et al, 1999] starting from the initial condition we imposed that A > 0, B > 0 and added a L1 regularization term to the external inputs, not to infer pathologically large positive and negative inputs that balance each other. This favors solutions with smaller external inputs and non-zero reconstruction error with respect to ones with larger inputs and almostzero reconstruction error.
The external currents inferred with step 3 depend on our initial choice of of the couplings parameters j 0 , j A s , j B s , j a .
4. Using step 3, the value of the couplings parameters j 0 , j A s , j B s , j a is inferred by minimizing the cost function E tot = αE rec + E ext composed of two terms: the reconstruction error and a term that favors small external currents:

Simulations of the model
We simulated the dynamics of a finite network of N neurons. To each neuron i, we assigned the variables θ B as follows.
• In the θ A\B -space, we sampled N θ points {θ in such a way that their joint distribution matches the distribution of the data (4), as shown in Fig. 5-figure supplement 2.a.
• In the η A\B -space, we drew N η points {η at random from the empirical distribution ρ pA (η A )ρ pB (η B ).
• For i = 1, 2, . . . N θ , we assigned to a block of N η neurons the same coordinates θ The network dynamics we simulated is defined by the following stochastic differential equation: where . (34) W in (33) is a Wiener process, and the parameters γ = 75/s, σ n = 0.35Hz/ √ s set the magnitude of the noise fluctuations. The level of noise is chosen so that the results of the PCA analysis (see next section) match the data. Note that by setting the noise to zero and taking the limit N → ∞ we recover the mean-field equations (1). The results of Fig.  3 are obtained from a network of N = 16000 neurons; we simulate the network dynamics for 8 location of the external input Φ and 10 trials for each of the 8 conditions, i.e. 10 instances of the noisy dynamics. The order parameters shown in Fig. 3 are computed from single trial activity, and then averaged over trials. The correlation-based analysis, instead, is obtained from trial-averaged activity.

Correlation-based analysis
The correlation-based analysis explained in this section was performed both on the smoothed and trial-averaged spike trains from recordings, and on the trial-averaged activity rates from simulations. Out of the 16000 units in the simulated network model, we considered a subset of the same size N = 141 as the recorded neurons. In particular, for each neuron in the data, we chose the corresponding one from simulations with the closest value of parameters θ A , θ B , η A , η B -the one with smallest euclidean distance in the space defined by (θ A , θ B , η A , η B ). To compute signal correlations, we preprocessed the data as follows, both for the recordings and for the simulations. For each neuron, we normalized the activity by its standard deviation (computed across all times and all conditions); then, we meancentered the activity across conditions. The T = 300ms long preparatory activity for all C = 8 conditions was concatenated into a N ×T C matrix denoted by P , and the movementrelated activity was grouped into an analogous matrix M .
CCA: We used a canonical correlation analysis (CCA) to compare the population activity from data and simulations. First, we projected the activity onto the K PCA dimensions that captured 90% of the activity variance across conditions; for the preparatory epoch, K = 20 for both the data and the simulations; for the movement-related epoch, K = 14 for the data, and K = 10 for the simulations. We then applied CCA to look for common patterns in the activity matrices from data and simulations. In brief, CCA finds a set of linear transformations of the activity matrices, in such a way that the transformed variables (called canonical variables) are maximally correlated. Fig. 5.a shows a high correlation between the sets of canonical variables from data and simulations. We repeated the procedure for the simulations of the purely feedforward network. In this case, the movement-related activity is slightly higher dimensional than for the recurrent network (K = 12), and the average canonical correlation is a bit lower. Note that both activities are trial-averaged, and the same realization of the noise was used in both sets of simulations (feedforward and recurrent).
PCA: We obtained correlation matrices relative to preparatory and movement related activity by computing the correlations between the rows of the respective matrices. We then identified the prep-PCs and move-PCs by performing PCA separately on the matrices P and M . The degree of orthogonality between the prep-and move-subspaces was quantified by the Alignment Index A [21], measuring the amount of variance of the preparatory activity explained by the first K move-PCs: where E mov is the matrix defined by the top K move-PCs, C prep is the covariance matrix of the preparatory activity and σ prep (i) is the i-th eigenvalue of C prep . K was set to the number of principal components needed to explain 90% of the execution activity variance. Hence, the Alignment Index ranges from 0 (orthogonal subspaces) to 1 (aligned subspaces). As random test, we computed the Random Alignment Index between two sets of K dimensions drawn at random within the space occupied by neural activity, using the Monte Carlo procedure described in [21]. We performed the same analysis on both the data (K = 14) and the model (K = 10) trial averaged activity. The dimensionality of the preparatory and movement-related subspaces can be understood as follows. The activity of the ring model encoding the value of a singular angular variable is two-dimensional in the Euclidean space. Similarly, the activity of the double-ring model encoding two distinct circular maps is four-dimensional. Our model is an extension of the double-ring model, where both the connectivity matrix and the external fields (eq. 5 and 6) are a sum of several terms, each one composed of an η-dependent term multiplying a θ-dependent term. The connectivity matrix is still rank-four, but its eigenvectors are modulated by the η variables. Although the dynamics is four-dimensional, we have shown that during movement preparation the activity is localized only in map A, while during movement execution it is predominantly localized in map B: in either epoch, we expect only two eigenvalues to explain most of the activity variance, as we indeed see from simulations in absence of additive noise (Fig. 6-figure supplement 1). The noise term introduces extra random dimensions; as the dynamics gets higher dimensional, both the alignment index and the random alignment index get smaller.
jPCA: We quantified the rotational structure in the data by applying the jPCA [26] dimensionality reduction technique to both simulated activity and recordings during movement execution. jPCA is a technique to identify the dimensions that capture rotational dynamics in the data. Given the matrix of activity M defined above, jPCA finds the best fit for the real skewed-symmetric matrix R that transforms the activity M into its derivativė M :Ṁ = RM.
The plane capturing the strongest rotations is spanned by the eigenvectors of the matrix R associated with the two largest eigenvalues. Fig. 5-figure supplement 4 shows that the trajectories from simulations qualitatively resembled the ones from data when projected onto the jPCA subspace that capture most of the variance.

EMG signals
The EMG signals were recorded from macaque monkey Bx during the execution of a centerout reaching movement, as described in  Fig. 5.b were rectified and smoothed; signals from different trials were rescaled to match the average length of movement execution, and then averaged over trials.
The linear readout of the activity from simulations was defined as where r j is the trial-averaged activity from simulations, with the C = 8 conditions concatenated into a matrix of dimensions N × T C, where T the average duration of movement execution; z is a matrix with dimension m × T C, m = 13 being the number of recorded muscles. For each muscle k, the EMG signal z k (t) was regressed onto the activity rate {r j (t)} j to find the readout weight vector {W k,j } N j=0 . In Fig. 5.b, the input firing rates are the activity of the N = 141 units in the network model that are closest to the recorded neurons in the (θ A , θ B , η A , η B ) space.    s , j B s inferred from data, through minimization of a cost function composed of two terms: E tot = αE rec + E ext . E rec represents the reconstruction error of the order parameters, while E ext is the magnitude of the external inputs required to sustain the activity. α is an hyperparameter of the fitting algorithm. Different colors correspond to the result of using a different value of the hyperparameter α in the cost function. The black line represents the bifurcation surface in the space j A s , j B s , at j a = 0.      Trial-averaged activity rates during movement execution (from 50ms before the start of the movement, to 250ms after the start of the movement) projected onto the first two jPCA dimensions, defined as in [26]. In order to compare the data to the theory, we averaged the activity of all trials corresponding to the same condition -defining 8 groups (each denoted in a different color) corresponding to the 8 different conditions. However, hand kinematics corresponding to the activity within the same group can differ quite substantially, even if they end up at the same target location. This is a major difference with respect to the analysis of [26], where many more groups were defined, each group containing the activity corresponding to very similar hand trajectories. Right. Same as in (a), but for activity from simulations. (b) Activity rates averaged across trials for condition 1, and z-scored across time, from data (left) and simulations (right). Neurons are sorted according to the time of the peak of activity.   Fig.7.g (strong couplings) when no noise is added to the dynamics. (b) Alignment index quantifying the degree of orthogonality between the preparatory and movement-related subspaces from simulations, for the same solution as in (a): strong couplings, no noise added to the dynamics. In absence of noise, the dynamics is confined onto a four-dimensional space; the random alignment index quantifies the alignment between two two-dimensional subspaces drawn at random within the four-dimensional space occupied by neural activity. Note that the alignment index is significantly smaller than the random test.  . Simulations of a model without η A , η B generate dynamics that is identical for all neurons with the same preferred direction -a part from the noise terms; this also yields tuning curves that have the same shape for all neurons. (b) Alignment index. In a model with no η A , η B , the preparatory and movement-related subspace overlap more than in the data.