Introduction

Working memory (WM) is a fundamental building block of intelligence (Baddeley, 1992). Storing items temporarily in WM enables animals to decouple motor actions from direct sensory input, facilitating flexible and adaptive behaviour and the learning of associations between temporally discontiguous events (Miller and Cohen, 2001). The neural basis of WM has long been thought to rely on persistent activity, with neurons firing both stably and selectively to specific memoranda throughout the delay between sensory input and motor response (Fuster and Alexander, 1971; Fuster, 1973; Funahashi et al., 1989; Goldman-Rakic, 1995). However, recent studies focusing on population-level neural representations during WM tasks have instead found highly dynamic neural activity patterns throughout the delay (Stokes et al., 2013; Stokes, 2015; Sreenivasan et al., 2014; Lundqvist et al., 2018; Stroud et al., 2024a; Daie et al., 2023). These findings suggest that the persistent activity hypothesis, and the mechanistic models inspired by it, are in need of revision.

Classic attractor models of WM dynamics, inspired by the persistent activity hypothesis, postulate that WM contents are stored in stable attractor states of population activity throughout the delay (Seung, 1996; Wang, 2001). In contrast, functionally-feedforward models propose that WM contents are stored in sequences of neural activity spanning the delay (Ganguli et al., 2008; Goldman, 2009; Harvey et al., 2012; Daie et al., 2023). However, neither attractor nor feedforward models can account for dynamic coding properties of WM circuits (Stroud et al., 2024a). While a hybrid model with transient feedforward amplification of external input into a stable attractor accounts for many aspects of dynamic coding (Stroud et al., 2023), attractor mechanisms are highly suboptimal in the presence of noise; recall precision deteriorates rapidly after stimulus offset due to noisy drift of the internal stimulus representation (Ganguli et al., 2008; Lim and Goldman, 2012; Burak and Fiete, 2012).

We asked whether the complex and dynamic activity patterns observed during WM tasks could be explained by a recurrent network architecture optimised for 1) robustness to neuronal noise and 2) energetic efficiency. In particular, noise introduces errors in WM tasks, and accurate WM performance therefore requires noise-resilient dynamics (Orhan and Ma, 2023; Stroud et al., 2024b; Burak and Fiete, 2012; Panichello et al., 2019). Moreover, there is considerable pressure to minimise the energic cost of neural circuit computations (Laughlin, 2001; Padamsey et al., 2022). We therefore developed a novel method to identify solutions to WM tasks that perform optimally in the presence of ongoing noise fluctuations and with minimal metabolic cost.

Surprisingly, the solutions that emerged from this method departed substantially from all previous models (Seung, 1996; Ganguli et al., 2008; Goldman, 2009; Stroud et al., 2023). In particular, we found that rotational dynamics, which are widely observed in a variety of non-WM tasks (Churchland et al., 2012; Aoi et al., 2020; Libby and Buschman, 2021), substantially outperformed the attractor and feedforward WM solutions in terms of both noise-robustness and energetic cost. These dynamics were integrated into high-dimensional feedforward-rotational motifs that closely resembled, but substantially outperformed, state-of-the-art State Space Models (SSMs) recently developed for analysis of time series data (Voelker et al., 2019; Gu et al., 2020, 2023). Moreover, these solutions recapitulated a wide range of functional properties of WM circuits, including dynamic coding at the single-neuron and population level (Spaak et al., 2017), and were present in large-scale recordings from prefrontal cortex of monkeys performing a classic spatial WM task (Panichello et al., 2024).

Our findings suggest that normative pressures on noise-robustness and energetic cost shape the complex dynamics of WM circuits. Given their similarity to SSMs, we suggest that the feedforward-rotational dynamics employed by WM circuits forms a high-dimensional substrate for flexible computation on temporally distributed data (Voelker et al., 2019; Gu et al., 2020, 2023). Beyond WM tasks, the method we develop to optimise noisy RNNs enables integration of the computation-through-dynamics framework with efficient coding theories (Averbeck et al., 2006; Sussillo, 2014; Vyas et al., 2020), allowing rigorous investigation of noise-robust and metabolically efficient dynamics across diverse cognitive tasks.

Results

Optimal Dynamics for Working Memory Tasks

We sought to determine the optimal dynamics for a WM task in which one of several possible input stimuli must be decoded from the response of a network after a delay, with ongoing noise inputs corrupting the stimulus representation over time (Figure 1, Figure 2A,B). For linear networks driven by Gaussian noise, the performance of an optimal decoder can be quantified as a signal-to-noise ratio (SNR, see Methods). We derived this SNR for a set of canonical linear WM models both analytically and numerically (Supplementary Figure S1, Supplementary Material). To test the optimality of these linear networks, we derived a Bayesian upper bound that holds for arbitrary nonlinear networks (Supplementary Figure S2, Supplementary Material).

Figure 1:

Figure 2:

We first considered the performance of two classic models for WM tasks - attractor and feedforward networks (Seung, 1996; Ganguli et al., 2008; Goldman, 2009). Attractor networks generate a set of slow, non-interacting dynamical modes that can be used to load and store information over the delay (Figure 1A, Supplementary Figure S1A). We found analytically that attractor networks are highly susceptible to noise (Figure 1D top, Supplementary Figure S1B-D), with the signal-to-noise ratio for discrimination of different input stimuli decaying exponentially after stimulus offset (Supplementary Material). In contrast, feedforward networks utilise a directed chain of dynamical modes to store a stimulus in an activity sequence that spans the delay (Figure 1B, Supplementary Figure S1E-G) (Ganguli et al., 2008; Goldman, 2009). We found that, in contrast to attractor networks, feedforward networks maintain a near-optimal SNR for a period of time after stimulus removal (Figure 1D middle, Supplementary Figure S1I-K). Moreover, for feedforward networks, but not attractor networks, the SNR increases as the number of neurons is increased (Figure 1E top vs middle, Supplementary Figure S1J). However, feedforward networks require strong, amplifying feedforward connectivity and an exponential increase in firing rates to achieve these SNR improvements (Figure 1F top vs middle, Supplementary Figure S1F-H). Thus, attractor networks are energetically cheap but highly susceptible to noise, while feedforward networks are energetically costly but substantially more noise-robust.

We asked whether solutions exist which outperform these hand-crafted attractor and feedforward networks in terms of noise-robustness and energetic cost. To this end, we optimised the recurrent connectivity of networks to achieve maximal response SNR while imposing a penalty on firing rates, using a novel method that performs gradient descent on these quantities directly (see Methods). Surprisingly, the optimised networks departed significantly from classic attractor and feedforward networks. Instead, they exhibited complex oscillatory dynamics distributed over a range of frequencies (Figure 1C). Remarkably, these networks simultaneously outperformed feedforward networks in terms of SNR and attractor networks in terms of energetic cost. As in feedforward networks, the SNR increased with the number of neurons in the network. However, unlike attractor or feedforward networks, their energetic cost simultaneously decreased (Figure 1F bottom). Moreover, the SNR of the optimised networks approached the theoretical upper bound for any (potentially nonlinear) network, suggesting that the assumption of linear dynamics was not a limiting factor.

Taken together, we find that networks optimised for WM performance exhibit a novel solution that is both noise-robust and energetically efficient, substantially outperforming classic attractor and feedforward models. We next sought to understand the dynamical mechanisms employed by these optimised networks.

Efficient Working Memory Encoding via Amplifying Rotational Dynamics

To gain insight into the dynamical solutions employed by the optimised networks, we considered a minimal model of the WM task comprising a network of two neurons receiving input from one of two stimuli presented as an instantaneous pulse (Figure 2A,B). The network was driven to stationary state by Gaussian noise before stimulus onset and optimised to maximise the SNR for discrimination of the two input stimuli after a delay, with or without a penalty on energetic cost (see Methods).

Three distinct phases were observed as the network performance progressively converged toward an optimised solution (Figure 2C, Supplementary Figure S3). In the first phase of optimisation, an approximate line attractor was formed (Figure 2C second panel). During the second phase there was a transition from attractor to feedforward dynamics (Figure 2C third panel). The final phase showed a transition to rotational dynamics (Figure 2C final panel). Similar dynamics emerged in networks with and without a penalty on energetic cost (Supplementary Figure S3).

We asked whether this rotational solution was a specific consequence of our choice of loss function, optimisation procedure, and task setup, or a general feature of optimal dynamics for WM tasks. To test whether these findings were sensitive to the loss function and optimisation method, we investigated the solutions learned by backpropagation-through-time (BPTT) with a squared error loss, using an analytical solution in the limit of large batch size (see Methods). Remarkably, this method produced near-identical weight trajectories to those produced by optimisation of the SNR (Supplementary Figure S4). Moreover, we showed that optimising the SNR is mathematically equivalent to optimising the probability of correct choice, suggesting that incremental learning rules based on trial-to-trial feedback should lead to similar learning trajectories (Supplementary Material). To exclude other factors that may bias the networks towards rotational solutions, we investigated the influence of 1) the range of delays the network was optimised for 2) the initial network weights 3) the initial state covariance at the time of stimulus onset 4) the length of the stimulus period relative to the delay. Rotational dynamics emerged in all settings (Supplementary Figures S4-S7), including the limiting case where the stimulus was presented continuously until the decision time, i.e. an evidence integration task (Supplementary Figure S7).

