The first task considered here is a simple delayed nonmatch-to-sample problem (Figure 1). In every trial, we present two brief successive inputs to the network, with an intervening delay. Each input can take either of two values, labelled A and B respectively. The task is to determine whether the two successive inputs are identical (AA or BB), in which case the network should output −1; or different (AB or BA), in which case the network should output 1 (see Materials and methods for a detailed description).

Figure 1

Download asset Open asset

This simple task exhibits several interesting features. First, it is arguably the simplest possible flexible (that is, context-dependent) decision task: on sensing the second stimulus, the network must produce a different response depending on the identity of the first stimulus. Second, because the intervening delay is much longer than the neural time constant (τ = 30 ms, see Materials and methods), the network must maintain some memory of the first stimulus before the second stimulus arises. Third, to solve this task, some neurons in the network must necessarily possess some form of nonlinear mixed selectivity (note that the problem is in essence a delayed exclusive-or problem), a hallmark of neural activities in prefrontal cortices (Rigotti et al., 2013).

The networks consistently learn to perform the task with high accuracy. Figure 1b shows the time course of the median error over 20 training runs, each starting with a different randomly initialized network. The ‘error’ for a trial is the mean absolute difference between the output neuron’s response and the correct response, i.e. 1 or −1, over the last 200 ms of the trial (see Materials and methods). Shaded area indicates 1 st and third quartile over the 20 runs. To define a measure of successful performance, we set a criterion of 95 trials with error lower than 1 (indicating the mean response is closer to the correct than to the incorrect response, i.e. of the correct sign) over 100 successive trials (p<10⁻²⁰ under random choice, binomial test). The median time to criterion across 20 runs is 843 trials (inter-quartile range: 692–1125). Response error reliably converges towards a very low residual value.

How does the network represent and maintain traces of incoming stimuli? One possibility is that certain neurons encode stimulus identity by maintaining a stable ‘register’ value over time, such that the firing rate of certain cells directly specify stimulus identity in a relatively time-independent manner. By contrast, physiological studies suggest that neural coding during a working memory task is highly dynamic, with stimulus identity being represented by widely fluctuating patterns of neural responses, in such a way that the tuning of individual neurons significantly changes over the course of a trial (Meyers et al., 2008; Barak et al., 2010; Stokes et al., 2013). As shown in Figure 1c, trained networks do exhibit highly dynamic responses over the course of a trial. This suggests that the networks might make use of dynamic representations.

To analyze the encoding and maintenance of stimulus identity over time in the network, we used a cross-temporal classification approach (Meyers et al., 2008; Stokes et al., 2013; Dehaene and King 2016). We trained a maximum-correlation classifier to decode various task-relevant features (identity of first and second stimulus, and final response), based on whole-population activity at any given time, and then used these time-specific classifiers to try and extract the same task-relevant features at all possible points in time. This method can detect not only whether the network encodes a certain task-relevant variable, but also whether the representation of this variable changes over time (see Materials and methods).

The results in Figure 2 suggest a highly dynamic representation of stimuli by the network. For example, the identity of the first stimulus can be successfully decoded during both first and second stimulus presentations, as well as during the intervening delay, as shown by high classification accuracy values on the diagonal during this entire period (Figure 2, left panel). However, the cross-temporal classification performance between these two periods, as seen on the off-diagonal areas (for example, in the areas at 0–200 ms on one axis and 400–600 on the other, corresponding to training the classifier based on data from one stimulus presentation and testing it on data from the other stimulus presentation) is essentially at chance level (accuracy ~0.5), or even below chance (dark patches). This suggests that while the network reliably encodes information about 1st-stimulus identity across the first 800 ms of the trial, the way in which this identity is represented changes widely between successive periods within the trial. Similarly, the 2nd-stimulus identity is maintained from its onset until the beginning of the response period, but in a dynamical manner (low off-diagonal, cross-temporal accuracy between the 400–600 ms period and the 600–800 ms period, in comparison to the high diagonal accuracy over the entire 400–800 ms period).

Figure 2

Download asset Open asset

Figure 3

Download asset Open asset

Another feature of these plots is that the accuracy of stimulus identity decoding strongly decreases over the course of the ‘response’ period (low values along the diagonal for the 800–1000 ms in first and second panel of Figure 2). This suggests that the network largely stops maintaining information about the specific identity of previous stimuli, and instead encodes solely the actual response, as shown by the very strong classification accuracy in the upper-right portion of the third panel.

To test this interpretation, following (Stokes et al., 2013), we produce Multi-dimensional scaling (MDS) plots of population activity at different time points and for different stimulus conditions (Figure 3). MDS attempts to find a two-dimensional projection such that the distance between any two data points is as similar as possible to their actual distance in the full-dimensional space: nearby (distant) population states should thus produce nearby (distant) points on the MDS plot. Early in the trial, all possible stimulus identity combinations generate different, consistent trajectories, indicating stimulus-dependent encoding. By the late response period, however (1000 ms), the trajectories have essentially merged into two clusters, corresponding to the network response (‘same’ or ‘different’) and largely erasing any distinction based on the specific identity of either first or second stimulus. Thus, during the response period, the network moves from a stimulus-specific representation to a response-specific representation: the stimulus-specific response is flexibly routed to the appropriate, context-dependent response state, as previously observed in cortical activity during a flexible association task (Stokes et al., 2013).

