Optimal control of transient dynamics in balanced networks supports generation of complex movements

Abstract

Populations of neurons in motor cortex engage in complex transient dynamics of large amplitude during the execution of limb movements. Traditional network models with stochastically assigned synapses cannot reproduce this behavior. Here we introduce a class of cortical architectures with strong and random excitatory recurrence that is stabilized by intricate, fine-tuned inhibition, optimized from a control theory perspective. Such networks transiently amplify specific activity states and can be used to reliably execute multidimensional movement patterns. Similar to the experimental observations, these transients must be preceded by a steady-state initialization phase from which the network relaxes back into the background state by way of complex internal dynamics. In our networks, excitation and inhibition are as tightly balanced as recently reported in experiments across several brain areas, suggesting inhibitory control of complex excitatory recurrence as a generic organizational principle in cortex.

Figure 1. Dynamical systems view of movement planning and execution.

(a) A typical delayed movement generation task starts with the instruction of what movement must be prepared. The arm must then be held still until the go cue is given, upon which the movement is performed. (b) During the preparatory period, model neurons receive a ramp input (green). Following the go-cue, that input is withdrawn, leaving the network activity free to evolve from the initial condition set up during the preparatory period. Model neurons then exhibit transient oscillations (black) which drive muscle activity (red). (c) Black-box view on movement generation. Muscles (red, right) are thought to be activated by a population of motor cortical neurons (“neural dynamical system”, middle). To prepare the movement, this network is initialized in a desired state by the slow activation of a movement-specific pool of neurons (green, left).

Figure 2. Random inhibition-stabilized networks (SOCs).

(a) Schematic of a SOC. A population of rate units is recurrently connected, with strong and intricate excitatory pathways (red) that would normally produce unstable, chaotic activity. Stabilization is achieved through fine-tuned inhibitory feedback (blue). (b) Eigenvalue spectrum of the connectivity matrix of a SOC (black), and that of the chaotic random network from which it is derived (gray). Stability requires all eigenvalues to lie to the left of the dashed vertical line. Note the large negative real eigenvalue, which corresponds to the spatially uniform activity pattern. (c) Matrices of synaptic connectivity before (unstable) and after (SOC) stability optimization through inhibitory tuning. By design, the excitatory weights are the same in both matrices. Matrices were thinned out to 40×40 for visualization purposes. The bottom row shows the strengths of all the inhibitory input synapses to a single sample neuron, in the unstable network (gray) and in the corresponding SOC (black). (d) Distribution of inhibitory weights in the unstable network (10% connection density, gray peak at W_ij ≈ -3.18) and in the stabilized version (40% connection density, black). The mean inhibitory weight of all possible synapses is the same before and after optimization (≈ -0.318, gray and black arrowheads). (e,f) Spontaneous activity in the unstable network (e) and in the SOC (f), for four example units. Note the difference in firing rate scales.

Figure 3. Transient amplification in SOCs.

(a) The energy E evoked by N=200 orthogonal initial conditions (a₁, …, a_N) as the network evolves linearly (Δr_i=x_i) with no further input according to Equation 1. The energy (Equation 4) is normalized such that it equals 1 for an unconnected network (W=0) irrespective of the initial condition (dashed horizontal line). Each successive initial condition a_i is defined as the one that evokes maximum energy, within the subspace orthogonal to all previous input patterns a_j<i (Experimental Procedures). The black arrowhead indicates the mean, or the expected evoked energy E₀ when the neurons are initialized in a random activity state. (b) Dynamics of the SOC in the linear regime. Top: time-evolution of ‖Δr‖/√N - which measures the momentary spread of firing rates in the network above or below baseline - as the dynamics unfold from any of the 10 best or 10 worst initial states (same color code as in (a)). Initial states have a standard deviation of σ=1.5 Hz across the population. The dashed gray line shows σ × exp(-t/τ), i.e. the behavior of an unconnected pool of neurons. Bottom: sample firing rate responses of 10 randomly chosen neurons following initialization in state a₁ or a₁₉₉. The red line indicates the momentary population-averaged firing rate. (c,d) Same as in (b), now with the nonlinear gain function shown in (e). Unlike in the linear case, the dynamics now depend on the spread σ of the initial firing rates across the network (1.5 Hz in (c) as in (b), 2 Hz in (d)). The larger this spread, the longer the duration of the population transient. When σ>3 Hz, the network initiates self-sustained chaotic activity (not shown). (e) Single-unit input-output nonlinearity (solid line, Δr_i=g(x_i) given by Equation 2) and its linearization (dashed line, Δr_i=x_i).

Figure 4. SOCs agree with experimental data.

(a) Experimental data, adapted with permission from Churchland et al., 2012. Each trace denotes the trial-averaged firing rate of a single cell (two sample cells are shown here) during a delayed reaching task. Each trace corresponds to one of 27 different movements. Vertical scale bars denote 20 spikes/sec. The go cue is not explicitly marked here, but occurs about 200 ms before movement onset. (b) Time-varying firing rates of two neurons in the SOC, for 27 “conditions”, each characterized by a different collective steady state of preparatory activity (see text). (c) Experimental data adapted from Churchland et al., 2012, showing the first 200 ms of movement-related population activity projected onto the top jPC plane. Each trajectory corresponds to one of the 27 conditions mentioned in (a). (d) Same analysis as in (c), for the SOC.

Figure 5. Generation of complex movements through SOC dynamics.