The robust emergence of rotational dynamics suggests that these solutions are more noise-resilient than attractor or feedforward mechanisms. We therefore investigated how networks integrate noisy sensory input to form an efficient output population code. We plotted the distribution of noisy network responses to each stimulus at the optimised decision time (Figure 2D). Response SNR is high if these response distributions have minimal overlap (Averbeck et al., 2006). At all stages of optimisation, recurrent integration and amplification of sensory input generated elliptical response distributions with low and high variance dimensions (Figure 2D, 2nd-4th panels). Moreover, for attractor and feedforward networks, common amplification of signal and noise caused the stimulus-coding direction to be aligned to the principal noise direction (Figure 2D, 2nd and 3rd panels). Thus, attractor and feedforward networks had strong information-limiting correlations (Moreno-Bote et al., 2014). In contrast, alignment between signal and noise was reduced in the rotational network (Figure 2D, rightmost panel). Thus, the rotational network acted to both amplify the input signal and minimise corruption of the signal by amplified noise (Supplementary Figure S8).

To determine how the optimal decoder varied over the delay in each network, we implemented a cross-temporal decoding analysis (Figure 2E). For attractor networks, information was encoded along a stable coding direction that was independent of delay time. In contrast, the feedforward and rotational networks exhibited signatures of dynamic coding (Stokes et al., 2013) - decoding performance depended on the disparity between the times the decoder was trained and tested on. Thus, rotational dynamics achieves optimally noise-robust WM performance through a dynamic population code.

Higher-Dimensional Networks Employ Feedforward-Rotational Sequences

Having characterised the optimal WM dynamics for two-dimensional networks, we next investigated the properties of higher-dimensional networks optimised to solve the task. Classic work has shown that feedforward networks are optimal for WM tasks in discrete-time networks, and can also outperform attractor networks in continuous-time networks (Ganguli et al., 2008; Goldman, 2009; Lim and Goldman, 2012). More recently, continuous-time State Space Models (SSMs) such as the Legendre Memory Unit (LMU) have been developed, which compute an online approximation of their recent input history in terms of a set of temporal basis functions and achieve state-of-the-art performance on memory tasks (Voelker et al., 2019; Gu et al., 2020, 2023). However, the optimal solution for WM tasks has not been resolved.

Numerically optimised networks exhibited a combination of feedforward and rotational dynamics. For optimised N-dimensional networks, all eigenvalues were complex, meaning that the network contained N/2 rotational planes (Figure 3A,D, Supplementary Figure S9). Although distributed across neurons, task-relevant input targeted a single dynamical mode, which elicited a feedforward cascade of activity through these rotational planes (Figure 3A,D). The LMU was based on a fundamentally different architecture, but exhibited rotational structure similar to the optimised network (Figure 3B, Supplementary Figure S9). Both the LMU and optimised networks departed substantially from classic feedforward models (Figure 3C). However, despite the similarity of the optimised networks to the LMU, the optimised networks substantially outperformed both feedforward networks and LMUs on the task (Figure 3E).

Figure 3:

To understand the relative performance of these networks, we compared them to the theoretically optimal ideal observer solution (Supplementary Figure S2). This optimal solution requires the network-readout cascade to implement a delay filter, assigning a positive weight to inputs at the time of stimulus presentation and zero weight to inputs at all other times to avoid integration of pre- and post-stimulus noise. We found that only the LMU and optimised networks could successfully approximate this optimal filter, with attractor networks weighting all times equally and feedforward networks assigning positive weights to a broad range of times before and after stimulus presentation (Figure 3F). Moreover, the optimised network better approximated this ideal observer solution than the LMU, with a more sharply peaked filter around the time of stimulus presentation that enabled reduced weighting of input noise onto the readout. This sharply peaked filter was formed as a weighted combination of oscillatory modes, which provide a more flexible basis for approximating complex input-output functions than purely decaying modes.

These findings suggest that a novel dynamical mechanism combining feedforward and rotational motifs enables optimally noise-robust WM maintenance. In particular, the optimised networks compute an online approximation of their recent input history in terms of a class of oscillatory basis functions, enabling accurate approximation of a wide range of input-output functions by varying the readout weights (Voelker et al., 2019; Gu et al., 2020, 2023).

Multiple Stimuli are Encoded in Orthogonal Rotational Subspaces

So far we have considered the optimal linear dynamics for a two-stimulus WM task (Voitov and Mrsic-Flogel, 2022; Daie et al., 2023). We next investigated WM tasks with arbitrary numbers of stimuli. As in classic delayed match to sample WM tasks, we considered a circular variable θ ∈ [0, 2π), for example the spatial location of a stimulus in the visual field or the colour of stimulus (Goldman-Rakic, 1995; Panichello et al., 2019). We optimised a minimal four-neuron network with cosine-tuned stimulus input (Figure 4A) for maximal response SNR after the delay.

Figure 4:

Both the stimulus inputs and network responses could be viewed within a four-dimensional state space, with the stimulus inputs constrained to a two-dimensional subspace (Figure 4B). Thus, each stimulus input could be described as a weighted sum of two orthogonal basis vectors (Figure 4B, coloured arrows). The optimised network comprised two orthogonal, non-interacting rotational subspaces (Figure 4C). Each subspace selected the input projection onto one of the two basis vectors and fed this into a private recurrent dimension via an amplifying rotational dynamics identical to that of networks optimised for two-stimulus discrimination (Figure 4D). Thus, each of the two orthogonal subspaces selectively maintained information about the magnitude of input along one of the two input dimensions, and a linear readout could reconstruct the stimulus by linearly weighting the responses within the two subspaces (see Supplementary Material).

Remarkably, this simple four-dimensional network exhibited a range of properties observed in prefrontal cortex and other brain areas during WM tasks. For example, a subset of neurons exhibited “switching selectivity”, whereby their tuning preference dynamically changed over the delay period towards their anti-preferred stimulus (Figure 4E) (Spaak et al., 2017; Libby and Buschman, 2021; Cavanagh et al., 2018). Moreover, cross-temporal pattern similarity analysis revealed dynamic coding early in the delay but stable coding throughout the remainder of the delay (Figure 4F) (Spaak et al., 2017; Stroud et al., 2024a). Finally, this dynamic code unfolded within a stable representational geometry, as estimated by the pairwise distances between population responses to different stimuli (Figure 4G)(Spaak et al., 2017).

Thus, a simple four-dimensional network optimised to perform a multi-stimulus WM task replicates wide range of properties of WM circuits via rotational dynamics in a set of orthogonal, non-interacting subspaces.

Prefrontal Cortex Exhibits Multi-Dimensional Rotational Dynamics in a Working Memory Task

The above findings suggest that rotational dynamics is optimal for WM tasks and accounts for features of dynamic coding found in the prefrontal cortex (PFC) and other brain regions. We next asked whether signatures of rotational dynamics could be found directly in large-scale recordings of neural circuit activity during WM tasks. To this end, we analysed neuropixel recordings from PFC of monkeys performing a classic spatial WM task (Panichello et al., 2024). On each trial of the experiment, a brief (50 ms) visual stimulus was presented at one of eight spatial locations. After a delay of 1.4-1.6 seconds, the monkeys reported the remembered stimulus via an eye movement to the corresponding location.

To uncover task-related neural dynamics in the data, we implemented a novel variant of Targeted Dimensionality Reduction (TDR) (Mante et al., 2013; Aoi et al., 2020). Specifically, we modelled activity of neuron i at time t on trial k as

where θ^(k) was the stimulus presented on trial k. The stimulus-dependent factors f(θ), g(θ) and temporal factors w_p,i(t) were fitted directly to the observed neural responses using least squares regression (Methods). The resulting model offered a compact description of the neural population responses in terms of a set of stimulus-dependent factors shared across neurons and temporal factors shared across stimuli.

To quantify the quality of model fit, we compared the post-stimulus time histograms (PSTHs) of each neuron under the TDR model to their true PSTHs. Despite being heavily underparameterised (37.5% as many TDR parameters as PSTH data points), the TDR model accurately fit the diverse PSTHs of neurons in the data (R² = 0.73 ± 0.06, mean±std over 25 recording sessions, Figure 5A). Moreover, using cross-validation to predict PSTHs on held out data, we found that a TDR model with two stimulus-dependent subspaces (f(θ) and g(θ)) was optimal (Figure 5B red bars), and performed substantially better than a baseline model in which the PSTHs formed on the training data were used to predict the PSTHs on the validation data (Figure 5B, grey bar). This suggested that the structure of the TDR model successfully captured shared stimulus- and time-dependent population structure in the data.

Figure 5:

Given that the TDR model provided an excellent fit to the data, we next sought to compare the population structure it identified to that expected under competing models for WM dynamics. We considered classic attractor models (Seung, 1996), feedforward models (Ganguli et al., 2008; Goldman, 2009), and task-optimised models with rotational dynamics (Figures 1-4). A TDR model with two stimulus subspaces and no stimulus-independent subspace (w_0,i(t) = 0) yielded an exact fit to the PSTHs generated by each WM model (R² = 1). The stimulus factors f(θ), g(θ) were identical across all WM models considered, reflecting the circular structure of the stimulus set used in the task (Figure 5C). The temporal factors w_p,i(t) evolved in low-dimensional subspaces captured by a small number of principal components (PCs) (Supplementary Figure S10B left). TDR applied to neural data revealed a similar low-dimensional encoding (Figure 5D, Supplementary Figure S10B right), albeit with an additional stimulus-independent subspace w_0,i(t) not present in classic WM models (Machens et al., 2010; Cueva et al., 2020).