An important aspect of cognitive control is the ability to attend selectively to specific portions of the sensory input, while ignoring the rest, in a flexible manner. Recently Mante, Sussillo and colleagues have studied the neural basis of this ability in macaque monkey prefrontal cortex (Mante et al., 2013). Monkeys were trained to report either the dominant color or the dominant motion direction of randomly-moving colored dots. Thus, the same stimulus could entail different appropriate responses depending on current context (i.e. which modality - color or motion - was relevant for this trial). Furthermore, due to the noisy stimulus, the task required selective temporal integration of the relevant sensory input. In addition to neural recordings, Mante and colleagues also trained a recurrent neural network to perform the same task, using supervised learning based on Hessian-free optimization. By analyzing the trained network, they identified mechanisms for selective integration of context-dependent inputs in a single network (Mante et al., 2013). This task was also used as an example application by Song and colleagues for their recurrent network training framework (Song et al., 2016).

We trained a network to perform the same task, using our proposed plasticity rule (see Figure 4). Our settings are deliberately similar to those described by Mante, Sussillo and colleagues. The network has two ‘sensory’ inputs (representing the two stimulus modalities of motion and color), implemented as random (Gaussian) time series, with a randomly chosen mean for each trial; the mean of each time-series thus represent the ‘value’ of the corresponding modality for this trial. In addition, two ‘context’ inputs specify which modality is relevant for each trial. The network must produce output 1 if the context-indicated sensory input has a positive mean, or −1 if it has negative mean (see Materials and methods for a detailed description).

Figure 4

Download asset Open asset

Figure 4b shows the psychometric curves of a fully-trained network, that is, the mean response as a function of stimulus value. For either modality, we show separate psychometric curves for when this modality was the relevant one and when it was irrelevant. When trials are sorted according to the value of the relevant modality, responses form a steep sigmoid curve with a relatively sharp transition between −1 and +1 centered roughly at 0. By contrast, when trials are sorted according to the value of the irrelevant modality, responses are evenly distributed across the entire range. Thus, the network accurately responds to the relevant signal, while largely ignoring the irrelevant one in each context (Compare to Figure Extended Data 2 in Mante et al., 2013). This indicates that the network has learned not only to perform temporal integration of an ambiguous, stochastic input, but also to flexibly ‘attend’ to different input streams depending on context.

How is information represented in the network over time? We use Mante and Sussillo’s orthogonal decoding procedure, which seeks to extract independent measures of how various task features (stimulus values, context, decision) are encoded in the network (see Materials and methods). Briefly, this method consists in using multiple linear regression to measure how much certain task-relevant features are being independently represented by the network at any time (see Materials and methods for a detailed description). The results are shown in Figure 5 (compare to Figures 2 and 5 in Mante et al., 2013). These trajectories plot the evolution of network information over time, for various combinations of context (relevant modality) and task features. Each trajectory represents the average population activity, at successive points in time, of all correct trials that have the same bias value for the averaging modality and resulted in the same final choice. Trajectories are colored according to the bias value of the averaging modality, ranging from −0.25 (bright red) to 0.25 (bright green). We project these averaged population trajectories along the orthogonal feature dimensions extracted by orthogonal decoding, which tells us how strongly the network encodes this particular feature at a given time. We then plot the resulting trajectories in feature dimension subspaces (‘final choice’ dimension is always used as the horizontal axis, while the vertical axis may be either of the two sensory modality dimensions).

Figure 5

Download asset Open asset

As observed in cortical recordings (Mante et al., 2013), these trajectories reveal that both the relevant and the irrelevant modality are actually represented in the network: the trajectories for varying value of either modality form an ordered progression in the corresponding ‘modality’ dimension (y-axis), even when that modality is irrelevant (bottom-left and top-right panels in Figure 5); however, only the relevant modality correlates with representation of final choice (compare panels where trajectories are separated by value of the relevant vs. irrelevant modality), in accordance with physiological observations (Mante et al., 2013). This confirms that the network has learnt to selectively integrate the context-indicated variable while discarding the irrelevant one for each trial.

In both of the previous tasks, the network output was a single response channel. However, flexible behavior often requires coordinating multiple outputs, especially during movement. To test whether our plasticity rule can produce coordinated multiplexed responses, we trained a network to control a biomechanical model of the human arm. The model is a custom modification of the one described in Saul et al. (2015) (itself an extension of Holzbaur et al., 2005) and uses the Thelen muscle model (Thelen, 2003). The model implements the human upper skeleton, with 4 degrees of freedom (three at the shoulder, one at the elbow), actuated by 16 muscles attached to the shoulder, chest, and upper and lower arm bones. Each of the 16 muscles is controlled by a specific network output cell. The task is to reach towards one of two spherical targets, located in front of the body on either side of the sagittal plane. The appropriate target ball for each trial is indicated by two input channels, set either to 1 and 0 (left-side target) or to 0 and 1 (right-side target) respectively for the entire duration of the trial (700 ms). No other inputs are provided to the system. The error at the end of each trial is measured by the absolute distance between the tip of the hand and the center of the target ball, plus a small penalty for total muscle activation over the entire trial. Note that while the target balls are symmetrically arranged with regard to the body, they are not symmetrical with regard to the right arm (which is the one we model): the right-side ball is closer than the left-side one, and thus reaching either target requires qualitatively different movements.

