Models for working memory dynamics.

A: Attractor networks generate a set of slow, non-interacting dynamical modes that store stimuli in a stable code after stimulus removal. B: Feedforward dynamics involve directed interactions over a sequence of dynamical modes, and store stimuli in a dynamic sequence of neural activity. C: Optimised networks exhibit complex, high-dimensional rotational dynamics. D: Information maintanence in each network in the presence of noise (with 10 neurons, optimised for stimulus readout at t = 10). Attractor networks rapidly lose information after stimulus offset due to accumulation of noise (top). Feedforward and optimised networks can achieve close to theoretically optimal decoding performance up to the optimised readout time. E: Stimulus encoding improves as the number of neurons (or dynamical modes) increases in feedforward and optimised networks, but not in attractor networks. F: The energetic cost (measured as the squared magnitude of the network response integrated over time) increases exponentially with the number of neurons for feedforward networks (note log scale on y axis), is independent of the number of neurons for attractor networks, and decreases with the number of neurons for optimised networks.

Optimisation of a two-neuron network for a WM task.

A: Networks were optimised to maximise the discrimination of two noisy input stimuli after a delay. B: Responses of the optimised network to each stimulus. C: Dynamics of the network at four points during optimisation: the initial network (left), the final network before feedforward dynamics emerged (middle left), the final network before rotational dynamics emerged (middle right) and the network at convergence (right). Trajectories show the mean network response for each stimulus, with solid lines up to the optimised decision time and dashed lines thereafter. Insets show the decomposition of the network into dynamical modes (Schur decomposition). D: Response distributions under each stimulus at the decision time. E: Performance of a decoder optimised for time t and tested on time t′. Dashed lines show the optimised decision time.

Numerically optimised higher-dimensional networks.

A-C: Left: Functional connectivity between dynamical modes of an optimised 6-dimensional network (A), Legendre Memory Unit (B) and feedforward network (C). In each network, the stimulus input targets a single source mode and is cascaded through the other modes. Right: The response of each mode to two input stimuli. D: Left: The inputs (top) and weight matrix (bottom) of an optimised 10-dimensional network. Right: The Schur basis of the weight matrix, which represents interactions between the dynamical modes (as in A). E: The SNR of network output at the optimised decision time, compared to that of an optimal feedforward network, a Legendre Memory Unit and an attractor network. F: Comparison of the filter applied by the network to task inputs vs that of the Bayesian ideal observer solution. The ideal observer applies a delay filter which assigns a non-zero weight only to inputs at the time of stimulus presentation while filtering out pre- and post-stimulus noise input. The LMU and and optimised networks approximate this filter using a set of oscillatory basis functions. Attractor networks can only implement constant or exponentially decaying filters.

Working memory for multiple stimuli.

A: Tuning curves of feedforward input to 4 neurons for a set of 8 stimuli. B: The inputs are constrained to a circle within a two-dimensional plane spanned by two basis vectors. C: The optimised network comprised two orthogonal, non-interacting subspaces, each of which integrated inputs along separate basis vector. D: Each plane exhibited rotational dynamics identical to that of networks optimised to discriminate two stimuli. Note that overlapping traces were shifted slightly for visualisation purposes. E: The response tuning curves changed dynamically over the delay, with two neurons switching to their anti-preferred stimulus. F: Cross-temporal pattern similarity analysis revealed dynamic coding early in the delay and stable coding throughout the rest of the delay. G: The representional geometry of the population response remained stable throughout the delay.

Rotational dynamics in monkey PFC during a spatial WM task.

A: PSTHs of two neurons for each of the 8 stimuli in the task (coloured lines). The TDR fit is shown in black. B: Cross- validation performance of the TDR model with different numbers of stimulus components (red bars). A baseline model in which the PSTH formed on the training data is used to predict the PSTH on the test data is shown for comparison (grey bar). Error bars show mean ± SEM over 12 sessions (sessions were excluded if the R2 for the baseline model was negative). Grey traces show individual sessions (note that the optimal number of subspaces was 2 for all included sessions). C: The stimulus factors f(θ), g(θ) learned from a TDR fit to simulated data under the models. D: The stimulus-loadings of a TDR fit to experimental data (a single recording session with 586 neurons). E: The population PSTHs (principal components of w1(t)) for each of the models. Note that the population PSTHs of w2(t) were identical to those of w1(t) for all models. F. As in E, but for experimental data. The second subspace w2(t) is shown in Supplementary Figure S10C. G: The alignment between the subspaces spanned by w1(t), w2(t) and w1(t′), w2(t′) for each of the models. H: Left: As in G, but for experimental data. I: The SNR of responses of the three models (computed on antipodal stimulus pairs). SNRs were computed assuming stationary state covariance, in contrast to Figure 1 which assumed a zero variance initial state. J: The cross-validated SNR of experimentally-recorded responses projected onto the TDR Principal Component subspace (mean±SEM over 4 antipodal stimulus pairs).