To investigate the structure of population activity in the low-dimensional stimulus-dependent subspaces, we plotted the top PCs of the temporal factors w_p,i(t), which can be understood as low-dimensional “population PSTHs” (Figure 5E). The neural data exhibited highly rotational trajectories, which closely matched those of the rotational model and departed substantially from the attractor and feedforward models (Figure 5F). To quantify dynamic coding across models and data, we computed the cross-temporal alignment between the planes spanned by w₁(t), w₂(t) and w₁(t^′), w₂(t^′) (Figure 5G,H, Supplementary Figure S10D,E). When applied to the optimised rotational model and experi-mental data, this revealed classic signatures of dynamic coding, with the coding subspace at stimulus time and early delay period being orthogonal to that of the late delay. Attractor and feedforward models failed to replicate this dynamic code. While a recently proposed hybrid feedforward-attractor model accounted for the cross-temporal subspace alignment (Stroud et al., 2023), it exhibited lower-dimensional dynamics and population PSTHs that did not match those of the data (Supplementary Figure S10B-E).

Finally, we computed the response SNR for discrimination of antipodal stimulus pairs at each time point (Figure 5I,J). When computing the SNR we assumed that that the networks were driven to stationary state by noise input before stimulus onset. As a consequence, the SNR of an attractor network was low at all time points due to excessive integration of pre-stimulus noise along the slow attractor mode. In constrast, feedforward and optimised networks achieved substantially higher SNR, with an initial decrease after stimulus onset followed by a stable plateau (Figure 5I). Computing the cross-validated SNR of experimentally measured responses, projected into the six-dimensional TDR PC subspace, revealed a remarkably similar time course to that of the optimised model (Figure 5J).

Taken together, only the optimised model could account for the population PSTHs, dynamic coding properties, and SNR dynamics of the experimental data. Thus, neural activity in PFC during WM tasks is inconsistent with previous models involving attractor or feedforward dynamics. In contrast, rotational dynamics, emerging from directly optimised models, accurately accounts for the joint stimulus- and time-dependent neural population structure observed of PFC circuits during WM tasks.

Discussion

We developed a novel method for optimisation of continuous-time RNNs driven by noisy inputs to solve cognitive tasks and leveraged this framework to find the optimally noise-robust and energetically efficient dynamics for working memory (WM). The solutions we identified combine functionally-feedforward and rotational dynamics, and substantially outperform all previously proposed linear models of WM (Seung, 1996; Ganguli et al., 2008; Stroud et al., 2023; Voelker et al., 2019). Moreover, these solutions provide a parsimonious explanation of a wide range of functional properties of WM circuits, including the widely observed phenomenon of dynamic coding (Stokes et al., 2013; Stokes, 2015; Spaak et al., 2017; Stroud et al., 2024a; Spaak et al., 2017; Cavanagh et al., 2018; Libby and Buschman, 2021). Thus, we suggest that WM circuits are adapted for efficient computation.

Although our theoretical analysis was restricted to linear networks, several lines of argument suggest that our conclusions apply to nonlinear networks. First, the nonlinear dynamics of PFC during simple tasks has been shown to be well approximated by a linear dynamical system (Machens et al., 2010; Soldado-Magraner et al., 2024), and the dynamical solutions of task-trained nonlinear RNNs in WM and other tasks can often be qualitatively captured by linear approximations (Mante et al., 2013; Stroud et al., 2023). Second, we derived a non-dynamical Bayesian ideal observer analysis that places strong constraints on the optimal dynamics of any (linear or nonlinear) dynamical solution to the WM task, and showed that optimal linear networks achieve this upper bound. Third, we analysed large-scale recordings from monkey PFC during a WM task using a non-dynamical TDR analysis, and showed that the population structure was consistent with optimal rotational dynamics but not classic attractor or feedforward models (Seung, 1996; Ganguli et al., 2008; Goldman, 2009) or recently proposed hybrid feedforward-attractor models (Stroud et al., 2023).

Although rotational dynamics have not previously been considered for WM tasks, they have been observed across a wide range of brain areas and tasks, including in PFC in contextual decision-making tasks (Aoi et al., 2020) and in auditory cortex in implicit short-term memory tasks (Libby and Buschman, 2021). Moreover, high-dimsensional dynamics distributed across multiple brain regions have been observed in WM tasks (Voitov and Mrsic-Flogel, 2022; Daie et al., 2023). Indeed, although we analysed experimental data collected from to PFC, we speculate that the theoretically optimal rotational dynamics we identify may be distributed across a wide network of regions of the brain that participate in the task (Christophel et al., 2017; Voitov and Mrsic-Flogel, 2022).

While, to our knowledge, feedforward-rotational dynamics have not been previously considered in the neuroscience literature, their properties can be understood in light of State Space Models (SSMs) recently developed for machine learning applications. SSMs such as the Legendre Memory Unit, MAMBA and HiPPO are built upon linear dynamical systems which are mathematically optimised to compute an online approximation of their recent input history in terms of a set of orthogonal polynomials (Voelker et al., 2019; Gu et al., 2020, 2023). We found that networks optimised for WM tasks exhibit qualitatively similar properties to these SSMs, but substanitally outperform them in the WM tasks we consider. This suggests that WM circuits may compute an online approximation of their input history in terms of an orthogonal basis of oscillatory functions. Given that SSMs now achieve state-of-the-art performance on a plethora of time-series tasks, the high-dimensional rotational solutions implemented by WM circuits may be suitable for a wide range of complex tasks beyond simple delayed recall paradigms typically studied experimentally. In particular, arbitrary functions of an SSM’s recent input history can be reconstructed through simple linear readouts, providing a flexible substrate for memory-based computation (Gu et al., 2023). This use of flexible high-dimensional representations mirrors those used by PFC for complex cognitive tasks (Fusi et al., 2016), but extends them to the time domain. Thus, we suggest that PFC dynamically computes a high-dimensional, optimally compressed representation of its input history for general-purpose computation on temporally distributed data.

Methods

Task Setup

We considered continuous-time recurrent linear networks driven by noisy inputs:

The input comprised a stimulus at time t = 0 and ongoing Gaussian noise u_i(s, t) = ũ(t)u_i(s) + n_i(t) where . The stimulus timecourse was either a delta pulse ũ(t) = δ(t) (for Figures 2-4, S1-S6, S8) or a boxcar function ũ(t) = Θ(t) − Θ(t − T_cue) where Θ(t) is the Heaviside step function and T_cue is the duration of stimulus presentation (for Figures 1, 5, S2, S7, S10).

For all networks, the input noise was set to σ_n = 1. The inputs were modelled as u_i(s_j) = u₀ cos(ŝ_i−s_j), where the stimuli were s_j ∈ {0, π/2} for two-stimulus tasks and s_j = 2π(j − 1)/N_stim for multi-stimulus tasks, and the stimulus preferences were ŝ_i ∈ {0, π/2} or ŝ_i = 2π(j − 1)/N. We set u₀ = σ_n/ ∥u(s₁) − u(s₂)∥ for two-stimulus tasks and u₀ = σ_n/ ∥u(s₁) − u(s_N/2)∥ for multi-stimulus tasks in order to ensure that the input signal-to-noise ratio was normalised to 1. For multi-stimulus tasks, the inputs could be decomposed into basis vectors u(s) = c₁(s)b₁ + c₂(s)b₂, where b₁ = cos ŝ, b₂ = sin ŝ spanned a two-dimensional input space embedded in N-dimensional state space.

We analysed the SNR of an optimal linear readout of the stimulus, given by SNR where is the mean response (over trials) of the network to stimulus s at time t and Σ_x is the covariance of x(s, t) (which was independent of s). Note that, for the linear-Gaussian networks we considered, this linear readout is a Bayes optimal decoder of the stimulus from x(s, t) and the SNR is monotonically related to the probability of correct choice as , where Φ is the cumulative function of the standard normal distribution. Moreover, for continuous estimation tasks (rather than discrete classification tasks) the SNR is closely related to the Fisher information for estimation of the stimulus s from the response x. In such tasks, the SNR places a lower bound on the mean squared decoding error.

The SNR was computed directly, without sampling noisy trajectories of the system, using a Lyapunov equation to compute the covariance matrix and matrix exponentials to compute the mean response to each stimulus (see Supplementary Material). We also analysed the response energy of the network, defined as , which can be computed directly in a similar manner to the SNR.

Attractor and Feedforward Network Models

Hand-crafted networks were constructed via their Schur form A = QT Q^⊺ where Q is a unitary matrix and T is a lower triangular matrix. For two-stimulus tasks, Q was constructed such that Q_:,1 = b and all other columns of Q were orthogonal, with b = (u(s₂) − u(s₁))/∥u(s₂) − u(s₁)∥ the input linear discriminant. For multi-stimulus tasks, Q contained the two orthogonal basis vectors spanning the input space b₁, b₂. Specifically, we set Q_:,1 = b₁ and Q_:,N/2+1 = b₂, with the remaining columns orthogonal. In this case, T was divided into four equal-sized blocks, T_MultiStim = [T, 0; 0, T], with T equal to that of the two-stimulus model (described for each network below). Responses of networks were often plotted in the Schur basis x_Schur(t) = Q^⊺x(t).