Results are shown in Figure 6. Initially, as expected, the untrained network performs random, aimless movements, resulting in high initial error (Figure 6b). Performance improves almost from the start of the training process, reaching a low residual error after about 3000 trials. To visualize the impact of training on the dynamics of population activity, we project the population activity over its first three principal components at successive points in time over the course of each trial, using 16 trials for either target context, both before and after training (i.e., 64 trials in total). The fully trained network correctly reaches the adequate target according to context (Figure 6c).

Figure 6

Download asset Open asset

The rule proposed here is a more biologically plausible implementation of the so-called node-perturbation rule (Fiete et al., 2006, 2007), which is itself a variant of the classical REINFORCE algorithm (Williams, 1992). Node-perturbation (like REINFORCE) consists in applying random perturbations to neural responses, then modifying the weights so as to make future responses more similar (resp. less similar) to the perturbation-induced responses, if the perturbed outputs led to a better-than-expected (resp. worse-than-expected) reward. Like the Exploratory-Hebbian (E-H) method of Legenstein, Hoerzer and Maass (Legenstein et al., 2010; Hoerzer et al., 2014), our rule extracts exploratory perturbations from outputs by subtracting a running average from ongoing neural responses. Unlike the E-H method, however, it can actually learn from sparse, delayed rewards at the end of each trial, without requiring a continuous, real-time reward signal. As we explain below, this is due to the supralinear amplification of plasticity increments.

The node-perturbation rule (Fiete et al., 2006, 2007) is a reinforcement learning rule that can train neural networks using only sparse, delayed reward signals to guide the learning. The intuitive mechanism of the rule is to apply random exploratory perturbations to neural output, then modify the weights so that future, unperturbed outputs will be more similar to this perturbed output if this perturbed output turned out to elicit a ‘good’ reward - or conversely, less similar if it produced a ‘poor’ reward (note that this is the central idea of the REINFORCE algorithm [Williams, 1992]).

The node-perturbation method can be summarized as follows (Fiete et al., 2006):

1- During each episode, apply random exploratory perturbations ξ(t) to the neuron’s output y_i(t).

2- At each synapse, compute the so-called eligibility trace e_i,j by accumulating the products of output perturbations by current input at this synapse at the time of perturbation: e_i,j = Σ_tξ _i(t) x_j(t)

3- At the end of each episode, compute the reward R for this episode.

4- Add to each synaptic weight the accumulated eligibility trace, multiplied by the net (baseline-subtracted) reward for this episode: Δw_i,j = η * e_i,j * (R - R₀).

In step 4, η is a learning rate parameter, and R₀ is the predicted reward for this episode in the absence of perturbations; thus, R-R₀ reflects whether the trajectory of outputs produced during this episode was ‘better’ or ‘worse’ than usual.

The effect of step 4 is to make future trajectories more similar to the just-experienced perturbed trajectory (in response to the same inputs), if this perturbed trajectory produced a ‘good’ reward; or less similar, if the trajectory produced a ‘bad’ reward. This is because adding ξ(t)x_j(t) to w_i,j will move future responses to the same input x_j(t) in the direction of ξ(t), and thus tend to reproduce the stochastically perturbed response. This can be easily shown analytically (notice that $\frac{\partial y\left(t\right)}{\partial w}=x\left(t\right)$, at least in the limit of a linear neuron). It also makes intuitive sense: if the perturbation was positive (resp. negative), this rule would add a positive (resp. negative) multiple of x(t) to the synaptic weight, which would make the weights more (resp. less) correlated with input x(t) and thus produce a higher (resp. lower) response to future presentations of the same input x(t).

While node perturbation was introduced by Fiete and Seung (Fiete et al., 2006), the original REINFORCE paper describes several rules that implement node-perturbation, i.e. modifying weights by the accumulated product of inputs by perturbations, multiplied by net rewards (see e.g. Eq. 11 in Williams, 1992).

A difficulty with the node-perturbation rule is that it is hard to reconcile with existing models of synaptic plasticity. It requires that each synapse maintains a distinction between the ‘actual’ inputs and the ‘perturbation’ input, and then perform learning by multiplying these two different types of inputs. The biological mechanism for such learning is not obvious. This is also in contrast with standard models of plasticity based on Hebbian learning, that is, a product between inputs and outputs rather than between two forms of input.

Legenstein, Hoerzer and Maass have proposed a method, which they call the Exploratory Hebbian (E-H) method, to implement node-perturbation in a biologically plausible, Hebbian manner, using information local to the synapse (Legenstein et al., 2010). The central idea is that if perturbations are fast and strong enough, they can be extracted from ongoing neuron output by subtracting a fast running average, which should isolate fast fluctuations in neural responses. The eligibility trace is then computed as the product of inputs by these fast fluctuations in output (which are deemed to represent mostly the exploratory perturbations), producing a Hebbian, synapse-local rule which implements the node-perturbation rule.