(a) Firing rates versus time for 10 neurons of the SOC, as the system prepares and executes either of the two target movements (snake, left or butterfly, right). Five test trials are shown for each neuron. The corresponding muscle trajectories following the go cue are shown for the same five test trials (thin traces) and compared to the target movement (black trace and dots). (b) Same as in (a), for a weakly connected (untuned) random balanced network (Experimental Procedures).

Figure 6. Precise balance of excitation and inhibition in SOCs.

The network is initialized in state a₁ (left), a₁₀ (middle) or a₁₀₀ (right), and runs freely thereafter. The amplitude of the initial condition is chosen weak enough for the dynamics of amplification to remain linear (cf. Figure 3). (a) Temporal evolution of the Pearson correlation coefficient r_EI between the momentary excitatory and inhibitory recurrent inputs across the population. (b) Corresponding time course of the correlation coefficients between the network activity and the initial state, calculated from the activity of the entire population (black), of the excitatory subpopulation (red), and of the inhibitory one (blue). (c) Temporal evolution of the correlation coefficient between the network activity when initialized in state a_i (i=1 (left), 10 (middle) or 100 (right)) and when initialized in a different state a_j (j ≠ i, j<100). Solid lines denote the average across j, and the dashed flanking lines indicate one standard deviation. Small values indicate that the responses to the various initial conditions a_i are roughly decorrelated. (d) Black: spontaneous fluctuations around baseline rate of a sample unit in the network. The corresponding rate distribution is shown on the right (black), and compared to the distribution obtained if the unit were not connected to the rest of the network (gray). The green line denotes the momentary population average rate, which fluctuates much less. (e) Histogram of pairwise correlations between neuronal firing rates estimated from 100 seconds of spontaneous activity. The black triangular mark indicates the mean (~0.014). (f) Excitatory (red) and inhibitory (blue) inputs taken in the same sample unit (top) or in a pair of different units (bottom), and normalized to z-scores. The corresponding Pearson correlation coefficients are indicated above each combination, and computed from 100 seconds of spontaneous dynamics. (g) Brown: lagged cross-correlogram of excitatory and inhibitory inputs to single units, each normalized to z-score (cf. (f), top row). The solid line is an average across all neurons; flanking lines denote ± 1 std. Inhibition lags behind excitation by a few milliseconds. Cross-correlating the E input into one unit with the I input into another unit (cf. (f), bottom row) yields the black curve, which is an average over 1’000 randomly chosen such pairs in the SOC.

Figure 7. Transient dynamics in a spiking SOC.

(a) The network is initialized in a mixture of its top two preferred initial states during the preparatory period. Top: raster plot of spiking activity over 200 trials for three cells (red, green, blue). Middle: temporal evolution of the trial-averaged activity of those cells (same color code), and that of the overall population activity (black). Rate traces were computed over 1'000 trials and smoothed with a Gaussian kernel (width 20 ms), to reproduce the analysis of Churchland et al., 2012. Bottom: sample voltage trace of a randomly chosen neuron. (b) Fast (black) and slow (brown) synaptic PSPs, corresponding to random and structured connections in the spiking circuit, respectively. (c) Distribution of average firing rates (top) and ISI CVs (bottom) during spontaneous activity. (d) Trial-averaged firing rate traces for a single sample cell, when the preparatory input drives the SOC into one of 27 random mixtures of its first and second preferred initial conditions. Averages were computed over 1'000 trials, and smoothed as described in (a). (e) First 200 ms of movement-related population activity, projected onto the top jPC plane. Each trajectory corresponds to a different initial condition in (d), using the same color code. See also Figure S1.

Figure 8. Spontaneous activity in spiking SOCs.

(a) Top: raster plot of spontaneous spiking activity in the SOC. Only the neurons in the first 5 subgroups (300 neurons) are shown. Bottom: momentary activity of the whole population (black), and of the second (green) and third (magenta) subgroups. Traces were smoothed using a Gaussian kernel of 20 ms width. (b) Pairwise correlations between instantaneous subgroup firing rates in the SOC, as empirically measured from a 1'000 second-long simulation (x-axis) versus theoretically predicted from a linear stochastic model (y-axis). Rate traces were first smoothed using a Gaussian kernel (20 ms width) as in (a). Distributions of pairwise correlations are shown at the top, for the SOC (black) and for a control random network with equivalent synaptic input statistics (brown; Experimental Procedures). (c) Distributions of pairwise spike correlations in the SOC (top) and in the control random network the random network (bottom), between pairs of neurons belonging to the same subgroup (blue), or to different subgroups (black). Spike trains were first convolved with a Gaussian kernel of width 100 ms. Gray curves were obtained by shuffling the ISIs, thus destroying correlations while preserving the ISI distribution. (d) Distributions of subthreshold membrane potential correlations. Colors have the same meaning as in (c). Voltage traces were cut off at the spike threshold. Gray curves were obtained by shuffling the time bins independently for each voltage trace. (e) Distributions of pairwise correlations between the slow excitatory and inhibitory currents, taken in the same cells (red) or in pairs of different cells (black). (f) Full lagged cross-correlograms between the slow excitatory and inhibitory currents, taken in the same cells (red) or in pairs of different cells. Thick lines denote averages over such E/I current pairs across the network, and thin flanking lines denotes ± 1 std. The peak at negative time lag corresponds to E currents leading I currents. See also Figure S2.