For attractor networks, T was diagonal, with T_ii = λ_i = −1/τ_i where λ_i are the eigenvalues and τ_i are the decay time constants of the dynamical modes. The network had a slow dynamical mode τ₁ = τ_slow and fast dynamical modes τ_i>1 = τ_fast, which with the above choice of Q ensured that the input linear discriminant was loaded onto the slowest dynamical mode (Chadwick et al., 2023).

Feedforward networks were constructed as T_ij = τ ⁻¹δ_ij + ωδ_i,j+1, where τ is the time constant of each mode and ω is the feedforward weight from mode i − 1 to mode i. This ensured that the input linear discriminant was loaded onto the source mode of T.

The hybrid feedforward-attractor model in Supplementary Figure S10 differed from the feedforward model in two ways. First, we included only one feedforward connection. Second, this feedforward connection fed from a rapidly decaying mode into a slow attractor mode (Stroud et al., 2023). In particular, we set T₂₂ = −1/τ_slow, with T_ii = −1/τ_fast for i ≠ 2. The feedforward connection was set as T₂₁ = ω.

For Figure 1, we considered a two-stimulus task. Networks had a fixed initial state x(s_i, 0) = 0 (i.e., with zero initial state covariance) and a stimulus was presented for T_cue = 1 (in arbitrary units of time). The SNR for stimulus discrimination was calculated at the decision time t_d = 10. Attractor networks had τ_slow = 10⁴. Feedforward networks had ω = 5 and τ selected to maximise SNR.

For Figure 3, we considered a two-stimulus task in which networks were driven to stationary state by noise input before stimulus onset and received a delta pulse input, with t_d = 50. The SNR of attractor and feedforward networks was computed analytically using their optimal parameters (see Supplementary Material). However, because the optimal feedforward network has infinitely strong feedforward weights, we instead plotted the trajectories of this network with τ = 6.25, ω = 0.5 (where τ was selected as the optimal value in the limit ω → ∞).

For Figures 5 and S10, we considered an eight-stimulus task, with T_cue = 1 and t_d = 10. Attractor networks had N = 6 neurons, with slow time constants τ_slow = 100 and fast time constants τ_fast = 1. Feedforward networks had N = 12 neurons with τ = 5 and ω = 0.5. Hybrid feedforward-attractor models had N = 6 neurons with τ_fast = 1, τ_slow = 100 and ω = 1. Note that the results we present are independent of N for attractor and hybrid feedforward-attractor models (if N ≥ 2 and N ≥ 4 respectively), but not feedforward models.

Optimised Networks

We minimised the following loss function:

where the first term quantifies the decoder performance L(t) = Σ _i>j 1/SNR(s_i, s_j, t) weighted over delays by l(t) and the second term quantifies the energy of responses weighted by the regularisation parameter β. The weighting function was set as l(t) = δ(t − t_d) (where t_d is the decision time), with the exception of Supplementary Figure S4 where we set l(t) = 1 for 0 ≤ t ≤ t_d and 0 otherwise. We set β = 10⁻⁶ for Figures 1 and S3 and β = 0 otherwise. For delta input tasks we set t_d = 50, and for finite cue tasks, we set T_cue = 1, t_d = 10 (with the exception of Supplementary Figure S7, where T_cue was systematically varied).

We computed the gradient analytically and performed gradient descent numerically with a fixed step size η that was set manually for each network. Network weights were initialised as A_ij = δ_ij0.05(1 − 0.005i), except for Supplementary Figure S5 where A was intialised randomly using the eigendecomposition A = V DV ⁻¹ with V = [cos θ₁, cos θ₂; sin θ₁, sin θ₂] with θ_i ∼ Unif[−π, π] and D = [λ₁, 0; 0, λ₂] with λ_i ∼ Unif[−0.2, 0], i.e. with random (real) eigenvectors and eigenvalues.

To transform the optimised networks from neuron space to the space of dynamical modes, we performed a real Schur decomposition on the weight matrix A_Schur = Q^⊺AQ, u_Schur = Q^⊺u where Q is an orthogonal matrix and A_Schur is a block upper-triangular matrix. Given that the Schur decomposition is non-unique, we exploited a set of invariances to transform the Schur decomposition into a standardised form (see Supplementary Material for details).

Analysis of Learning Dynamics

We selected networks at four stages of optimisation to plot the flow fields, Schur modes, and responses:

1) the initial unoptimised network 2) the initial attractor network, before non-normality had increased 3) the final network before eigenvalues became complex 4) the network at convergence.

For each iteration of optimisation, we computed the SNR (defined above) and the projection of the mean and covariance of responses onto the unit length readout vector , given by and w^⊺Σ_xw respectively (e.g., Supplementary Figure S3B,C).

The decay time constants and rotational periods of dynamical modes were computed as as τ = −1/ℛ(λ) and ν = 2π/ℐ(λ) respectively, where λ_i = ℛ(λ) ± iℐ(λ) were the eigenvalues of A. The alignment of modes with the input discriminant b were computed using either eigendecomposition or Schur decomposition. For two-neuron networks with real eigenvalues, we computed the dot product between each unit length left eigenvector of A (defined as A^⊺v_i = λv_i) and the input linear discriminant b. For higher-dimensional networks or networks with complex eigenvalues, we computed the dot product between the Schur vectors (columns of Q) and b. The non-normality of the dynamics was computed as , where 0 ≤ H ≤ 1 and H = 0 for a normal (attractor) network.

To optimise networks via backpropagation through time, we used analytically derived gradients in the limit of infinite batch size (Supplementary Figure S4, see Supplementary Material for derivation). We used a squared error loss ℒ = Σ_s⟨∥w · x(s, t) + c − s∥²⟩ (where the expecation was taken over trials) and optimised A using gradient descent with a fixed step size η. The readout parameters w, c were obtained analytically at each iteration by setting and .

Analysis of Dynamic Coding

In Figure 2D, we plotted ellipses of one standard deviation of the response distribution for networks at each of the four stages of optimisation we considered.

For Figure 2E, we performed a cross-temporal decoding analysis on networks at different stages of optimisation using standard signal detection theory measures. Labelling s₁ = 0, s₂ = 1, the optimal decoder uses a binarised linear readout ŝ(c, w, t) = Θ(w(t)^⊺x(s, t) + c(t)), with and . Using an arbitrary (suboptimal) w and c results in the probability of correct choice for each stimulus and where Φ is the cumulative distribution function of the standard normal distribution. Assuming p(s₁) = p(s₂), the overall probability of correct choice is then . Using these formulae, we computed the probability of correct choice using the optimal decoder for time t on responses at time t^′.

For Supplementary Figure S8, the signal, noise, and signal-noise-alignment were computed for networks at each stage of optimisation as signal , noise = Trace(Σ_x), signal. The SNR was computed for optimised networks numerically (using a Lyapunov equation to compute the covariance) and the optimal SNRs of the attractor, feedforward and rotational networks were plotted using analytically-derived expressions (see Supplementary Material).

To replicate the analyses performed on experimental data by (Spaak et al., 2017), we computed the cross-temporal pattern similarity (the result was independent of the stimulus s for the models we considered). In addition, we computed the representational dissimilarity matrices R_ij = SNR(s_i, s_j, t) (Kriegeskorte and Wei, 2021).

Legendre Memory Unit

The Legendre Memory Unit is governed by the equations

where u(t) is a scalar input time series, τ is a global network time constant, A_ij are the recurrent network weights, b_i is an input loading vector and w_i are the readout weights (Voelker et al., 2019; Voelker, 2019). Classic results from control theory provide methods to optimise the parameters of the above system such that the output y(t) is related to the input u(t) by a desired transfer function f, i.e. y ≈ (u ⋆ f)(t) where ⋆ is convolution (Partington, 2004; Brockett, 2015). Voelker (2019) considered the case of a delay transfer function , which we show in the Supplementary Material is a critical component of the Bayesian ideal observer solution to the workingmemory task we consider (Supplementary Figure S2). In approximating the delay transfer function using the above SSM, it was further shown that the solution can be understood through the lens of function approximation theory - each element of the network state x_i(t) in the optimised system computes the coefficient of a Legendre polynomial approximation to the input time series u(t^′) over an interval t − τ ≤ t^′≤ t, and subsequent work has shown how this approach can be extended to a wide range of approximating function classes (Gu et al., 2020, 2023). As a consequence, it is possible to reconstruct arbitrary functions of the input time series over this interval through an appropriate choice of readout weights w_i, as described below.

The parameters of the LMU are found by forming a Padé approximation of the desired transfer function in the Laplace domain (Partington, 2004; Brockett, 2015; Voelker et al., 2019). The resulting system can be written in multiple equivalent bases, all of which implement an identical transfer function. In the standard Legendre basis of the LMU, the recurrent weight matrix is given by:

and the input vector is given by b_i = (2i)(−1)ⁱ⁻¹. In this basis, the activations x_i(t) represent the delayed inputs as where is the ith shifted Legendre polynomial defined on the interval (Voelker et al., 2019). As a consequence, an arbitrary filtering of the input time series can be approximately recovered by setting the readout weights as where .

The controllable canonical form provides an alternative representation of the LMU, where the weight matrix is given by A_i1 = −ν₁, A_{i,i − 1} = ν_i, with A_ij = 0 otherwise, and the inputs are given by b_i = [ν₁, 0, …, 0]^⊺ (where ν_i = (N + i −1)(N − i + 1)/i for i = 1…N). In this basis, the inputs target only the first neuron of the network, consistent with the Schur basis for the optimised and hand-crafted networks. We therefore presented results in controllable canonical form in the main text, and compared to results in the Legendre basis in the Supplementary Material. However, we note that the two forms of the LMU are fully equivalent, with identical eigenvalues and SNR.