The E-H rule is expressed formally as follows (Legenstein et al., 2010):

Where x_j is the input at the synapse, y_i is the neuron’s output, and the overbar denotes a short-term running average (i.e. ${\displaystyle \overline{Y}}$(t + 1) = α ${\displaystyle \overline{Y}}$(t) + (1-α) y(t) where α is a constant). Thus, y_i(t) - ${\displaystyle \overline{Y}}$_i(t) extracts fast fluctuations in neural output, while R(t) -$\overline{R}$(t) extracts fast fluctuations in the reward signal.

Hoerzer and colleagues (Hoerzer et al., 2014) showed that this rule still works if the R(t)-$\overline{R}$(t) term is replaced by a less informative quantity M(t) which is 1 if R(t)>$\overline{R}$(t) and 0 otherwise. They also showed that this rule can be used to successfully train chaotic recurrent networks for complex tasks.

A difficulty with the E-H rule is that it requires a continuous, real-time reward signal R(t): at every point in time, the system must know whether it is doing better or worse than before. This negates a central advantage of reinforcement learning: the ability to learn from sparse, delayed rewards.

One reason why real-time reward signals are needed in the E-H rule is that simply subtracting a running average from ongoing activity does not reliably isolate external perturbations, due to spurious relaxation effects. To take a maximally simplified example, consider the example trace in Figure 8. This represents the output of a single neuron receiving constant inputs, with a single positive perturbation at time T = 100 (top panel). To extract fluctuations, we subtract a running average from the output (bottom panel).

Figure 8

Download asset Open asset

When the perturbation is applied, initially the output is larger than the running average, and thus the difference reflects the perturbation (‘Perturbation effect’ gray area). However, the running average then increases to include the recent perturbed outputs, and now the difference between decaying ongoing activity and running average switches to negative (‘Relaxation effect’ gray area). These negative relaxation terms will be accumulated into the Hebbian product and counteract the positive, perturbation-related initial terms. In fact, for the simple case shown in Figure 8 (constant input and a single perturbation) the total sum of all $y\left(t\right)-\overline{y}\left(t\right)$ terms over a full episode will converge to zero, and so will the eligibility trace (being a product of constant inputs by fluctuations summing to zero).

In more complex cases such as recurrent networks considered here, with time-varying inputs, the cancellation will only be partial, depending on the amount of autocorrelation in the inputs. The eligibility trace will include both the product of the ‘perturbation effect’ by the inputs at the time of the perturbation (which is what we want), and also the product of immediately subsequent, unpredictable inputs by the ‘relaxation effect’ (which is undesirable). The latter part adds uncontrollable (‘garbage’) components to the eligibility trace. Furthermore, if there is any temporal correlation in successive inputs, then these ‘garbage’ components will tend to partially cancel out the perturbation-related terms. This, we suggest, causes the inability of the E-H rule to learn without continuous, real-time reward signal (we show experiments in support of our interpretation below, in Figure 9).

Figure 9

Download asset Open asset

The Exploratory Hebbian method negates this undesirable effect by using a real-time reward signal R(t) which is also high-passed by subtracting its own running average $\overline{R}$(t). Now the same relaxation effects will occur both in the y and the R traces, and thus the product of relaxation terms in both traces results in a positive value, which will actually reinforce the perturbation-related terms rather than cancel them (assuming that inputs have any significant temporal autocorrelation).

This analysis immediately suggests an alternative solution to the problem. Notice that the perturbation-related fluctuation is large, while the subsequent countering relaxation terms are small (‘Perturbation effect’ vs. ‘Relaxation effect’ in Figure 8). Thus, if we impose a supralinear function on plasticity increments, the large perturbation-related terms will be amplified, while the small relaxation-related terms will be suppressed. This leads to the plasticity rule proposed in the present paper.

Note that this is conceptually related to recently-proposed thresholded Hebbian rules, whereby plasticity is only triggered by events in which the Hebbian product reaches a certain threshold (Soltoggio and Steil, 2013). A supralinear amplification offers a smoother amplification of larger Hebbian events, by comparison to the all-or-nothing effect of a threshold; however, the overall effect is similar: ignore small, possibly incidental correlations of input and output, but retain the larger ones, which are more likely to be informative.

Node-perturbation increases weights by a product of (net) rewards, inputs and perturbations of outputs. Why use the perturbations, rather than the raw outputs, in the product? Why not just use the standard Hebbian product of inputs by outputs (rather than perturbations), and multiply this by the reward? This is because the full Hebbian products (inputs-outputs) are dominated by reward-unrelated terms. The Hebbian products tend to simply reinforce all existing co-occurrences of inputs and outputs, regardless of how this increase would affect performance. By contrast, the perturbation-based products ensure that the weights are modified specifically to reproduce the perturbation-caused change in trajectory. The two have no reason to correlate, and the former is usually larger than the latter if perturbations are sparse.

