Predictive coding of dynamical variables in balanced spiking networks

Abstract

Two observations about the cortex have puzzled neuroscientists for a long time. First, neural responses are highly variable. Second, the level of excitation and inhibition received by each neuron is tightly balanced at all times. Here, we demonstrate that both properties are necessary consequences of neural networks that represent information efficiently in their spikes. We illustrate this insight with spiking networks that represent dynamical variables. Our approach is based on two assumptions: We assume that information about dynamical variables can be read out linearly from neural spike trains, and we assume that neurons only fire a spike if that improves the representation of the dynamical variables. Based on these assumptions, we derive a network of leaky integrate-and-fire neurons that is able to implement arbitrary linear dynamical systems. We show that the membrane voltage of the neurons is equivalent to a prediction error about a common population-level signal. Among other things, our approach allows us to construct an integrator network of spiking neurons that is robust against many perturbations. Most importantly, neural variability in our networks cannot be equated to noise. Despite exhibiting the same single unit properties as widely used population code models (e.g. tuning curves, Poisson distributed spike trains), balanced networks are orders of magnitudes more reliable. Our approach suggests that spikes do matter when considering how the brain computes, and that the reliability of cortical representations could have been strongly underestimated.

Figure 1. Spike-based implementation of linear dynamical systems.

(A) Structure of the network: the neurons receive an input , scaled by feedforward weights , which is internally processed through fast and slow recurrent connections, and , to yield firing rates that can be read out by a linear decoder with weights to yield estimates of the dynamical variables, . Connections: red, excitatory; blue, inhibitory; filled circle endpoints, fast; empty diamond endpoints, slow. (B) Exemplary, effective postsynaptic potentials between neurons from two different networks. (C) Sensory integrator network for (perfect integrator). Top panel: Sensory stimulus (blue line). Before , the neurons integrate a slightly noisy version of the stimulus, , where is unit-variance Gaussian noise. At s (downward pointing arrow) all inputs to the network cease (i.e. , ). Middle panel: Raster plot of 140 model units for a given trial. Top 70 neurons have negative kernels (), and bottom 70 neurons have positive kernels (). Each dot represents a spike. Thin blue line: state . Thick red line: Network estimate . Bottom panel: Mean firing rate (over 500 presentations of identical stimuli , but with different instantiations of the sensory noise ) for the population of neurons with positive kernels (magenta) or negative kernels (green). (D) Same as C but for Hz. Parameters in A–D: , for , for , , Hz, Hz, , , (in C) and (in D). Simulation time step (Euler method) msec. The noise parameters, and , represent the standard deviation of the noise injected in each ms time step.

Figure 2. Response properties of the sensory integrator.

(A) Tuning curves to variable for the network with uniform kernels. Plain line: . Dashed line: . Parameters are as in Fig. 1 C. Tuning curves were obtained by providing a noiseless () sensory input of various strength during the first 250 ms, then measuring sustained firing in the absence of inputs during the next 1000 ms. The response shown is averaged over 500 trials. (B) Example tuning curves for the inhomogeneous network (Plain lines: all components of positive. Dashed lines: all elements of negative). Parameters are , , for , for , is a uniform distribution within , is a binomial distribution, , Hz, Hz, , ,  = 0, based on a simulation with the Euler method and time step msec. (C) Inter-spike interval distribution for a typical unit (inhomogeneous network with ). ISI distribution is measured during “persistent activity” in the absence of sensory stimulation (firing rate is constant at 5 Hz). Red lines show the prediction from a Poisson process with the same rate. (D) Mean cross-correlogram for a pair of units with the same kernel sign (inhomogeneous network). Probability of a spike in unit 1 is plotted at different delays from a spike in unit 2.

Figure 3. Response of the inhomogeneous integrator network.

Same format as in Fig. 1 C. The network is entirely deterministic (). Top panel: sensory input (blue line). Before , the sensory signal is corrupted by sensory noise with variance . Sensory input and sensory noise stop after s, at which point the network is entirely driven by its own internal and deterministic dynamics. The network is simulated twice using exactly the same initial conditions and input . Up to s, the two simulations give exactly the same spike train as represented by the dots (deterministic network with identical inputs). At s, a small perturbation is introduced in the second simulation (a single spike is delayed by 1 ms). The subsequent spike trains are completely re-shuffled by the network dynamics (First simulation: dots. Second simulation: circles). Simulation parameters are , , for , for , is a uniform distribution within , is a binomial distribution, , Hz, Hz, , , , based on a simulation with the Euler method and time step msec. Bottom panel shows the mean instantaneous firing rate for the population of neurons with positive kernels (magenta) and negative kernels (green) measured in an exponential time window with width 100 ms.