We simulated responses of the LMU for varying numbers of units N and varying time constant τ. Inputs comprised a delta pulse with ongoing Gaussian noise u(s, t) = δ(t)u(s) + n(t). We set n(t) ∼ N(0, 1) and u(s) = ±u₀, with u₀ chosen to normalise the input SNR to 1. We computed the mean and stationary state covariance of responses and the . For each N, we chose the value of τ that maximised this SNR. We note that the task faced by the LMU is simpler than that faced by the other models we considered; for the LMU, input noise targets only the first unit, since b ∝ [1, 0, …, 0]^⊺. Indeed, when simulating LMUs with noise targetting all units (as was the case for all other networks), performance was substantially degraded and decreased with N (not shown).

Analysis of Electrophysiological Data

We analysed publicly available data from (Panichello et al., 2024). Stimuli were presented for 50 ms, followed by a delay of 1400-1600 ms after which monkeys were cued to report the remembered stimulus. We analysed data during the period from stimulus onset (0 ms) to the end of the shortest delay (1400 ms). As preprocessing steps, the raw spiking data were converted into spike counts within 25 ms bins using a sliding window of step size 5 ms. We discarded error trials and subtracted the mean spike count of each neuron (averaged over both time bins and trials). All analyses were applied to individual experimental sessions and repeated over multiple sessions to verify consistency. For Figure 5F,H,J, we used session 22-10-21. For Supplementary Figure S10, we used sessions 22-10-24, 22-10-21 and 22-10-19.

We fitted the targeted dimensionality reduction (TDR) model

to the mean-subtracted spike counts using alternating least squares, iteratively solving for w_p,I and f_p until convergence. We set f₀(θ) = 1 for the stimulus-independent subspace. In this notation, f₁(θ) = f(θ), f₂(θ) = g(θ), f₀(θ) = 1 recovers the P = 2 model described in the main text. Given the 8 stimuli in the experimental task, the f_p(θ) parameters were given by 8-dimensional vectors f_p = [f_p(θ₁), …, f_p(θ₈)]^⊺ that were inferred directly from the data. We initialised f_p(θ) ∼ N(0, 1) (for p > 0) and began the optimisation algorithm by updating w_p,i(t). For N_stim = 8 stimuli, T time points and N neurons in a given experimental session, the model had (P + 1)NT + PN_stim parameters, with NT parameters for each of the w_p,i(t) components and N_stim parameters for each of the f_p components. In contrast, the PSTHs of the data contained N_stimNT datapoints. This reduction in parameters reflects the shared structure of the TDR model across stimuli, which incorporates the assumption that the population PSTHs for different stimuli should have a similar temporal profile.

For cross-validation, we split the data into two subsets of randomly selected trials of equal size. The TDR model was fit to each split and used to predict the remaining split. We reported the cross-validated coefficient of determination

where is the predicted PSTH formed by the TDR fit to the train-ing data, is the PSTH formed on the test data and . Given that the TDR model was underparameterised, we also reported this measure without cross-validation (where were formed on the full dataset). We repeated this cross-validation for 25 random two-way splits of the data (with random TDR parameter initialisations) and averaged R² both over splits within each recording session (grey traces in Figure 5B) and over recording sessions (red bars in Figure 5B). We compared this to a PSTH only model in which the predicted PSTH was .Note that a TDR model with P = N_stim and w_0,i(t) = 0 is equivalent to this PSTH only model.

For application of TDR to classic WM models, we simulated the PSTHs (without noise) of each model under each stimulus and fit a TDR model to this set of PSTHs. We observed that a TDR model with no time component w_0,i was sufficient to fit these simulated PSTHs exactly, and therefore omitted this term.

Having fit the TDR model to data, we computed the Principal Component (PC) subspaces by taking the singular value decomposition of the matrix (W_p)_i,t = w_p,i(t). Writing W_p = USV ^⊺, where S is a diagonal matrix containing the singular values in descending order, the columns of U contained the N-dimensional basis vectors in neural state space and the rows of SV ^⊺ contained the population PSTHs. Before computing these subspaces, we first exploited several invariances of the model to transform the fitted parameters w_p,i(t) and f_p(θ) into a standard form (Supplementary Material, cf Aoi et al. (2020)). We plotted these population PSTHs within the top three PCs of W₁ and W₂, and plotted the singular value spectra (diagonal of S) to assess the dimensionality of each subspace p ∈ {0, 1, 2}.

Subspace alignments were computed using random two-way splits over trials, as for the cross-validation analysis. Given the vectors obtained from the TDR model fit to split s, we computed the subspace angle between the plane spanned by and the plane spanned by using the Matlab function subspace. We plotted the cosine of this subspace angle, averaged over 25 random two-way splits of the data, as a function of t and t^′. The cosine of the angle between and was computed in a similar fashion.

To compute the SNR of the data in the TDR Principal Component Subspace, we used a similar two-fold cross-validation split as above. First, a TDR model was fitted to the full dataset (without cross-validation). Given this fit, we formed the matrix W_concat = [W₁, W₂] to obtain an N × 2T matrix describing the set of neural activity patterns generated over stimuli and time points. We per-formed a singular value decomposition W_concat = USV ^⊺ and projected the binned spiking data onto the top six PCs as . We then computed the mean and covariance of x^(k)(t) over trials separately within each split of the data, conditioned on stimulus and time bin. This gave the set of linear discriminant vectors computed on split s,. The cross-validated SNR was computed by applying the de-coder trained on one split of the data to the responses on the second split, i.e. . We then averaged the resulting SNR over 100 ran-dom two-way splits.

Data availability

The current manuscript involves only computational modelling and previously published data. Modelling code are available on GitHub and previously published data are available on Dryad (http://doi.org/10.5061/dryad.kkwh70sct).

Supplementary Material

A Network Model

We considered a linear dynamical system governing the evolution of the network state x(t) via the differential equation

which had as parameters the connectivity matrix A and the stimulus-dependent input u(s, t) where the stimulus was labelled s. The input

consisted of a stimulus-dependent vector u(s) indicating the direction in state space along which the network was driven by the stimulus s at the time of stimulus presentation t_s, and temporally uncorrelated Gaussian white noise n(t) ∼𝒩(0, Σ_n), ⟨n(t)n^⊺(t^′)⟩ = Σ_nδ(t − t^′). We also considered extensions to temporally extended inputs u(s, t) = u(s)ũ(t) + n(t) where ũ(t) was a boxcar function (as in Figure 1).

Given the initial condition x₀ = x(t₀), Equation 1 has the general solution for t ≥ t₀

Inserting Equation 2 into Equation 3 gives (for t, t_s ≥ t₀)

where Θ(t) is the Heaviside step function. Using the shorthand for the stimulus-conditioned and marginal means of a variable a as and , the stimulus-conditioned mean and covariance of x(t) are

and

where we assumed that the initial condition is distributed as and is independent of the subsequent input noise n ∼ 𝒩(0, Σ_n). The terms and capture the response covariance generated by propagating the initial state and subsequent noise input through the recurrent dynamics, respectively. The response at time t under each stimulus is normally distributed over trials, with , and the covariance Σ_x(t) is independent of the stimulus s.

B Loss Functions and Their Derivatives

We consider the case where one of M possible input stimuli s ∈ {0, …,M −1} is presented, and optimised the performance of an optimal decoder of s from the network response x(t). Since the response statistics are Gaussian with stimulus-independent covariance, the performance of an optimal decoder when M = 2 is given by a linear readout ŝ = Θ(w_LD · x(t) + c), where is the linear discriminant vector and c is a bias term (Chadwick et al., 2023). The performance of an arbitrary linear decoder w can be described by a signal-to-noise ratio given by:

In particular, the probability of correct choice is , where is the cumulative distribution function of the standard normal distribution (Stanislaw and Todorov,1999). Although we consider discrete classification tasks, the SNR also places a lower bound on the decoding error for continuous estimation tasks. In particular, the Fisher information is given by , and the Cramer-Rao lower bound is ⟨(s − ŝ)²⟩ ≥ 1/I_F (s).

For two-stimulus tasks, we therefore defined the loss function , where and l(t) is an arbitrary weighting function. We chose to minimise the reciprocal of the SNR rather than directly maximising the SNR, as this avoids solutions where a low SNR at time point t is compensated by a high SNR at t^′.

The derivative of this loss function with respect to the weights is

where M · N = Trace(MN^⊺) = Σ_i,j M_ijN_ij denotes a sum over elements. Given the derivatives and , the above equation allows for a direct computation of the gradient of the loss function with respect to A. We describe how these derivatives can be computed in Section C.

For the general setting of networks with N neurons, M stimuli, and a penalty on response energy, we optimised the loss

where and .

C Derivatives of the Network Statistics with Respect to the Weights

To perform gradient descent on this loss function, we computed the derivatives of the mean , covariance Σ_x(t) and energy E with respect to the weight matrix A.

Derivative of mean response

The derivative of the mean response (Equation 4) is

The loss functions we consider depend on the difference in mean response to the two stimuli , which has derivative:

The derivative of the matrix exponential can be computed numerically (Najfield and Havel, 1995; Al-Mohy and Higham, 2008), and a solver is available via the expm_frechet function in SciPy’s Linear Algebra library.

Derivative of response covariance

The derivative of the covariance (Equation 5) requires the solution of two Lyapunov equations, with numerical solvers widely available. The first Lyapunov equation is the result of integrating by parts

where . The second Lyapunov equation results from differentiating the first with respect to A_ij to obtain

where , with and . Since this is another Lyapunov equation it can be solved in terms of A and Ψ to obtain . The derivative of the network covariance can then be evaluated by handling as described here and using a numerical solver for the derivative of the matrix exponential, as described above for the derivative of the mean (assuming that the initial state covariance is independent of the weights A).

Limiting cases for response covariance

Two special cases allowed us to neglect the role of the initial condition-driven covariance the scenario where the network was initialised at stationary state with respect to the input noise Σ_n, which can be found by setting t₀ → −∞, so that Φ = Σ_n and , and the scenario where , which corresponds to a fixed (zero-variance) initial condition. We focused on these special cases in order to reduce the space of free parameters in the model and to gain insight into the behaviour of the model in these limiting situations.

Derivative of the network energy term

The energy E and its gradient could be computed directly without numerical integration using , with the inte-gral and its derivative having direct solutions via Lyapunov equations, as was the case for the network covariance.

D Relationship to Backpropagation and Other Learning Rules

The loss function and optimisation procedure described above are based on the trial-averaged statistics of the network response. How does this relate to learning rules that operate based on trial-to-trial feedback, such as reward-based Hebbian learning, backpropagation through time, or reinforcement learning? To address this question, we show that 1) in the limit of slow learning rate, any learning rule that maximises the fraction of correct choices based on trial-by-trial feedback will converge to the same weight updates as our method 2) in the limit of large batch size, backpropagation through time on a squared error loss function produces similar weight updates to our method and converges to the same solution.