Behaviour of attractor (normal) and functionally-feedforward (non-normal) networks on the WM task with delta pulse stimulus input and stationary state covariance.

A: Normal networks can be characterised by a set of independent eigenmodes. Each eigenmode applies an independent filter to the projection of network input onto its left eigenvector . A downstream readout could in principle combine information from multiple such eigenmodes to obtain information about different projections of the recent history of sensory input over different timescales. B: A single eigenmode can be optimised for readout at a specific delay time td by setting its time constant as τk = 2td. C: Readout performance along a single eigenmode decays exponentially after stimulus presentation. D: Even when optimised for a given delay, the readout performance decays as (etd)−1 for attractor networks. E: Functionally-feedforward networks are better characterised by a Schur decomposition, in which each mode has a decay time constant τ and a feedforward weight onto the next mode ω (also called a homogeneous delay line). F, G: Activity propagates along the delay line, and depending on the value of ωτ, can either decay (ωτ < 1, F) or grow (ωτ > 1, G) exponentially along the delay line. H: Despite being asymptotically stable when τ > 0, activity can grow to arbitrarily large values as ω and N are increased. I-K: The performance of the functionally-feedforward network grows as ω and N are increased and τ is decreased, with SNR → ∞ when N, ω → ∞ and τ → 0. Thus, the delay line can substantially outperform normal networks (A-D), but these improvements require an exponential increase in the firing rate of neurons to achieve (H).

Bayesian deal observer solutions to WM and perceptual decision-making tasks.

A: Stimuli can be optimally discriminated based on network input by first projecting these inputs onto their linear discriminant, then convolving them with a delay filter matched to the readout time td and applying a threshold operation (not shown). B: Signal-to-noise ratio of network input vs time for the WM task with delta pulse stimulus input. C: A linear time-varying linear filter f(td, τ) = δ(td − ts − τ) achieves the optimal response SNR. D: The SNR of filter output is infinite at the optimised readout time td. The optimal time-varying filter achieves perfect performance at all times t > ts, while a static filter can achieve optimal performance only at the optimised delay td. E: The SNR of filter output vs time for static filters optimised for each readout time. F: In a perceptual decision-making task in which the stimulus is presented continuously, the SNR of network input is zero before stimulus onset and constant after stimulus onset, requiring temporal integration to achieve a high output SNR. G: The optimal time-varying filter for this task integrates the history of network input for all times t > ts, but assigns zero weight to inputs at times t < ts. H: Response SNR for the optimal time-varying filter grows linearly after stimulus onset, while an optimal static filter has bounded performance. J-M: When the stimulus appears for a fixed interval of time Tcue, the optimal filter integrates network input during the stimulus period (with an integration window Tcue) and then delays this input for a period td − ts − Tcue.

Influence of energetic penalty on learned dynamics.

A: As in Figure 2C, but showing the linear discriminant and two dynamical modes (Schur modes) at each stage of optimisation (without energy penalty). B: The mean separation and variance of responses along the readout vector at the decision time, at each iteration of the optimisation process for networks with and without a penalty on energetic cost. C: The SNR of responses for the two networks. D: The total energy of the network response (time-integrated squared norm of firing rate vector). E: The time constants of the two dynamical modes. Before the transition to rotational dynamics, the time constants are determined by the (real) eigenvalues. After the transition, they are given by the real and imaginary parts of a complex-conjugate pair of eigenvalues. F: The angle between each input mode (left eigenvector) and the input linear discriminant. After the transition to rotational dynamics the input Schur mode is shown. G: Non-normality of the dynamics matrix A (normalised between 0 and 1, see Methods).

Comparison of networks optimised using 1) SNR loss vs backpropagation through time (BPTT) with mean squared error (MSE) loss, 2) loss evaluated at a single time point vs all time points over the delay, and 3) a fixed decoder vs a time-varying decoder.

