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

  1. Ludovica Bachschmid-Romano  Is a corresponding author
  2. Nicholas G Hatsopoulos
  3. Nicolas Brunel  Is a corresponding author
  1. Department of Neurobiology, Duke University, United States
  2. Department of Organismal Biology and Anatomy, University of Chicago, United States
  3. Committee on Computational Neuroscience, University of Chicago, United States
  4. Department of Physics, Duke University, United States
  5. Duke Institute for Brain Sciences, Duke University, United States
  6. Center for Cognitive Neuroscience, Duke University, United States
6 figures and 1 additional file

Figures

Figure 1 with 2 supplements
Tuning to the direction of motion during movement preparation and execution in a delayed reaching task.

(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 (Equation 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.

Figure 1—figure supplement 1
Variance of preferred directions across time.

(a) Cumulative density function of the circular variance of preferred directions. We binned time into time bins of length 160 ms, and computed neurons’ preferred directions in each time bin. We then considered the delay period only (3 time bins) and computed the circular variance of preferred directions across time bins, for each neuron. Its cumulative distribution is plot in blue. We repeated the same procedure for the movement-related epoch (3 time bins, red) and for the whole duration of the experiment (9 time bins, gray). (b) Bootstrap analysis (see Methods for details). Cumulative density function of the circular variance of preferred directions as in (a), but considering trials from the bootstrap distribution. (c) Thick line: median (across neurons) circular variance of preferred directions from data, for the delay period. Histogram: distribution of median circular variances of preferred directions from bootstrap samples, for the delay period. The dotted line represents the 95% quantile of the distribution.

Figure 1—figure supplement 2
Distribution of tuning parameters from data.

(a) Distribution of θA-θB computed from recorded data (orange) and data generated from the distribution: 1+2cos(θA-θB)/3 (gray). (b) Variables ηA,ηB plotted as a function of θA-θB. The circular correlation is r=0.13 (P=0.3), for both ηA and ηB. (c) Variables ηA,ηB plotted as a function of θA. The circular correlation is r=0.08 (P=0.6) for ηA and r=0.06 (P=0.7) for ηB. (d) Variables ηA,ηB plotted as a function of θB. The circular correlation is r=0.01 (P=0.9) for ηA and r=0.1 (P=0.2) for ηB. (e) Scatter plot of the amplitude of a cosine function fitted to the tuning curve of the preparatory activity vs the R2 coefficient of the fit. The neurons with larger tuning amplitude are the ones with higher R2. (f) Same as in (c) but for movement-related activity.

Spatially localized activity of the network model.

(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.

Figure 3 with 1 supplement
Dynamics of the order parameters and phase diagram of the model.

(a) Dynamics of the order parameters rA (degree of spatial modulation of the activity in map A) and rB (degree of spatial modulation in map B) computed from data (solid line; shaded area: ± SEM across trials). Black dots on the x-axis represent the trial-averaged time of: target onset, go cue, start of movement and end of movement. (b) Phase diagram of the model, shown as a function of the parameters jsA,jsB, 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 ja, 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.

Figure 3—figure supplement 1
Results are independent on chosen duration of preparatory and movement-related epochs.

For two different definitions of the preparatory (blue segment) and movement-related (red segment) epochs: (a) vs (b), we show the corresponding distributions of θA,θB: (c) vs (e); the distributions of ηA,ηB: (d) vs (f); the dynamics of the order parameters: (g) vs (i); and the results of the PCA analysis: (h) vs (l).

Figure 4 with 3 supplements
Inferred dynamics of the external currents required to sustain the observed dynamics of the order parameters for different sets of couplings parameters jsA,jsB,ja.

(a–c) Scenario where neurons are connected by uniform inhibitory connections, in the absence of direction-specific couplings. (a) Couplings parameters jsA,jsB,ja are set to zero (orange dot on the phase diagram). The black line represents the bifurcation surface in the space jsA,jsB, at ja=0. (b) Dynamics of the external inputs inferred from the data. Gray line: homogeneous (untuned) input; light blue/red line: input that is specific to map A/B, but untuned to direction; input that is directionally tuned, and specific to map A/B. Black dots on the x-axes represent the trial-averaged time of: target onset, go cue, start of movement and end of movement. (c) Dynamics of the order parameters r0,rA,rB computed from data (thin line; shaded area: ± SEM across trials) and model (thick line). (d-i) Both couplings parameters and external inputs are inferred from data. (d) Solution where the coupling parameters are slightly above the bifurcation line (α=0.17). (d) Dynamics of the external inputs inferred from the data. Note that tuned inputs are drastically lower than in b, especially during movement execution. (f) Dynamics of the order parameters r0,rA,rB computed from data (thin line) and model (thick line). (g) Solution where the coupling parameters are below but close to the bifurcation line (see Figure 4—figure supplement 1, α=0.5). (h) Dynamics of the external inputs inferred from the data. (i) Dynamics of the order parameters r0,rA,rB computed from data (thin line) and model (thick line).

Figure 4—figure supplement 1
Value of the couplings parameters jsA,jsB inferred from data, through minimization of a cost function composed of two terms: Etot=αErec+Eext.

Erec represents the reconstruction error of the order parameters, while Eext 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 jsA,jsB, at ja=0.

Figure 4—figure supplement 2
Reconstruction error vs cost function.

Left: Reconstruction error Erec computed at different values of jsA and jsB, for ja=0. The reconstruction error quantifies the squared difference between the dynamics of the order parameters computed from data and the one predicted by the model for a given value of coupling parameters jsA and jsB. In a large region of the space, the reconstruction error has the same value. Right: Example of the cost function Etot=αErec+Eext, for a given value of the hyperparameter α, yielding one of the solutions of Figure 4—figure supplement 1.

Figure 4—figure supplement 3
Analysis based on data recorded from monkey Rj.

(a) Couplings parameters inferred from data; different points correspond to solutions for different values of the hyperparameter α of the fitting procedure. The black line represents the bifurcation surface in the space jsA,jsB, at ja=0. (b) Dynamics of the external inputs inferred from data. (c) Dynamics of the order parameters (solid line: data, shaded area: ± standard error across the population; dotted line: analytical prediction).

Figure 5 with 5 supplements
Comparisons between model and data, for parameters of Figure 4d.

(a) CCA analysis to compare the population activity from simulations and from recordings. First, we projected the activity onto the PCA dimensions that captured 90% of the activity variance, for both the data and the simulations; then, we applied CCA to look for common patterns in the activity matrices from data and simulations. Top: canonical correlations, activity during preparation. Bottom: canonical correlations, activity during execution. (b) Blue lines: electromyographic (EMG) signals from 13 muscles, recorded during a center-out reaching movement. Each panel corresponds to a different condition (location of the target on the screen). Black lines: patterns of muscle activity predicted by a linear readout of the activity of 1000 units drawn at random from the 16000 units in the network model (cross- validated NMSE = 0.0066). (c) For each neuron in the data, we chose the corresponding one from simulations with the closest value of parameters θA,θB,ηA,ηB. Top row: examples of trial averaged activity from data; bottom row: corresponding neurons from simulations. Different shades of blue correspond to the 8 different conditions.

Figure 5—figure supplement 1
Dynamics of the order parameters from simulations, for the coupling parameters indicated in red in panel (a), above but close to the bifurcation line.

(b) Dynamics of the order parameters r0,rA,rB computed from data (thin line; shaded area: ± SEM across trials) and numerical simulations (thick line) of the dynamics of a finite-size network model with additive noise. (c) Location of the bump of activity in map A from trial-averaged activity from simulations, for 8 conditions corresponding to different locations of the target on the screen. (d) Location of the bump of activity in map B.

Figure 5—figure supplement 2
Network architecture for simulations.

(a) To build the network architecture, we assigned to 16,000 neurons coordinates ηA,ηB (top) and θA,θB (bottom) drawn from their respective empirical distribution. (b) Scatter plot of the coordinates initially assigned to the neurons (x-axes) vs the corresponding variables that we computed from simulations of the network activity during the preparatory and movement-related epochs (y-axes). Points clustered on the diagonal show self-consistency of the model. (c) R2 coefficient of the cosine fit for tuning curves computed from simulation during preparation (top) and execution (bottom) as a function of ηA and ηB, respectively.

Figure 5—figure supplement 3
Examples of tuning curves during the preparatory (blue) and movement-related (red) epochs computed from simulations (a) and from data (b).
Figure 5—figure supplement 4
Rotational dynamics and sequential activity.

(a) Left. 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 Churchland et al., 2012. 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 Churchland et al., 2012, 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.

Figure 5—figure supplement 5
Simulations of a network with no recurrent connections.

(a) r0,rA,rB computed from data (thin line; shaded area: ± SEM of the population) and numerical simulations (thick line). (b) CCA analysis to compare the population activity from simulations and from recordings. Top: canonical correlations, activity during preparation. Bottom: canonical correlations, activity during execution. (c) For each neuron in the data, we chose the corresponding one from simulations with the closest value of parameters θA,θB,ηA,ηB. Top row: examples of trial averaged activity from data; bottom row: their corresponding neurons from simulations. Different colors correspond to the 8 different conditions.

Figure 6 with 2 supplements
Orthogonality of the preparatory and movement-related subspaces.

(a) Percentage of variance of the preparatory (blue) and movement-related (red) activity from data explained by the first 11 principal components calculated from preparatory (top) and movement-related (bottom) trial-averaged activity. (b) As in a, activity from simulations. (c) The alignment index quantifies the degree of orthogonality between two subspaces. Top: alignment index between the preparatory and movement-related activities computed from data, compared to the randomized test (random alignment index, distribution in light gray and average in dark grey). Bottom: alignment index computed from simulations.

Figure 6—figure supplement 1
PCA analysis of network dynamics in various conditions.

(a) PCA analysis on the results of simulations for the solution of Figure 4d. (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. (c) PCA analysis on the results of simulations for the solution of Figure 4d. (strong couplings) and (d) of Figure 4a (zero couplings).

Figure 6—figure supplement 2
Network with uniform degrees of tuning, ηA=ηB=1 for all neurons.

(a) Examples of trial averaged activity from data (top panel) and simulations (bottom panel). 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) PCA analysis and (c) alignment index computed from the results of simulations with ηA,ηB. The preparatory and movement-related subspace overlap more than in the data.

Additional files

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Ludovica Bachschmid-Romano
  2. Nicholas G Hatsopoulos
  3. Nicolas Brunel
(2023)
Interplay between external inputs and recurrent dynamics during movement preparation and execution in a network model of motor cortex
eLife 12:e77690.
https://doi.org/10.7554/eLife.77690