Relationship to learning from trial-by-trial feedback

The objective we optimise L = 1/SNR is directly related to the probability of a correct decision p_correct using an optimal decoder, since where 0.5 ≤ Φ(x) ≤ 1 is monotonically increasing (as described above). Thus, and . The factor depends only on the SNR, so that the two gradients differ only up to a scaling factor, and therefore produce identical weight trajectories. In the limit of a small learning rate, any learning rule that optimises the fraction of correct choices using trial-by-trial feedback will converge to the gradient , and will therefore follow the same learning dynamics as our method.

Relationship to backpropagation through time

Methods which do not explicitly optimise the probability of correct choice may nevertheless yield very similar learning dynamics to those produced by our method. To illustrate this finding, we derive analytical expressions for the updates generated via backpropagation through time (BPTT) and show that numerical implementation of these updates generates near-identical weight trajectories to those produced by our method. Specifically, we consider BPTT on the continuous-time linear recurrent neural network with u(s, t) = u(s)δ(t − t_s) + n(t). On each trial, we consider a trajectory from initial time point t₀ to decision time t_d with a loss L that depends only on the final time point x(t_d). The continuous-time BPTT updates are then given by

We consider the squared error loss function L(t_d) = (w · x(t_d) + c − s)² = (ŝ − s)². The terms in Equation (7) are the . Therefore

We next consider the stationary state limit t₀ → −∞, with

Inserting the expressions derived above for and we obtain

where .Using ⟨n(t)⟩ = 0 and ⟨n(t)n^⊺(t^′)⟩ = Σ_n δ(t− t^′), we then compute the expected gradient over all possible realisations of the noise term n (for fixed s):

Each of the above integrals I₁, I₂, I₃ can be reformulated as a Lyapunov or Sylvester equation, with numerically efficient and stable solvers available. In particular, I₂ = lyap(A^⊺, ww^⊺) and I₃ = lyap(A, Σ_n). Integrating I₁ by parts gives , which is a Sylvester equa-tion of the form PI₁ + I₁Q = C. However, this Sylvester equation does not have a unique solution (Sylvester equations have a unique solution if P and −Q do not share a common eigenvalue. Here, however, P = −Q, violating this condition). Instead, a unique solution for I₁ can be found by taking the eigendecomposition , which gives

where Q = (V ^⊺wu^⊺V ^−⊺) ⊙ S, with and (for i ≠ j), with ⊙ the Hadamard (or element-wise) product. Thus, each update of A during the BPTT algorithm can be achieved by solving two Lyapunov equations (for I₂ and I₃) and performing an eigendecomposition on A to compute I₁.

To complete the BPTT algorithm, we require updates for the readout weights w and threshold c. These updates are given as follows:

These readout parameters could be updated via gradient descent alongside the recurrent weights A. However, if we assume that the updates of the recurrent weights are slow relative to the readout weights, we can set

where , and we included a sum over time points to handle the case of a fixed (time-independent) readout over a range of delays. Combining these two expressions gives:

with w obtained by inserting this expression into Equation (12). Taken together, the full BPTT algorithm involves iteratively using Equations (12-13) to update readout weights followed by Equations (8-10) to update the recurrent weights (summed over stimuli and time points). An implementation of this algorithm is shown in Supplementary Figure S4.

E Optimal Solution to Task: An Ideal Observer Analysis

We consider an ideal observer of the network input u(s, t) = u(s)ũ(t) + n(t) tasked with reporting the stimulus identity s ∈ {0, 1} at time t_d. Assuming that the input noise is Gaussian, stimulus-independent, and temporally uncorrelated, the optimal solution involves two steps: 1) projecting the inputs onto the linear discriminant vector ) applying a temporal filter f that selects and delays inputs over time. Here, we formalise this solution mathematically, showing that the optimal solution is a form of matched filtering (Turin, 1960).

We consider a linear projection and temporal filtering of the input time series u(s, t ≤ t_d) onto a vector w and filter with

In particular, the above SNR is maximised when and f(t, t_d) ∝ ũ(t_d − t). These filtering and projection operations can be understood as linear discriminant analysis over both neurons and time, and can also be derived as the maximum likelihood estimate of s given the full input time series up to time t_d, i.e. by maximising p(u(s, t ≤ t_d) s), with . Thus, Equation (14) quantifies the performance of a Bayes optimal ideal observer of the network input (assuming equal prior likelihood of each stimulus).

The above solution requires a filter that depends on the optimised readout time, also known as a time-varying linear filter. In contrast, a fixed filter optimised for a single readout time t_d will have a suboptimal output SNR for t ≠ t_d

To illustrate these solutions, we provide examples from three tasks: first, a working memory task with instantaneous inputs, ũ_δ(t) = δ(t − t_s); second, an evidence integration task where ũ_Θ(t) = Θ(t − t_s); third, a working memory task where the stimulus appears for a fixed interval of time, modelled as a “boxcar” function ũ_∏(t) = Θ(t − t_s) − Θ(t − (t_s + T)). The optimal time-varying filters for these tasks are f_δ(t, t_d) = δ(t_d − t_s − t), f_Θ(t, t_d) = Θ(t_d − t_s − t), f_∏(t, t_d) = Θ(t_d − t_s − t) − Θ(t_d − t_s − T − t).

Following Equation (14), these filters give rise to the following SNRs:

These solutions are illustrated in Supplementary Figure S2, where we also show the performance of fixed filters optimised for a single readout time (Equation (15)).

F Analysis of Linear Dynamical System Task Performance

We next sought to determine the performance of linear network models on the working memory task. The above ideal observer analysis provides an upper bound on the SNR of network responses, and suggests that networks could be suboptimal due to two factors: 1) a failure to correctly combine inputs to different neurons via a suboptimal projection w, or 2) a failure to correctly combine inputs across time via the wrong temporal filtering f.

We show that attractor networks (those with orthogonal eigenvectors, also called normal networks) perform poorly on working memory tasks, with their SNR decaying at best as . We then show that (non-normal) functionally-feedforward networks can approximate the ideal observer solution to an arbitrary level of precision, but require exponential increases in firing rates to achieve this. Finally, we show that networks with rotational dynamics achieve higher SNR than both attractor and feedforward networks.

In the following we consider networks receiving delta pulse inputs (we set t_s = 0 for convenience). We consider the signal-to-noise ratio of network responses along a readout vector w defined as . We will often consider the optimal readout with .

F.1 SNR Along Individual Dynamical Modes

Previous work has suggested that alignment of dynamical modes with the input linear discriminant provides a neural mechanism for working memory and evidence integration tasks (Seung, 1996; Mante et al., 2013; Chadwick et al., 2023; Stroud et al., 2023; Pagan et al., 2024). We therefore first investigated the SNR along individual dynamical modes of a network.

To decompose the dynamics into individual modes, we use the eigendecomposition , where are left and right eigenvectors chosen such that ∥v∥ ⁱ = 1. Given Equations (4) and (5), the mean and variance of network output projected onto a left eigenvector are (for t ≥ t_s = 0 and t₀ ≤ t_s = 0)

Thus,

where and .

Optimal dynamics at stationary state

If the network is initialised at stationary state with respect to the input noise (i.e. if t₀ → −∞), the response SNR is

This solution is optimised when the left eigenvector is aligned to the input discriminant, i.e. . The SNR depends non-monotonically on the integration time constant τ_k = −1/λ_k, and has an optimal time constant for readout time t_d given by τ_k = 2t_d. Thus, a single eigenmode is optimised for a readout at time t_d if its left eigenvector is aligned to the input discriminant and its eigenvalue is . However, this optimal SNR decreases as a function of the optimised readout time according to SNR(b, t_d) = (et_d)⁻¹SNR_input(b). This is in contrast to the optimal ideal observer solution, whose performance is independent of the readout time provided t_d ≥ t_s.