Top row: Network optimised using BPTT with MSE loss at a single time point td = 50. Second row: Optimisation with the SNR loss, weighted equally at all delays . Third row: optimisation using BPTT+MSE with a separate linear decoder fit to each time point td ∈ [0, 50]. Bottom row: Optimisation using BPTT+MSE with a single decoder fit to minimise MSE across all time points td ∈ [0, 50].

Rotational solutions are insensitive to initial network weights.

Each row shows four phases of the optimisation of a network with different random weight initialisation.

Effect of initial state covariance on optimal dynamics.

A: We derived analytical solutions for a one-dimensional integrator with recurrent weight λ receiving inputs from one of two stimuli given by u(s0, t) = −δ(t − ts) + n(t) and u(s1, t) = δ(t − ts) + n(t). The integrator is inititialised at a fixed (zero variance) initial condition at a time ts −t0 and the readout occurs at time ts +td, where ts is the stimulus time. Thus, noise accumulates for a time window t0 before stimulus onset. B: The response SNR at the readout time depends on the values of t0 and λ. C, D: The optimal value of λ is positive for t0 < td and negative for t0 > td, reflecting a transition from leaky integrator to unstable amplifier as the pre-stimulus noise integration time window t0 is decreased. E: Numerically optimised networks exhibit rotational dynamics, but the damping rate exhibits the predicted transition from stable and

Effect of length of stimulus period and initial state covariance on optimal dynamics.

Top row shows effect of changing the length of the stimulus period Tcue. The readout time was td = 50 and networks were driven by noise for a period t0 = 100 before stimulus onset (at ts = 0). All optimised networks had complex eigenvalues, and the differences in trajectories were primarily due different input durations Tcue rather than changes in flow fields. Bottom row shows effect of varying the duration of pre-stimulus noise accumulation t0 for networks with the stimulus presented continuously from t = 0 to t = td = 50. The network optimised with t0 = 0 learned a classic line attractor solution to the task, which is known to be optimal in this setting (Gold and Shadlen, 2007). In contrast, for t0 > 0 a line attractor would generate substantial pre-stimulus variability, and so is no longer optimal. In this case, networks must trade off integration during stimulus presentation against avoidance of integration of pre-stimulus noise (Supplementary Figure S2F-I). All networks optimised with t0 > 0 had complex eigenvalues, suggesting that rotational dynamics are optimal for evidence integration tasks in which pre-stimulus noise influences task performance. Note that, although the trajectories during the stimulus presentation (solid lines) use the linear part of the curved rotational trajectories, noise inputs before stimulus onset would have been integrated for a longer period of time and therefore would be rotated through the curved part of the trajectories (dashed lines), thereby reducing their influence on the decoder.

A: The signal (squared norm of vector separating two stimulus-evoked mean responses), noise (total variance of responses to either stimulus) and signal-noise-alignment (see Methods) for each network as a function of time during the delay. Vertical dashed black line shows optimised decision time td. Note that networks were initialised at stationary state, so that noise statistics do not vary during the delay. B: The performance of the optimal readout at each time following stimulus onset (grey line). Red lines show the analytically computed optimal normal and two-dimensional non-normal or rotational networks.

Rotational structure in optimised networks vs the Legendre Memory Unit (LMU).

A: Representation of the LMU in the Legendre basis. B: Equivalent representation of the LMU in the Controllable Canonical Basis. C: The rotational periods (given by imaginary parts of eigenvalues) of the LMU. Note that the eigenvalues are the same in the two bases in A and B. D: The ratio of each rotational period with that of the output rotational plane. E: The damping times (given by real parts of eigenvalues) of the LMU. F-H: As in C-E, but for the numerically optimised networks.

A: Stimulus loadings f(θ), g(θ) in TDR models fit to data simulated from each WM model (left panel, all models yield identical results), and three experimental recording sessions (right three panels). B: Principal component spectrum of wp(t) matrices of TDR models fit to simulated data (left panel) and three experimental sessions (right three panels). For the simulated data, only w1(t) is shown, but w2(t) has an identical spectrum and w0(t) = 0. For the experimental recordings, all three matrices are shown, with singular values normalised by the maximum singular value of w1(t). C: Responses to each stimulus along the top principal components of w1(t) and w2(t) for the amplification into attractor model (left) and three recording sessions. Other models are shown in Figure 5. D: The alignment of the plane spanned by w1(t), w2(t) and the plane spanned by w1(t′), w2(t′) (left: amplification into attractor model, right: data). E: Alignment of w1(t) and w1(t′) (left: amplification into attractor model, right: data).