In theory, it is possible to average out the non-reward-related component of reward-modulated Hebbian learning, by making sure that the reward signal is absolutely zero-centered separately for each trial type (that is, for each possible combination of input range and target outputs) (Frémaux et al., 2010; Sprekeler et al., 2009). This corresponds to the role of a critic (a complex reward predictor) in reinforcement learning, and is implemented in the present paper by maintaining separate values of $\overline{R}$ for each trial type. However, this procedure leaves us with the problem of considerable variance, which suffices to make learning all but impossible in complex tasks (see the specific discussion of this point in Legenstein et al. (2010): as connection weights increases, the purely-Hebbian component will tend to take over and drown the reward-related component).

This problem is compounded in recurrent networks such as the one described here, because now a single perturbation can have significant, unpredictable effects on future neural responses over the rest of the trial. For example, a single positive perturbation can have an arbitrary indirect effect on the response of the perturbed neuron, even making it eventually lower than it would have been over an extended period of time. As a result, the accumulated product of inputs and outputs would not be dominated by the (single-point) response increase caused by the perturbation, but by the (extended) subsequent decrease in response, which would actually have unpredictable effects on weight modification and future trajectories. Note that this effect (Hebbian learning futilely learning the co-activations during the perturbed trajectories, rather than the perturbations that cause the changes in trajectories) cannot be averaged out by centering the rewards, since the trajectories are what determines the rewards.

By contrast, node-perturbation only modifies weights to incorporate products of perturbations and inputs at the time of the perturbation. As a result, it will tend to specifically reproduce the perturbation itself, that is, the event that caused the complex changes in trajectory, and thus in reward.

To better illustrate the effects of the rule, we run a simple experiment to compare the weight changes computed by our rule with those prescribed by several variants and other rules under the same conditions, isolating the effects of specific elements in the rules.

The experiment involves a recurrent network similar to those described in the previous sections. The task of the network is simply to determine whether a set of inputs have a positive or negative mean. In each episode, a vector of 10 randomly chosen inputs is presented for 100 ms, then (after a delay of 100 ms) the output neuron of the network must return 1 if the inputs have positive mean, or −1 if the inputs had negative mean, over the last 100 ms of the trial (thus each trial lasts 300 ms in total). The reward for each trial is computed as the mean negative absolute difference between neuron output and target response (−1 or 1) over the last 100 ms. At each episode, we apply a single perturbation (of random sign), at a fixed point in time within the response period. We then use various plasticity rules, including the one proposed in this paper, to compute the appropriate change in weights (hereafter called the ‘gradient’), and compare the gradients obtained from different rules. Note that no learning actually occurs: weights are randomly initialized at each trial, since we are only interested in comparing the gradients obtained under various rules.

Figure 9a plots the gradient obtained by the rule presented here, against the gradient obtained by the node-perturbation method, for many episodes (each with randomly generated weights). This reveals that the two rules produce largely similar gradients, with some warping introduced by the nonlinearity. Next we apply the rule described in this paper, but without supralinear amplification: the eligibility trace is now simply the accumulated produce of inputs by fluctuations, and this is then multiplied by the net reward. Note that this is equivalent to using the Exploratory-Hebbian method, but with a single reward value for each episode, rather than a continuous real-time reward signal. The computed gradients (Figure 9b) are now uncorrelated with the node-perturbation gradients, as expected from the reasoning presented above. This illustrates the crucial role of the supralinearity. To emphasize that this difference is caused by post-perturbation effects, Figure 9c uses the same method as Figure 9b (E-H with delayed rewards, no supralinear amplification), but now only accumulates the eligibility trace product over the 10 ms that follow the perturbation time. Since neurons have a time constant of 30 ms, this results in only partial relaxation. The gradients are now largely aligned with the node-perturbation gradients. If only 1 ms was used, the gradients would be identical to node-perturbation gradients and the points would fall exactly on the median line. As the time window is increased, relaxation effects eventually overwhelm the gradient and leave a residue dominated by noise, as in Figure 9b.

As a safety check, Figure 9d confirms that the full E-H method (with a real-time, continuous reward signal) correctly recovers the node-perturbation gradients in our experimental settings.

Finally, we emphasize that the choice of supralinear function is not crucial, as long as it is indeed supralinear. In Figure 9e, using a different supralinear function (the signed square function f(x) = x.|x|) produces largely similar result to Figure 9a. By contrast, in Figure 9f, using a sub-linear function (namely the square root) largely randomizes the obtained gradients, since the (large) perturbation effects are now suppressed with regard to (smaller) non-perturbation related effects.

For most of our experiments, the model is a fully-connected continuous-time recurrent neural network of N neurons, governed by the classical RNN equations (Sompolinsky et al., 1988; Sussillo and Abbott, 2009; Jaeger, 2001; Maass et al., 2002):