Optimal dynamics for fixed initial state

If the network is initialised at a fixed (i.e. zero-variance) initial condition x₀ at time t₀ ≤ t_s = 0, the response SNR is:

Plotting this equation numerically as a function of λ_k = −1/τ_k for different values of t₀ shows that stable dynamics (λ_k < 0) are optimal when |t₀| > t, neutrally stable dynamics (λ_k = 0) are optimal when |t₀| = t, and unstable dynamics (λ_k > 0) are optimal when |t₀| < t (Supplementary Figure S6B-D). The optimal λ_k scales monotonically with t₀, respectively approaching infinitely strong amplification λ_k → ∞ as t₀ → 0 and the stationary state solution as t₀ → −∞. The amplifying dynamics for |t₀| < t act to amplify the signal throughout the delay, while the decaying dynamics for |t₀| > t are required to avoid amplification of pre-stimulus noise input above the signal.

These solutions are highly suboptimal with respect to the ideal observer solution. This suboptimality arises because an individual eigenmode implements an exponential filter with time constant τ_k = −1/λ_k, and therefore integrates noise before and after stimulus presentation. The ideal observer solution has a narrow integration time window but a long memory time constant, and departs substantially from a simple exponential filter.

F.2 Attractor Networks

For attractor (normal) networks, the eigenvectors form an orthogonal basis with and v_i · v_j = δ_ij. As a consequence, their response SNR follows trivially from the above analysis of individual eigenmodes as we show below. The mean response of an attractor network to the inputs considered above is

and the response covariance is

where the last step assumed that and .

The SNR of network output at time t is

where SNR(v_k, t) is the single-mode SNR given in equation (16). The SNR of an attractor network is therefore maximised when a single eigenmode is aligned to the input discriminant, i.e. v_k = b, in which case the solution reduces to that of single modes considered above, with all other modes having zero SNR.

Thus, attractor networks perform poorly on the working memory task, with SNR decaying exponentially over time due to noise (and decaying as 1/t_d as a function of the optimised delay). We next show how functionally-feedforward (non-normal) and rotational dynamics can improve upon attractor mechanisms.

F.3 Preliminary: The Real and Complex Schur Decomposition

While the eigendecomposition of the matrix A is sufficient to characterise the dynamics and task-performance of a normal network, this description often obscures the transient dynamics of non-normal and rotational networks (Ganguli et al., 2008; Goldman, 2009; Murphy and Miller, 2009). Instead we turned to the Schur decomposition to characterise such networks. Here we briefly describe the real and complex forms of the Schur decomposition that we use extensively in the following sections.

Complex Schur Decomposition

Any matrix A ∈ ℝ^{N ×N} can be expressed in Schur form as A = ULU ^† where U ∈ ℂ^{N ×N} is a unitary matrix, i.e. U ⁻¹ = U ^† where U ^† := (U^⊺)^∗, and L ∈ ℂ^{N ×N} is an upper-triangular matrix, i.e. L_i,i+k = 0 for all k < 0. Given a connectivity matrix A, the off-diagonal elements of L describe “functionally-feedforward” weights between a set of orthogonal neural activity patterns (Schur modes, i.e. columns of U). However, when A has complex eigenvalues, the elements of both U and L become complex-valued, which complicates this interpretation.

Real Schur Decomposition

When some eigenvalues are complex, a more interpretable basis in which to analyse the dynamics is given by the real Schur decomposition A = QT Q^⊺, where Q and T are real-valued matrices, with Q⁻¹ = Q^⊺ (i.e., Q is an orthogonal matrix) and T a block upper-triangular matrix

Here, K is the number of real eigenvalues plus the number of pairs of complex conjugate eigenvalues and the diagonal blocks B_i are either 1×1 submatrices corresponding to real eigenvalues of A (i.e., B_i = λ_k) or 2 × 2 submatrices corresponding to complex conjugate eigenvalue pairs (with for some a, b ∈ ℝ such that the eigenvalues of B_i are equal to a pair of complex-conjugate eigenvalues of A, λ = λ_r ± iλ_ω). Since the columns of Q are orthonormal, the real Schur decomposition provides an orthonormal basis partitioned into one-dimensional non-rotational modes and two-dimensional rotational planes, with functionally-feedforward interactions between these one- and two-dimensional subspaces given by the off-diagonal blocks Ω_ij. Thus, the Ω_ij’s are either 1 × 1, 2 × 1, 1 × 2, or 2 × 2 matrices capturing functionally-feedforward interactions between these modes.

Two-Dimensional Case

A matrix with eigendecomposition and λ ∈ ℝ can be expressed in Schur form as , where ω is a functionally-feedforward weight from the input Schur mode to the output Schur mode q₂ := v₂. Thus, Q = [q₂, q₁] and

A matrix A ∈ ℝ² with complex eigenvalues λ_r ± iλ_ω can be written in real Schur form as A = QT Q^⊺ where

with and Q^⊺ = Q⁻¹ ∈ ℝ^2×2 a real orthogonal matrix. In the second equality, we introduced the parameter ϵ = a/λ_ω, which can be understood as the ellipticity (or eccentricity) of the rotation, and is closely related to the non-normality of the system. In this basis, the major and minor axes of the ellipse are aligned to the two basis axes (defined by the columns of Q). Thus, in real Schur form the matrix A is characterised by four parameters: the real and imaginary parts of the eigenvalues, λ_r and λ_ω, which characterise the damping rate and rotational frequency respectively; the ellipticity of the orbit ϵ; and the angle of the major axis of the ellipse in neural state space (determined by the orthogonal matrix Q). The matrix exponential of A is e^At = Qe^Tt Q^⊺, and the matrix exponential of T can be computed directly as

Importantly, these results for two-dimensional systems also apply to any 2 × 2 submatrix B_i in the N-dimensional real Schur decomposition.

Invariances of the Real Schur Decomposition

The Schur decomposition is non-unique: each ordering of eigenvalues leads to a different Q and T. For the real Schur decomposition, this can be factored into an ordering of blocks and an ordering within 2 × 2 blocks. Within a 2 × 2 block, the ordering of eigenvalues makes only a superficial difference to the resulting decomposition. In particular, if k, k + 1 corresponds to a rotational plane (i.e., if T_k:k+1,k:k+1 = B_i for some i), the change of Schur basis Q = [q₁, q₂, …, q_k, q_k+1, …, q_N] → [q₁, q₂, …, q_k+1, q_k, …, q_N] and T_:,k ↔ T_:,k+1, T_k,: ↔ T_k+1,: yields another real Schur decomposition with ϵ → −1/ϵ. A further invariance for both real and complex Schur decompositions is to reverse the sign of a mode, i.e. Q = [q₁, q₂, …, q_k, …, q_N] → [q₁, q₂, …, q_k, …, q_N], T_k,: → −T_k,:, T_:,k → −T_:,k with ϵ → −ϵ. Through this combination of reordering and sign-reversal of real Schur modes, it is possible to select a real Schur decomposition for which rotations are clockwise or anticlockwise (ϵ ≷ 0) and have their major axis aligned to the x- or y-axis (|ϵ| ≷ 1), and have either positive or negative loading onto the input discriminant (q_k · b ≷ 0). For numerically optimised networks, we utilised these invariances to define a unique Schur decomposition in which 1) the decay time constants τ = −1/λ_r are in descending order from input to output mode 2) elliptical orbits have their major axis aligned to the x axis (|ϵ| > 1) 3) rotations in each plane are clockwise (ϵ > 0) 4) the input mode to each plane has positive alignment with the input discriminant.

F.4 Two-Dimensional Feedforward Networks

We first considered the performance of feedforward (non-normal) networks with real eigenvalues on the working memory task. We restrict our analysis to the stationary state limit t₀ → −∞. Using the Schur basis discussed in Section F.3, the dynamics of a two-dimensional non-normal network with real eigenvalues decouples into two equations:

The mean and variance of responses along the first mode are given by the single-mode analysis in Section F.1:

For the second mode these are given by

and

where we omitted the explicit form of the 𝒪(ω) interaction term for brevity. In the non-normal limit ω → ∞, the response SNR along the output Schur mode becomes:

This SNR behaves non-monotonically as a function of time, reaching a peak at a finite time before decaying back to zero. This non-monotonic timecourse is due to the transfer of information from the input to output mode over time, and the fact that the readout is restricted to the output mode - the SNR of the optimal (time-varying) readout decreases monotonically due to a shift in weighting of the readout from input to output mode over time (Section F.6).

For a given choice of eigenvalues, the SNR is maximal at the delay t_d = log(λ₂/λ₁)/(λ₁ − λ₂). How should the eigenvalues then be chosen in order to optimise response SNR at a particular delay t_d? Writing λ₁ = λ, λ₂ = λ − δλ and taking the limit δλ → 0, we have t_d = log(1 − δλ/λ)/δλ ≈ −1/λ and

where the inequality gives the behaviour of a network with λ optimised for the delay t_d, and is the SNR of a normal network optimised for readout time t_d. Compared to the optimal solution of the normal network with stationary state covariance, the response SNR in this non-normal solution is scaled by a factor of 4e⁻¹ ≈ 1.47. We note the limits ω → ∞ and λ₁ ≈ λ₂ taken above are consistent with the observed behaviour of networks during optimisation before transition to rotational dynamics (Supplementary Figure S3E-F).