Figure 4. Membrane potential profiles for the integrator networks.

(A) Homogeneous network. Example profiles for two neurons with identical kernels. Vertical line represents a spike in the red unit, plain horizontal line represents the firing threshold, and dashed horizontal line the reset potential. Values are interpreted in mV after rescaling the membrane potential (mV and mV). These profiles are taken from the simulation shown in Fig. 1 C. (B) Inhomogeneous network. Membrane potential profiles for two neurons with strongly correlated kernels (i.e. large ) and no synaptic background noise (). These profiles are taken from the simulation shown in Fig. 3. (C) Schema explaining the distribution of spikes across neurons in a homogeneous network (see text). (D) Same two units as in (B) shown for a longer period of time.

Figure 5. Spike-based implementation of a 2-D arm forward model.

(A) Network response for a reaching arm movement. Top panel: Control variables (force exerted on the arm in and axis). Bottom panel: raster plot for a sub-population of 140 neurons. Thin lines: Real arm state ; Thick lines: network estimate . Thin and thick lines are perfectly superposed. Blue and green: positions and . Red and cyan: velocities and . (B) Tuning curve to direction for an example unit. Blue, Red and Magenta correspond respectively to arm speed of 0.2, 0.5, and 1 m/s, as represented by the inlaid schemata. (C) Tuning curves to direction (same neuron as in B) tested at 3 different arm starting position. Blue, Red and magenta correspond to arm position , and . (D) Direction tuning at 3 different arm positions for another example unit (same legend as C). (E) Schema explaining the tuning properties of model units. Dots in panels B and E represents the same arm state. Parameters in A–D: , , , normalization constraint , , , , Hz, Hz, , , .

Figure 6. Other example networks, same format as Fig. 3.

(A) Neural implementation of a “leaky differentiator”. The network tracks two dynamical variables with a state transition matrix . Top panel: command variable . (Note that is zero.) Bottom panel: network response and estimates. Thick red and purple lines: Network estimates and . Thin blue lines: Variables and . The variables and network estimates perfectly track each other, making the thin blue lines hard to see. Overlaid dots represent the corresponding output spike trains, with a different color for each neuron. (B) Neural implementation of a damped harmonic oscillator. The network tracks two dynamical variables with . Format as in A. Simulation parameters for A and B: 2-D vectors were generated by drawing each coordinate from a normal distribution and normalizing the vectors to a constant norm, so that . Other parameters were , ms, Hz, Hz, , . Dots represent spike trains, one line per neuron, shown in black to improve visibility.

Figure 7. Distinguishing spiking codes from Poisson rate codes.

(A) Example profile of total excitatory current (red) and inhibitory current (blue) in a single unit on two different time scales (time scale of the stimulus and time scale of an inter-spike interval). Currents were convolved with a 2 ms exponential time window. (B) Response of the homogeneous integrator network (same parameters as in Fig. 1 C). The input is shown in the top panel. (C) Spike trains (dots), true state (blue), and estimate (red) for a rate model with the same slow connections and input as in B. Fast connections were removed and the greedy spiking rule was replaced by a random draw from an equivalent instantaneous firing rate. Four different trials are shown (four thick red lines) to illustrate the variability in the rate model's estimate. (D) Spike trains (dots), state (blue) and estimate (red) when each spike train is recorded from a different trial of the network shown in (B). (E) Estimation error, , as a function of the number of recorded neurons, , for a spiking network with neurons. For the blue line, all neurons were recorded simultaneously, for the red line, each neuron is recorded in a different trial (red). The red line follows perfectly the prediction for independent Poisson processes. Data are from an homogeneous integrator network with and , other parameters as in Fig. 1 C. (F) Effective connectivity filters of two randomly chosen pairs in the network, as measured through a GLM analysis. (G) Consequence of suddenly inactivating half of the active neurons for the network shown in B. Blue bar: unit to inactivated (membrane potential set to ). Orange bar: units to inactivated. Other parameters as in Fig. 2 B–D,F.