where x_i is the excitation (or ‘potential’) of neuron i, r_i is its response (or ‘firing rate’ / activity), J_i,j is the connection weight from neuron j to neuron i, u_k(t) is the current value of each of the M external inputs to the network, and B_i,k is the connection weight from external input k to neuron i (τ is the relaxation time constant of neuron activation). J is initialized with weights taken from a normal distribution with mean 0 and variance g²/N, while the input weights B_k,i are fixed and taken from a uniform distribution over the [−1,1] interval. Activations x_i are initialized at the start of every trial with uniform noise in the [−0.1, 0.1] range. For the simulations reported here, N = 200 (400 for the motor control task), τ = 30 ms, and g = 1.5. Note that the latter value places the networks in the early chaotic regime, where the long-term behavior generally remains non-periodic (Sompolinsky et al., 1988).

These canonical networks have signed responses and mix positive and negative weights, which is unrealistic. However, recently, Mastrogiuseppe and Ostojic have studied the conditions under which chaotic ongoing activity can emerge in excitatory-inhibitory networks with nonnegative responses (Mastrogiuseppe and Ostojic, 2016). Using their results, we also implemented a network that produces ongoing chaotic activity with nonnegative responses, and separate populations of strictly excitatory and strictly inhibitory neurons (in accordance with Dale’s law). As we demonstrate in the last section of Results, the rule can still learn cognitive tasks in this more realistic network. In these simulations, there are 100 excitatory and 100 inhibitory neurons. The neuron response function is a nonnegative piecewise-linear function:

We set −m=−2 and M = 20. Note that the network is largely non-saturating: neural responses usually do not reach responses close to the maximum M, which is another similarity with cortical networks. The connection matrix is initialized as semi-sparse: each neuron receives connections from 50 randomly chosen excitatory neurons and 50 randomly chosen inhibitory neurons. All excitatory connections have initial weight 1. All inhibitory connections have initial weight −g_inhib = −1.2. As learning occurs, connection weights change from their initial values; however, connections from excitatory neurons are always clipped to 0 from below, while connections from inhibitory neurons are always clipped to 0 from above. The network otherwise operates as described above (in particular, the differential equation governing x_i(t) is unchanged).

In all simulations, four arbitrarily chosen neurons have a constant activation x = 1 and thus provide a bias input to other neurons. There is no separate feedback or output network. Instead, one or more neurons in the network are arbitrarily designated as the ‘output’ neurons, and their responses at any given time are used as the network’s response (these neurons are otherwise identical to all others).

Synapses between neurons are modified according to a novel form of reward-modulated Hebbian learning, which we now describe.

First, in order to produce exploratory variation in network responses across trials, each neuron in the network occasionally receives a random perturbation Δ_i(t) to its activation; these perturbations are not segregated from ‘normal’ inputs (in contrast to Fiete et al., 2006, 2007). Note that Δ_i(t) might also represent random noise, or a ‘teaching’ signal from a different area. Also, perturbations are applied in all simulations reported here, both at learning and testing/decoding time (perturbations correspond to the sharp spikes in the curves in Figure 1c, which would otherwise show only smooth curves; they also cause the spread of dots around their central values in Figure 3). In this paper, Δ_i(t) is taken from a uniform distribution within the [−0.5, 0.5] range, occurring randomly and independently for each neuron with a mean rate of 3 Hz (rates of 1 Hz or 10 Hz also give satisfactory results).

During a trial, at every time step, every synapse from neuron j to neuron i accumulates a potential Hebbian weight change (also called eligibility trace [Izhikevich, 2007]) according to the following equation:

Remember that r_j represents the output of neuron j, and thus the current input at this synapse. x_i represents the activation of neuron i and $\overline{x}$_i represents a short-term running average of x_i, and thus $x\left(t\right)-\overline{x}$ tracks the fast fluctuations of the neuron’s output. Thus this rule is essentially Hebbian, based on the product of inputs and output (fluctuations). Crucially, S is a monotonic, supralinear function of its inputs; in other words, we posit that the plasticity mechanism is dominated by large increments, and tends to suppress smaller ones. The particular choice of S is not critical, as long as it is supralinear. In this paper we simply used the cubic function S(x)=x³. Sign-preserving squaring S(x) = x|x| also gives satisfactory results; however, simply using the identity function fails to produce learning. The supralinear amplification of co-occurrences allows our learning rule to successfully learn from instantaneous deviations of activity, using only sparse, delayed rewards, without requiring a continuous, real-time reward signal (Legenstein et al., 2010; Hoerzer et al., 2014); see Discussion and Analysis.

Note that the eligibility trace for any synapse is accumulated over the course of a trial, with each new timestep adding a small increment to the synapse’s eligibility trace / potential weight change.

At the end of each trial, a certain reward R is issued to the network, based on the network’s performance for this trial as determined by the specific task. From this reward, the system computes a reward prediction error signal, as observed in physiological experiments, by subtracting the expected reward for this trial $\overline{R}$ (see below for computation of$\overline{R}$ ${\displaystyle \overline{R}}$) from the actually received reward R. This reward-prediction signal is used to modulate the eligibility trace into an actual weight change:

where η is a learning rate constant, set to 0.5 for all simulations described here.