F.5 Two-Dimensional Rotational Networks

We next considered a two-dimensional rotational network for which the two eigenvalues and eigenvectors of the weight matrix A form a complex conjugate pair λ, λ^∗ and v, v^∗ (where x^∗ = Re(x) −iIm(x) is complex conjugation). Using the results of Section F.3, the mean and stationary state covariance of the network response can be computed in the real Schur basis as

where we defined ũ = Q^⊺u and assumed that .

The covariance matrix can be inverted directly as

The SNR of the rotational network takes on a simple form when the input linear discriminant is aligned to the one of the Schur modes. In particular, setting Δũ₁ = 0, the SNR is

The SNR when ϵ = 1 is identical to that of a normal network with real eigenvalue λ = λ_r (given by ), independent of the rotational frequency λ_ω. When ϵ → ∞, the SNR derived above matches that observed in optimised networks (Supplementary Figure S8B). This case corresponds to a highly elliptical rotation, with the input linear discriminant aligned to the minor axis of the ellipse. Indeed, when optimising numerically (without an energy penalty), we found that ϵ → ∞ with increasing numbers of iterations, albeit very slowly and with only marginally increasing SNR. The third limit, ϵ → 0, corresponds to a highly elliptical rotation with the input discriminant aligned to the major axis of the ellipse. In this case, noise input to the minor axis of the ellipse is strongly amplified, but the input signal is damped, leading to low SNR.

F.6 N-Dimensional Functionally-Feedforward Networks

We next considered delay line architectures, in which multiple feedforward modes are chained together in sequence (Ganguli et al., 2008; Goldman, 2009). Previous studies have shown that such architectures maximise the SNR of discrete-time linear dynamical systems, provided the number of neurons is equal to or greater than the number of time steps in the delay (Ganguli et al., 2008). However, the optimal dynamics for continuous-time networks, or for discrete-time networks with more time steps than neurons, was not determined. While Goldman (2009) considered a continuous-time network and showed that a delay line architecture can generate responses with memory time constant proportional to the number of neurons, they did not consider noise robustness or show that such architectures are optimal. We derive analytical expressions for the signal-to-noise ratio of a continuous-time delay line and show that, although SNR scales with the number of neurons, it is outperformed by the rotational solutions obtained via numerical optimisation (Supplementary Figure S1, Figure 3F).

F.6.1 Mean and Covariance of Network Response

We considered networks with A_ij = λδ_ij + ωδ_i,j−1 (we choose a lower-triangular rather than upper-triangular decomposition for this analysis as it leads to simpler indices in the following derivations). This network is known as a homogeneous delay line Ganguli et al. (2008). In this basis, the network is structurally feedforward with each element i corresponding to a single neuron. However, a rotation of the connectivity matrix A → QAQ^⊺ produces a functionally-feedforward network for which the feedforward structure acts on Schur modes rather than single neurons.

The matrix exponential for the homogeneous delay line can be computed directly from its power series, , where δ_i≥j is 1 for i ≥ j and zero otherwise. If the network is initialised in the state x_i(0) = δ_i1, then in the absence of input the response propagates down the delay line as . Thus, if the stimulus inputs u(s, t) = u(s)δ(t) + n(t) are such that Δu = ∥Δu∥δ_i1, the trial-averaged responses follow . The station-ary state covariance matrix for the delay line driven by Gaussian white noise input .

F.6.2 Response SNR of Individual Schur Modes

In this section, we place a tractable lower bound on the total response SNR by computing the SNR of the ith Schur mode. Writing , we have

In the limit ω → ∞ the first term in the sum dominates, giving

where SNR_normal(t) is the SNR of a normal network with a left eigenvector aligned to the input discriminant and Stirling’s approximation was used for i ≫ 1. This SNR is maximised for a delay t_d when λ = −(i − 1/2)/t_d, with

Thus, although the optimal response SNR decays as 1/t_d, it also grows with the number of neurons as , and an arbitrarily high SNR can be achieved for any delay t_d by increasing the length of the delay line. Note that the optimal time constant decreases as a function of the length of the delay line as as τ = −1/λ ∼ t_d/N.

F.6.4 Response SNR of Optimal Population Readout

The previous section considered the performance of a single-mode readout, which is in general suboptimal. We next derive a general expression for the SNR of the optimal readout of the delay line, and calculate this SNR explicitly for the two-dimensional case.

The SNR of the network response at time t along a readout w is

This can be written as:

where and . The SNR is maximised using the linear discriminant readout w ∝ P ⁻¹m(t), which gives

Thus, the factor m(t) · P ⁻¹m(t) measures the gain in performance introduced by non-normal integration.

To determine the behaviour of this SNR in the large ω limit, we define , which gives

In the limit ω → ∞, this becomes

where and and equality is achieved by setting .

Two-Dimensional Case

For a two-dimensional network, , so that where β = 2λt and we set . The extrema of this SNR for a given delay t_d are found by solving , which gives β³ + 4β² + 4β + 2 = 0. This equation has three real roots at . The first root (β = 1) has , which is equal to the SNR of a normal network optimised for t_d. The second and third roots have . The root has , while the root has .

We numerically plot the solution to Equation (17) for different values of N (with optimal λ) in Figure 3F), finding an approximately linear scaling of response SNR with the number of neurons N.

F.7 Encoding Multidimensional Stimulus Spaces Using Orthogonal Subspaces

For tasks in which multiple stimuli s ∈ {s₀, s₂, …, s_M−1} can be presented, we considered the SNR for each stimulus pair SNR(s_i, s_j, t). Here, we show how networks optimised for discrimination of two stimuli can naturally be composed to form networks that discriminate arbitrarily many stimuli within a P-dimensional space with SNR for any stimulus pair equal to that of an N-dimensional network optimised exclusively for that stimulus pair. The assumption that inputs span a linear subspace is true in particular for the case of cosine-tuned inputs considered in the main text (where P = 2), but we present results for the general case here.

− ∝ ∈

To construct such a network, we start with an N-dimensional network optimised for discrimination of two stimuli s₁, s₂, with A = QT Q^⊺ and u(s₂) u(s₁) b ℝ ^N. We then use this network to construct a network that discriminates a set of input stimuli that span a P-dimensional subspace of a PN-dimensional space, i.e. with . In particular, we set with

and with with and . Then we have:

We can then compute the signal-to-noise ratio for any stimulus pair

where SNR_2−Stim(s_i, s_j) is the SNR of the original two-stimulus network A for discrimination of a pair of inputs with .

Thus, by aligning an orthogonal subspace to each stimulus dimension , networks optimised for discrimination of a single stimulus pair can be naturally composed to form networks that discriminate arbitrarily many stimuli drawn from a P-dimensional space, with the number of neurons scaling only linearly with P.

G Targeted Dimensionality Reduction

Given a set of measured responses x^(k)(t) of N neurons, over K trials with stimuli stim(k) = θ^(k), and T time points per trial, we fit the model:

where w_p(t) are neuron- and time-dependent factors and f_p(θ^(k)) are stimulus-dependent factors. This equation can be written in matrix form as:

where

Thus, given an estimate of the set of P vectors , a least squares estimate of the factors w_p(t) can be obtained by solving the above linear equation via the Moore-Penrose pseudoinverse W = XF⁺.

To solve for f_p(θ) given estimates of w_p(t), note that:

where contains the PSTHs and is a vector of ones (with 1 ≤ θ ≤ N_stim the stimulus labels). Thus, we obtain . To fit the TDR model to data, we iterated between these two updates until convergence, initialising with f_p(θ) ∼ N(0, 1) (for 1 ≤ p ≤ P and 1 ≤ θ ≤ N_stim).

The TDR model has a number of invariances that we exploited to convert the model into a standard form. First, the stimulus factors can be written as where is the mean over stimuli. We removed the mean component for all p by setting .Second, there exists an invariance for arbitrary invertible M ∈ ℝ^{P ×P}. We selected a unique basis by taking the singular value decomposition and redefining and W_:,1:P → W_:,1:P US (note that this was called the Principal Component Basis by Aoi et al. (2020)). An advantage of this choice of basis is that the vectors f_p are orthonormal, and the stimulus-independent component f₀ = 1 is also orthogonal (since all other vectors have zero mean, and the dot product of a constant vector with a zero mean vector is zero). Finally, for the circular stimulus set in the working memory task, F_1:P,: had two leading principal components of roughly equal magnitude (the two basis vectors that spanned the circle), so that the Principal Component Basis was approximately invariant to a further rotation in this two-dimensional subspace. To ensure that all models (with P = 2) were plotted and compared in a comparable basis, we therefore performed an Orthogonal Procrustes transformation to obtain a 2 × 2 orthogonal matrix and set and W_1:2,: → W_1:2,:R^⊺.

Supplementary Figures

Figure S1:

Figure S2:

Figure S3:

Figure S4:

Figure S5:

Figure S6:

Figure S7:

Figure S8:

Figure S9:

Figure S10:

Acknowledgements

We thank Alex Cayco-Gajic, Arthur Pellegrino, Sina Tootoonian, Matt Nolan and Matt Panichello for valuable feedback. We are particularly grateful to Matt Panichello and Tirin Moore for provision of experimental data. This work was supported in part by the Biotechnology and Biological Sciences Research Council [BB/Y513957/1].

Additional information

Funding

UKRI | Biotechnology and Biological Sciences Research Council (AFRC) (BB/Y513957/1)

• Angus Chadwick