To compute the reward prediction error signal (R-${\displaystyle \overline{R}}$), we need to estimate the expected reward in the absence of perturbation, $\overline{R}$. Following (Frémaux et al., 2010), we simply maintain a running average of recent rewards for trials of the same type (where trial type is determined by the combination of inputs). As (Frémaux et al., 2010) pointed out, it is important that separate traces should be maintained for each trial type, so as to provide an accurate estimation of the expected reward $\overline{R}$ for each trial. Thus, after the n-th trial of a given type, $\overline{R}$ is updated as follows:

Where R(n) is the reward for this trial, and $\overline{R}$(n-1) was the expected reward after the previous trial of the same type. In all simulations, α_trace = 0.33.

To stabilize learning, we clip the weight modifications for each trial to have a maximum absolute value of 10⁻⁴ (across experiments, roughly 10% of all potential weight modifications exceed this value and are clipped).

The first task considered here is a simple delayed nonmatch-to-sample problem (Figure 1). In every trial, we present two brief successive inputs to the network, with an intervening delay. Each input can take either of two values, labelled A and B respectively. The task is to determine whether the two successive inputs are identical (AA or BB), in which case the network should output −1; or different (AB or BA), in which case the network should output 1. We specify the input stimuli by using two different input channels u1 and u2; the identity of the input stimulus is determined by which channel is activated (i.e., for stimulus A, u1 = 1 and u2 = 0; for stimulus B, u1 = 0 and u2 = 1; remember that each input channel u_k is transmitted to the network by its own independent set of weights - see Model Description above). In every trial, the first stimulus is presented for 200 ms, then after a 200 ms delay the second stimulus is presented for 200 ms. Outside of input presentation periods, both input channels are set to 0. The trial goes on for an additional 400 ms, thus each trial is 1000 ms long. The network’s overall response is determined by the activity of the arbitrarily chosen output neuron over the last 200 ms of the trial (the so-called ‘response’ period). The overall error for this trial is the average absolute difference between the network’s output (that is, the activity of the output neuron) and the target response (1 or −1 depending on presented stimuli), over these last 200 ms.

For details on how the network activity was analyzed for Figures 2 and 3, see below.

This task was introduced by Mante and colleagues (Mante et al., 2013). In this study, monkeys looked at randomly-moving colored dots, in which both the value and coherence of motion direction and dot color varied from trial to trial. Monkeys had to report the dominant motion direction, or the dominant color, according to current task conditions; thus, the same stimulus could entail different appropriate responses depending on current context. Furthermore, due to the noisy stimulus, the task required temporal integration of sensory input. Importantly, the authors showed that prefrontal neurons registered inputs from both the relevant and the irrelevant modality; however, inputs from the irrelevant modality had no long-term impact on neural activities, while inputs from the relevant modality were selectively integrated over time. Thus, only information from the relevant modality contributed to the final decision.

Our settings are deliberately similar to those described by Mante, Sussillo and colleagues in their neural network implementation of the task, and Song and colleagues in their own implementation (Song et al., 2016). The network has two ‘sensory’ inputs (representing the two stimulus modalities of motion and color) and two ‘context’ inputs (which specify which of the two modalities must be attended to). The sensory inputs are noisy time-series, centered on a specific mean which indicates the ‘value’ of this input for this trial. More precisely, each of the two sensory inputs is sampled at each time step from a Gaussian variable with variance 1 and a mean, or bias, randomly set at either −0.5 or 0.5 for each trial (this is for the learning phase; for the testing phase used to generate the psychometric curves in Figure 4, the bias is varied in increments of 0.1 from −0.5 to 0.5, inclusive). The mean/bias of the Gaussian (positive or negative) represents the ‘direction’ or ‘value’ of the corresponding sensory input (left vs. right, or red vs. green). The context inputs are set to 1 and 0, or 0 and 1, to indicate the relevant modality for this trial. The goal of the network is to determine whether the sensory input in the relevant modality has positive or negative mean.

Sensory inputs are presented for the first 500 ms of the trial, followed by a 200 ms response period during which all sensory inputs are set to 0. The expected response of the network is 1 if the relevant sensory input has a positive mean, and −1 otherwise; thus the same sensory input can entail different appropriate responses depending on the context. As for the previous task, the network’s response for a trial is the firing rate of the arbitrarily chosen output cell, and the error for a trial is the average absolute difference between the firing rate of this output cell and the appropriate response for this trial (either −1 or 1) over the 200 ms response period.

For details of network analysis (Figure 5), see below.

In the delayed nonmatch-to-sample task, we used a cross-temporal classification analysis (Meyers et al., 2008; Stokes et al., 2013) to investigate how fully trained networks encode information over time (Figure 2). The interpretation of these cross-temporal decoding accuracy matrices is that they tell us not only whether the network encodes a certain task-relevant variable, but also whether it uses similar representations to encode this variable at different points in time. If we train one such classifier using data at time t in some trials, and then use it to decode population activity from data at the same time t in other trials, then decoding accuracy measures how strongly the network encodes the feature at that time point t. However, when the decoder is trained on data at time t_learn and then applied to population activity data at time t_decode, the resulting accuracy measures the stability in the network’s ‘neural code’ for this feature across both time points, i.e., how similarly the decoded feature is represented by the network across these time points. If representations are similar across both time points (that is, if the network use similar patterns to represent each possible value of the feature across both time points), then classifiers successfully trained with population activities at time t_learn should also produce accurate decoding of population activities at time t_decode. By contrast, if the network uses different representations/encoding of task features at these two time points, cross-temporal accuracy should be poor; this should be represented as ‘bottlenecks’ of high accuracy on the cross-temporal decoding plots, whereby information is high along the diagonal (i.e. the feature is indeed encoded by the network at that given time), but away-from-diagonal (cross-temporal) decoding accuracy is low. This is precisely what we observe in Figure 2.

We follow the maximal-correlation classifier approach described in Meyers et al. (2008) as closely as possible. Briefly, we want to measure how well a certain task-relevant feature (identity of first presented stimulus, or identity of second presented stimulus, or final response) can be predicted by observing network activity at time t₁, using a classifier trained on network activity at time t₂. First, we sample the activation of each neuron, every 10 ms, for each trial. This data is stored in a matrix of 100 rows and 200 columns, indicating the activities (firing rates) of all 200 neurons at each of the 100 sampling times. We first generate 80 trials (20 per possible condition, where ‘condition’ is defined as one of the four possible stimulus combination: AA, AB, BA or BB) with a trained network. The time course of neural activity will differ somewhat between successive trials, even for identical conditions, due to noise. Then we iterate the following procedure. For each of all four possible conditions, we randomly choose half the trials as ‘training’ trials, and the other half as ‘testing’ or ‘decoding’ trials. The training trials corresponding to the same category that we are trying to decode (for example, all stimuli having the same first presented stimulus) are averaged together, pointwise, for each neuron and each time point, giving a ‘prototype’ matrix of activation for each neuron at each timepoint under this category. This training data allows us to decode the category of each testing trial, at each point in time, using maximum-correlation classification, in the following way. We compute the Pearson correlation between each row of each ‘testing’ trial and each row of each ‘prototype’ trial. Each such correlation between row i of a testing trial and row j of a training category-average tells us how much the population activity at time i in the testing trial resembles the average population activity at time j for this particular category. We can then select the category for which this correlation is maximal, at each training/testing timepoint pair, as the ‘decoded’ category for each testing trial. For each testing trial, this provides a 100 × 100 matrix of decoded categories (one for each pair of training and testing timepoints). Of course, each testing trial belongs to only one category, so only one possible answer is correct, and thus we can compute another 100 × 100 matrix of binary values, indicating whether the decoded category at a given point in the decoding matrix (i.e., for any given pair of training and testing timepoints) is correct. The average of these ‘correctness matrices’, over all testing trials, provides the accuracy in cross-temporal decoding of this category for every training/testing pair of timepoints. We iterate this whole procedure 100 times and average together the resulting ‘correctness’ matrices. The resulting 100 × 100 matrix indicates at each row i and column j the proportion of times that the decoded category for population activity at timepoint j was correct, using training data from timepoint i. This is the matrix shown in each of the panels in Figure 2 (one for each of the three categories to be decoded).

For the selective integration task, we used the analysis method introduced by Mante, Sussillo and colleagues (Mante et al., 2013), and also used by Song and colleagues (Song et al., 2016) (see Figure 5). Intuitively, the purpose of this method is to estimate how much information the network encodes about different task feature (input value, context, final choice, etc.) independently from each other.

After generating multiple trials under various conditions (context -- that is, relevant modality -- and bias for each modality) with a fully trained network, we regress the activity of each neuron over the values of features of interest (context, value of each modality, and final choice) for each trial. This gives us a set of weights for each neuron, one for each feature, representing how much each feature influences the neuron’s firing rate. We then ‘switch views’ by grouping together all such weights for any given feature (200 weights - one per neuron). This in turn produces vectors in neuron population space, along which the feature is in a sense maximally represented (notice that this is quite different from, and not equivalent to, the simpler idea of simply regressing each feature over the firing rates of the neurons across trials). We then orthogonalize these vectors using QR decomposition, to ensure that these representations are as independent from each other as possible. Projecting population activity at a given time over the resulting vectors approximates the network’s current estimate of the corresponding feature at that time. For successive time slices, we average network activity vectors corresponding to the same value of bias in a certain modality, a certain attended modality, and a certain final choice. We refer the reader to Mante et al. (2013) for a complete description of the method.

We project population activity, averaged within various groups of trials, at each point in time, over these decoding axes. The trials are grouped according to final choice, value of one modality (either modality 1 or modality 2), and current context (i.e., relevant modality), and the population activity at each point in time is averaged across all trials within each group. When the resulting averages are projected over the orthogonal feature vectors, they produce trajectories, indicating the network’s encoded value for each feature, at each point in time, for trials of this group. Only correct trials are used, and thus certain combinations are impossible (for example, positive value of modality 1 bias, while attending modality 1, with a final choice of −1); this is reflected in the top-left and bottom-right panels of Figure 6, which contain half as many trajectories as the top-right and bottom-left panels.

• Thomas Miconi

• Thomas Miconi