Elucidating the selection mechanisms in context-dependent computation through low-rank neural network modeling

The task paradigm we focused on is the pulse-based CDM task (Pagan et al., 2022), a novel rat-version CDM paradigm inspired by the previous monkey CDM work (Mante et al., 2013). In this paradigm, rats were presented with sequences of randomly-timed auditory pulses that varied in both location and frequency (Figure 1A). In alternating blocks of trials, rats were cued by an external context signal to determine the prevalent location (in the ‘LOC’ context) or frequency (in the ‘FRQ’ context). Note that compared to the continuous sensory input setting in previous works (e.g. Mante et al., 2013), this pulse-based sensory input setting allowed the experimenters to better characterize both behavioral and neural responses (Pagan et al., 2022). We will also demonstrate the unique advantage of this pulse-based input setting later in the present study (e.g. Figure 7).

Figure 1

Download asset Open asset

To solve this task, rats had to select the relevant information for the downstream evidence accumulation process based upon the context. There were at least two different mechanisms capable of performing this selection operation, i.e., selection vector modulation and input modulation (Mante et al., 2013; Pagan et al., 2022). To better introduce these mechanisms, we first reviewed the classical linearized dynamical systems analysis and the concept of selection vector (Figure 1B; Mante et al., 2013; Sussillo and Barak, 2013; Maheswaranathan et al., 2019; Maheswaranathan and Sussillo, 2020; Nair et al., 2023). In the linearized dynamical systems analysis, the neural dynamics around the choice axis (Figure 1B, red line) is approximated by a line attractor model. Specifically, the dynamics in the absence of external input can be approximated by the following linear equation $\frac{dr}{dt}=Mr$, where M is a matrix with one eigenvalue being equal to 0 and all other eigenvalues having a negative real part. For brevity, let us denote the left eigenvector of the 0 eigenvalue as $s$. In this linear dynamical system, the effect of any given perturbation can be decomposed along the directions of different eigenvectors: the projection onto the $s$ direction will remain constant while the projections onto all other eigenvectors will exponentially decay to zero. Thus, for any given input $I$, only the component projecting onto the $s$ direction (i.e. $I\cdot s$) can be integrated along the line attractor (see Methods for more details). In other words, $s$ serves as a vector selecting the input information, which is known as the ‘selection vector’ in the literature (Mante et al., 2013).

Two distinct selection mechanisms can then be introduced based upon the concept of selection vector (Pagan et al., 2022). Specifically, to perform the CDM task, the stimulus input (LOC input, for example) must have a larger impact on evidence accumulation in the relevant context (LOC context) than in the irrelevant context (FRQ context). That is, ${\displaystyle \mathit{I}\cdot \mathit{s}}$ must be larger in the relevant context than in the irrelevant context. The difference between these two can be decomposed into two components:

where the ${\displaystyle \mathrm{\Delta }}$ symbol denotes difference across two contexts (relevant – irrelevant) and the bar symbol denotes average across two contexts (see Methods for more details). The first component ${\displaystyle \mathbf{\Delta }I\cdot \overline{\mathit{s}}}$ is called input modulation in which the change of input information across different contexts is emphasized (Figure 1B, top). In contrast, the second component ${\displaystyle \overline{\mathit{I}}\cdot \mathrm{\Delta }\mathit{s}}$ is called selection vector modulation in which the change of selection vector across contexts is instead highlighted (Figure 1B, bottom).

While these two selection mechanisms were clearly defined, recent work showed that it is actually challenging to differentiate them through neural dynamics (Pagan et al., 2022). For example, both input modulation and selection vector modulation can lead to similar trial-averaged neural dynamics through targeted dimensionality reduction (Figure 1C; Pagan et al., 2022). Take the input modulation as an example (Figure 1C, top). One noticeable aspect we can observe is that the input information (e.g. location information) is preserved in both relevant (LOC context) and irrelevant contexts (FRQ context), which seems contradictory to the definition of input modulation. What is the mechanism underlying this counterintuitive result? As Pagan et al. pointed out earlier, input modulation is not the input change (${\displaystyle \mathrm{\Delta }\mathit{I}}$) per se. Rather, it means the change of input multiplied by selection vector (i.e. ${\displaystyle \mathrm{\Delta }\mathit{I}\cdot \overline{\mathit{s}}}$). Therefore, for the input modulation, while input information indeed is modulated by context along the selection vector direction, input information can still be preserved across contexts along other directions orthogonal to the selection vector, which explains the counterintuitive result and highlights the challenge of distinguishing input modulation from selection vector modulation in experiments (Pagan et al., 2022).

In this study, we sought to address this challenge using the low-rank RNN modeling approach. In contrast to the ‘black-box’ vanilla RNN approach (e.g. Mante et al., 2013), the low-rank RNN approach features both well-controlled model complexity and mechanistic transparency, potentially providing a fresh view into the mechanisms underlying the intriguing selection process. Specifically, the low-rank RNNs we studied here implemented an input-output task structure similar to the classical RNN modeling work of CDM (Figure 1D; Mante et al., 2013). More concretely, the hidden state ${\displaystyle \mathit{x}}$ of a low-rank RNN with $N$ neurons evolves over time according to

where ${\displaystyle J=\sum _{r=1}^R{\mathit{m}}_r{\mathit{n}}_r^T/N}$ is a low-rank matrix with $R$ output vectors ${\displaystyle {\mathit{m}}_r,r=1,\cdots ,R}$ and R input-selection vectors ${\displaystyle {\mathit{n}}_r,r=1,\cdots ,R}$, $\tau$ is the time constant of single neurons, $\varphi$ is the nonlinear activation function, $u_s\left(t\right),s=1,2$ embedded into the network through ${\displaystyle {\mathit{I}}_s}$ mimic the location and frequency click inputs, and $u_s^{ctx}\left(t\right),s=1,2$ embedded through ${\displaystyle {\mathit{I}}_s^{ctx}}$ indicate whether the current context is location or frequency. The output of the network is a linear projection of neural activity (see Methods for more model details). Under this general architectural setting, on one hand, through controlling the rank $R$ of the matrix $J$ during backpropagation training, we can determine the minimal rank required for performing the CDM task and reverse-engineer the underlying mechanism, which will be demonstrated in Figure 2. On the other hand, recent theoretical progress of low-rank RNNs (Beiran et al., 2021; Beiran et al., 2021; Dubreuil et al., 2022) enabled us to explicitly construct neural network models with mechanistic transparency, complementing the reverse-engineering analysis (Mante et al., 2013; Sussillo and Barak, 2013), which will be shown in Figure 3.

Figure 2 with 1 supplement see all

Download asset Open asset

Figure 3 with 1 supplement see all

Download asset Open asset

In the literature, it was found that rank-one RNN models suffice to solve the CDM task (Dubreuil et al., 2022). Here, we further asked whether selection vector modulation can occur in rank-one RNN models. To this end, we trained many rank-1 models (see Methods for details) and found that indeed having rank-1 connectivity (e.g. with the overlap structure listed in Figure 2A; for the detailed connectivity structure, see Figure 2—figure supplement 1) is sufficient to perform the CDM task, consistent with the earlier work. As shown in Figure 2B, in context 1, the decision was made based on input-1 evidence, ignoring input-2 evidence, indicating that the network can effectively filter out irrelevant information. To answer what kind of selection mechanisms underlie this context-dependent computation, we computed the selection vector in two contexts through linearized dynamical systems analysis (Sussillo and Barak, 2013). Cosine angle analysis revealed that selection vectors kept invariant across different contexts (Figure 2C, left), indicating no selection vector modulation. This result was preserved across different hyperparameter settings (such as different regularization coefficients or activation functions). Note that this result actually can be mathematically proved (Pagan et al., 2023; Pagan et al., 2025). Therefore, our modeling result reconfirmed the limitation of rank-one models on selection vector modulation.

While the selection vector was not altered by contexts, the direction of input representations changed significantly across different contexts (Figure 2C, left; see Methods for the definition of input representation direction). Further analysis revealed that the overlap between the input representation direction and the unchanged selection vector is large in the relevant context and small in the irrelevant context, supporting the input modulation mechanism (Figure 2D). These results indicate that while a rank-1 network can perform the task, it can only achieve flexible computation through input modulation (Figure 2E). Importantly, when applying a similar targeted dimensionality reduction method to this rank-1 model, we found that the irrelevant sensory input information was indeed well-represented in neural activity state space (Figure 2F), supporting the conclusion made in the Pagan et al. paper that the presence of irrelevant sensory input in neural state space cannot be used as a reliable indicator for the absence of input modulation (Figure 1C).

In summary, we conclude that to study the mechanism of selection vector modulation, instead of limiting to the simplest model of CDM task, it is necessary to explore network models with higher ranks.

To study the mechanism of selection vector modulation, we designed a rank-3 neural network model, with one additional rank for each sensory input feature (i.e. ${\displaystyle {\mathit{m}}_{i{v}_1}{\mathit{n}}_{i{v}_1}^T}$ for input 1 and ${\displaystyle {\mathit{m}}_{i{v}_2}{\mathit{n}}_{i{v}_2}^T}$ for input 2; Figure 3A, left). Specifically, we ensured that ${\displaystyle {\mathit{I}}_1({\mathit{I}}_2)}$ has a positive overlap with ${\displaystyle {\mathit{n}}_{i{v}_1}({\mathit{n}}_{i{v}_2})}$ and zero overlap with ${\displaystyle {\mathit{n}}_{dv}}$, while ${\displaystyle {\mathit{m}}_{iv1}({\mathit{m}}_{i{v}_2})}$ has a positive overlap with ${\displaystyle {\mathit{n}}_{dv}}$ (Figure 3A, right; see Figure 3—figure supplement 1 and Methods for more details). Moreover, our rank-3 network relies on a multi-population structure, consistent with the notion that higher-rank networks still require a multi-population structure to perform flexible computations (Dubreuil et al., 2022). This configuration implies that the stimulus input 1 (2) is first selected by ${\displaystyle {\mathit{n}}_{i{v}_1}({\mathit{n}}_{i{v}_2})}$, represented by ${\displaystyle {\mathit{m}}_{i{v}_1}({\mathit{m}}_{i{v}_2})}$, and subsequently selected by ${\displaystyle {\mathit{n}}_{dv}}$ before being integrated by the accumulator. In principle, this sequential selection process enables more sophisticated contextual modulations.

We confirmed that a model with such a connectivity structure can perform the task (Figure 3B) and then conducted an analysis similar to that performed for rank-1 models. Unlike the rank-1 model, the selection vector for this rank-3 model changes across contexts while the input representation direction remains invariant (Figure 3C). Further analysis revealed that the overlap between the selection vector and the unchanged input representation direction is large in the relevant context and small in the irrelevant context (Figure 3D), supporting a pure selection vector modulation mechanism (Figure 3E) distinct from the input modulation counterpart shown in Figure 2D. When applying a similar targeted dimensionality reduction method to this rank-3 model, as what we expected, we found that both relevant and irrelevant sensory input information was indeed well-represented in neural activity state space (Figure 3F), which was indistinguishable from the input modulation counterpart (Figure 2F).

Together, through investigating these two extreme cases—one with pure input modulation and the other with pure selection vector modulation, we not only reconfirm the challenge of distinguishing input modulation from selection modulation based on neural activity data (Pagan et al., 2022) but also point out the previously unknown link between selection vector modulation and network connectivity dimensionality.

What is the machinery underlying this link between selection vector modulation and network connectivity dimensionality? One possible way to address this issue is through linearized dynamical systems analysis: first computing the selection vector and the sensory input representation direction through reverse-engineering (Sussillo and Barak, 2013) and then calculating both selection vector modulation and input modulation according to Equation 1. However, the connection between the network connectivity dimensionality and the selection vector obtained through reverse-engineering is implicit and in general non-trivial (Pagan et al., 2023), hindering further investigation of the underlying machinery. Here, by combining recent theoretical progress in low-rank RNNs (Mastrogiuseppe and Ostojic, 2018; Dubreuil et al., 2022) and linearized dynamical systems analysis (Sussillo and Barak, 2013), we introduced a novel pathway-based information flow analysis approach, providing an explicit link between network connectivity, neural dynamics, and selection mechanisms.

To start with, the low-rank RNN dynamics (i.e. Equation 2) can be described by an information flow graph, with each task variable as a node and each effective coupling between task variables as an edge (Mastrogiuseppe and Ostojic, 2018; Dubreuil et al., 2022). Take the rank-1 RNN in Figure 2 as an example. A graph with three nodes, including two input variables ${\kappa }_{p_s}\left(t\right),s=1,2$ and one decision variable ${\kappa }_{dv}\left(t\right)$, suffices (Figure 4A; see Methods for more details). In this graph, the dynamical evolution of the task variable ${\kappa }_{dv}$ can be expressed as:

where the effective coupling ${\displaystyle E_{in{p}_s\to dv}}$ from the input variable ${\displaystyle {\kappa }_{in{p}_s}}$ to the decision variable ${\kappa }_{dv}$ is equal to the overlap between the input representation direction ${\tilde{I}}_s$ (each element is defined by ${\tilde{I}}_{s,i}={\varphi }^{`}\left(x_i\right)I_{s,i}$) and the input-selection vector $n_{dv}$. More precisely, ${\displaystyle E_{in{p}_s\to dv}=⟨{\stackrel{\sim }{I}}_s,n_{dv}⟩}$, where $⟨a,b⟩$ is defined as $\frac{1}{N}{\sum }_1^N{a}_i{b}_i$ for two length-$N$ vectors. Since the input representation direction ${\tilde{I}}_s$ depends on the single neuron gain ${\varphi }^{`}$ and the context input can modulate this gain, the effective coupling ${\displaystyle E_{in{p}_s\to dv}}$ is context-dependent. Indeed, as shown in Figure 4A, Figure 4—figure supplement 1A, ${\displaystyle E_{in{p}_1\to dv}}$ exhibited a large value in context 1 but was negligible in context 2, while ${\displaystyle E_{in{p}_2\to dv}}$ exhibited a large value in context 2 but was negligible in context 1. In other words, information from the input variable can arrive at the decision variable only in the relevant context, which exactly is the computation required by the CDM task.

Figure 4 with 2 supplements see all

Download asset Open asset

To get a more intuitive understanding of the underlying information flow process, we discretized the equation and followed the information flow step by step. Specifically, by discretizing Equation 3 using Euler’s method with a time step equal to the time constant $\tau$ of the system, we get

Take context 1 as an example (Figure 4B, top panel). Initially, there is no information in the system. In step 1, pulse inputs of size $A_1$ and $A_2$ are placed in the $p_1$ and $p_2$ slots, respectively. The information from these slots, after being multiplied by the corresponding effective coupling, then flows to the $dv$ slot. In context 1, ${\displaystyle E_{in{p}_2\to dv}\approx 0}$, meaning only the content from the $p_1$ slot can arrive at the $dv$ slot. Consequently, in step 2, the information content in the $dv$ slot would be ${\displaystyle E_{in{p}_1\to dv}{A}_1}$. The following steps will replicate step 2 due to the recurrent connectivity of the $dv$ slot. The scenario in context 2 is similar to context 1, except that only the content from the $p_2$ slot arrives at the $dv$ slot (Figure 4B, bottom panel). The continuous dynamics for each task variable given pulse input is displayed in Figure 4—figure supplement 2A–C.

The same pathway-based information flow analysis can also be applied to the rank-3 model (Figure 4C and D). In this model, similar to the rank-1 models, there are input variables (${\displaystyle k_{in{p}_1}}$ and ${\displaystyle k_{in{p}_2}}$) and decision variable ${\displaystyle (k_{dv})}$. Additionally, it includes intermediate variables ($k_{i{v}_1}$ and $k_{i{v}_2}$) corresponding to the activity along the $m_{i{v}_1}$ and $m_{i{v}_2}$ axes. In this scenario, instead of flowing directly from the input to the decision variable, the information flows first to the intermediate variables and then to the decision variable. These intermediate variables act as intermediate nodes. Introducing these nodes does more than simply increase the steps from the input variable to the decision variable. In the rank-1 case, the context signals can only modulate the pathway from the input to the decision variable. However, in the rank-3 case, context signals can modulate the system in two ways: from the input to the intermediate variables and from the intermediate variables to the decision variable.

Take the rank-3 model introduced in Figure 3 as an example. Context signals did not alter the representation of input signals, leading to constant effective couplings (i.e. constant ${\displaystyle E_{in{p}_1\to i{v}_1}}$ and ${\displaystyle E_{in{p}_2\to i{v}_2}}$) from input to intermediate variables across contexts. Instead, it changed the effective coupling from the intermediate variables to the decision variable (i.e. large $E_{i{v}_1\to dv}$ in context 1 and near zero $E_{i{v}_1\to dv}$ in context 2; Figure 4C, Figure 4—figure supplement 1B). Consider the context 1 scenario in the discrete case. In step 1, pulse inputs of size $A_1$ and $A_2$ are placed in the $p_1$ and $p_2$ slots, respectively. In step 2, information flows to the intermediate slots, with ${\displaystyle E_{in{p}_1\to i{v}_1}{A}_1}$ in the ${iv}_1$ slot and ${\displaystyle E_{in{p}_2\to i{v}_2}{A}_2}$ in the ${iv}_2$ slot. In step 3, only the information in the ${iv}_1$ slot flows to the $dv$ slot, with the information content being ${\displaystyle E_{in{p}_1\to i{v}_1}{E}_{i{v}_1\to dv}{A}_1}$ (Figure 4D, top panel). The scenario in context 2 is similar to context 1, except that only the content from the ${iv}_2$ slot reaches the $dv$ slot in the third step, with the content being ${\displaystyle E_{in{p}_2\to i{v}_2}{E}_{i{v}_2\to dv}{A}_2}$ (Figure 4D, bottom panel). The continuous dynamics for each task variable given pulse input is displayed in Figure 4—figure supplement 2D, E.

Together, this pathway-based information flow analysis provides an in-depth understanding of how the input information can be routed to the accumulator depending on the context, laying the foundation for a novel pathway-based information flow definition of selection vector and input contextual modulations.

Based on the understanding gained from the pathway-based information flow analysis, we now provide a novel definition of input modulation and selection vector modulation distinct from the one in Equation 1. To begin with, we first considered a model with mixed input and selection vector modulation (Figure 5, left), instead of studying extreme cases (i.e. one with pure input modulation in Figure 2 and the other with pure selection vector modulation in Figure 3). In this more general model, input information can either go directly to the decision variable (with effective coupling $E_{inp\to dv}$) or first pass through the intermediate variables before reaching the decision variable (with effective coupling $E_{inp\to iv}$ and $E_{iv\to dv}$ respectively). Applying the same information flow analysis, we see that a pulse input of unit size will ultimately reach the $dv$ slot with a magnitude of $E_{inp\to dv}+E_{inp\to iv}{E}_{iv\to dv}$ (Figure 5A, right; see Methods for more details). In other words, the total effective coupling $E_{tol}$ from the input to the decision variable is equal to $E_{inp\to dv}+E_{inp\to iv}{E}_{iv\to dv}$. Now, it is straightforward to decompose the context-dependent modulation of $E_{tol}$ in terms of input or selection vector change:

in which the first component $\Delta E_{inp\to dv}+{\Delta E}_{inp\to iv}{\overline{E}}_{iv\to dv}$ stands for the change of input representation (termed as input modulation) and the second component ${\overline{E}}_{inp\to iv}{\Delta E}_{iv\to dv}$ is the one without changing the stimulus input representation (termed as selection vector modulation).

Figure 5

Download asset Open asset

We then asked if this pathway-based definition is equivalent to the one in Equation 1 based on linearized dynamical system analysis (Figure 5; Pagan et al., 2022). To answer this question, we numerically compared these two definitions using a family of models with both input and selection vector modulations (Figure 5C; see Methods for model details) and found that these two definitions produced the same proportion of selection vector modulation across a wide range of parameter regimes (Figure 5D). Together with theoretical derivation of equivalence (see Methods), this consistency confirmed the validity of our pathway-based definition of contextual modulation decomposition.

Having elucidated the pathway-based definitions of input and selection vector modulation, we next provide a novel pathway-based definition of selection vector. For the network depicted in Figure 5A, the total effective coupling $E_{inp\to dv}+E_{inp\to iv}{E}_{iv\to dv}$ can be rewritten as $⟨\tilde{I},n_{tol}⟩$, where $\tilde{I}$ is the input representation direction and $n_{tol}=n_{dv}+⟨{\tilde{m}}_{iv},n_{dv}⟩n_{iv}$. This reformulation aligns with the insight that the amount of input information that can be integrated by the accumulator is determined by the dot product between the input representation direction $\tilde{I}$ and the selection vector. Thus, $n_{tol}$ is the selection vector of the circuit in Figure 5A.

To better understand why selection vector has such a formula, we visualized the information propagation from input to the choice axis using a low-rank matrix (Figure 6A). Specifically, it comprises two components, each corresponding to a distinct pathway. For the first component, input information is directly selected by $n_{dv}$. For the second component, input information is sent first to the intermediate variable and then to the decision variable. This pathway involves two steps: in the first step, the input representation vector is selected by $n_{iv}$, and in the second step, to arrive at the choice axis, the selected information has to be multiplied by the effective coupling $E_{iv\to dv}=⟨{\tilde{m}}_{iv},n_{dv}⟩$. By concatenating these two steps, information propagation from the input to the choice axis in this pathway can be effectively viewed as a selection process mediated by the vector $⟨{\tilde{m}}_{iv},n_{dv}⟩n_{iv}$ (termed as the second-order selection vector component). Therefore, $n_{tol}$ provides a novel pathway-based definition of selection vector in the network. We further verified the equivalence between this pathway-based definition and the linearized-dynamical-systems-based classical definition (Mante et al., 2013) in our simple circuit through theoretical derivation (see Methods) and numerical comparison (Figure 6B).

Figure 6

Download asset Open asset

To visualize the pathway-based selection vector in neural activity state space, we found that a minimum of three dimensions is required, including the input representation direction, the decision variable representation direction, and the intermediate variable representation direction (Figure 6C, left). This geometric visualization highlighted the role of extra dimensions beyond the classical two-dimensional neural activity space spanned by the line attractor and selection vector (Figure 6C, right) in accounting for the selection vector modulation. This is simply because only the second-order selection vector component, which depends on the existence of the intermediate variable, is subject to contextual modulation. In other words, without extra dimensions to support intermediate variable encoding, there will be no selection modulation.

Together, this set of analyses provided a parsimonious pathway-based understanding for both selection vector and its contextual modulation.

The new insights obtained from our new framework enable us to generate testable predictions for better differentiating selection vector modulation from input modulation, a major challenge unresolved in the field (Pagan et al., 2022).

First, we predict that it is more likely to have a large proportion of selection vector modulation for a neural network with high-dimensional connectivity. To better explain the underlying rationale, we can simply compare the number of potential connections contributing to the input modulation with those contributing to the selection vector modulation for a given neural network model. For example, in the network presented in Figure 5, there are three connections (including light blue, dark blue, and green ones) while only one connection (i.e. the green one) supporting selection vector modulation. For a circuit with many higher-order pathways (e.g. Figure 7A), only those connections with the input as the sender is able to support input modulation. In other words, there exist far more connections potentially eligible to support the selection vector modulation (Figure 7B), thereby leading to a large selection vector modulation proportion. We then tested this prediction on vanilla RNNs trained through backpropagation (Figure 7C; Mante et al., 2013; Song et al., 2016; Yang and Wang, 2020). Using effective dimension (see Methods for a formal definition; Rudelson and Vershynin, 2007; Sanyal et al., 2020) to quantify the dimensionality of the connectivity matrix, we found a strong positive correlation between the effective dimension of connectivity matrix and the selection vector modulation (Figure 7D, left panel and Figure 7—figure supplement 1A, see Method for details).

Figure 7 with 4 supplements see all

Download asset Open asset

We then asked if we could generate predictions to quantify the proportion of selection vector modulation purely based on neural activities. To this end, as what has been performed in Pagan et al., we took advantage of the pulse-based sensory input setting to calculate the single neuron response kernel (see Methods for details). For a given neuron, the associated response kernel is defined to characterize the influence of a pulse input on its firing rate in later times. For example, for a neuron encoding the decision variable, the response kernel should exhibit a profile accumulating evidence over time (Figure 7E, top left). In contrast, for a neuron encoding the sensory input, the response kernel should exhibit an exponential decay profile with time constant $\tau$ (Figure 7E, top right). For each RNN trained through backpropagation, we then examined the associated single neuron response kernels. We found that for the model with low effective dimension (denoted by a star marker in Figure 7D), there were mainly two types of response kernels, including sensory input profile and decision variable profile (Figure 7E, top). In contrast, for the model with the highest effective dimension (denoted by a square marker in Figure 7D), aside from the sensory input and decision variable profiles, richer response kernel profiles were exhibited (Figure 7E, bottom). In particular, there was a set of single neuron response kernels with peak amplitudes occurring between the pulse input onset and the choice onset (Figure 7E, bottom right). These response kernels cannot be explained by the combination of sensory input and decision variables. Instead, the existence of these response kernels signifies neural dynamics in extra dimensions beyond the subspaces spanned by the input and decision variable, the genuine neural dynamical signature of the existence of selection vector modulation (Figure 6C).

While single neuron response kernels are illustrative in highlighting the model difference, they lack explanatory power at the population level. Therefore, we employed the singular value decomposition method to extract the principal dynamical modes of response kernels at the population level (see Methods for details). We found that similar dynamical modes, including one persistent choice mode (gray) and three transient modes (blue, orange, green), were shared across both the low and high effective dimension models (Figure 7F, left). The key difference between these two models lies in the percentage of explained variance (PEV) of the second transient mode (orange): while there is near-zero PEV in the low effective dimension model (Figure 7F, top right), there is substantial PEV in the high effective dimension model (Figure 7F, bottom right), consistent with the single neuron picture shown in Figure 7E. This result led us to use the PEV of extra dynamical modes (including the orange and green ones; see Methods for details) as a simple index to quantify the amount of selection vector modulation in these models. As expected, we found that the PEV of extra dynamical modes can serve as a reliable index reflecting the proportion of selection vector modulation in these models (Figure 7G, right panel and Figure 7—figure supplement 1B and C). Similar results for vanilla RNNs trained with different hyperparameter settings are displayed in Figure 7—figure supplement 2.

Together, we identified novel neural dynamical signatures of section vector modulation at both the single neuron and population level, suggesting the potential great utility of these neural dynamical signatures in distinguishing the contribution of selection vector modulation from input modulation in experimental data.

In their seminal work, Mante, Sussillo, and their collaborators developed a numerical approach to compute the selection vector for trained RNN models (Mante et al., 2013). Based on this concept of selection vector, recently, Pagan et al. proposed a new theoretical framework to decompose the solution space of context-dependent decision-making, in which input modulation and selection vector modulation were explicitly defined (i.e. Equation 1). Here, taking the theoretical advantage of low-rank RNNs (Mastrogiuseppe and Ostojic, 2018; Dubreuil et al., 2022), we went beyond numerical reverse-engineering and provided a complementary pathway-based definition of both selection vector and selection vector modulation (i.e. Equation 5). This new definition gained us a novel geometric understanding of selection vector modulation, revealed a previously unknown link between extra dimensions and selection vector modulation (Figure 6C), and eventually provided us with experimentally identifiable neural dynamical signature of selection vector modulation at both the single neuron and population levels (Figure 7, E-G).

One hallmark of higher cognition is individual variability, as the same higher cognition problem can be solved equally well with different strategies. Therefore, studying the neural computations underlying individual variability is no doubt of great importance (Hariri, 2009; Parasuraman and Jiang, 2012; Keung et al., 2020; Nelli et al., 2023). Recent experimental advances enabled researchers to investigate this important issue in a systematic manner using delicate behavioral paradigms and large-scale recordings (Pagan et al., 2022). However, the computation underlying higher cognition is largely internal, requiring discovering novel neural activity patterns as internal indicators to differentiate distinct circuit mechanisms. In the example of context-dependent decision-making studied here, to differentiate selection vector modulation from input modulation, we found the PEV in extra dynamical modes is a reliable index for a wide variety of RNNs (Figure 7D, Figure 7—figure supplement 2). However, cautions have to be made here as we can conveniently construct counter-examples deviating from the picture depicted by this index. For instance, manually introducing additional dimensions that do not directly contribute to the computation can disrupt the index (Figure 7—figure supplement 4A and B). In the extreme scenario, we can construct two models with distinct circuit mechanisms (selection vector modulation and input modulation, respectively) but having the same neural activities (Figure 7—figure supplement 3), suggesting that any activity-based index alone would fail to make this differentiation. Then, why did the proposed index work for the trained vanilla RNNs shown in Figure 7D? Our lesion analysis suggests that the underlying reason is that the major variance in neural activity of vanilla RNNs learned through backpropagation is task-relevant (Figure 7—figure supplement 4C, see Methods for details). However, it is highly likely that task-irrelevant neural activity variance exists in higher brain regions, meaning the proposed index may not perform well in neural recordings. Therefore, our modeling work suggests that, to address the intricacy of individual variability of neural computations underlying higher cognition, integrative efforts incorporating not only large-scale neural activity recordings but also activity perturbations, neuronal connectivity knowledge, and computational modeling may be inevitably required.

While we mainly focused on context-dependent decision-making tasks in this study, the issue of whether input or selection vector modulation prevails is not limited to the domain of decision-making. For instance, recent work by Chen et al., 2024 demonstrated that during sequence working memory control (Botvinick and Watanabe, 2007; Xie et al., 2022), sensory inputs presented at different ordinal ranks first entered into a common sensory subspace and then were routed to the corresponding rank-specific working memory subspaces in monkey frontal cortex. Here, similar to the decision-making case (Figure 1C), the same issue arises: where is the input information selected by the context (here, the ordinal rank)? Can the presence of a common sensory subspace (similar to the presence of location information in both relevant and irrelevant contexts in Figure 1C) preclude the input modulation? The pathway-based understanding of input and selection vector modulation gained from CDM in this study may be transferable to address these similar issues.

In this study, we linked the selection vector modulation with transient dynamics (Aoi et al., 2020; Soldado-Magraner et al., 2023) in extra dimensions. While the transient dynamics in extra dimensions are not necessary in context-dependent decision-making here (Dubreuil et al., 2022; Figure 2), more complex context-dependent computation may require its presence. For example, recent work by Tian et al., 2024 found that transient dynamics in extra subspaces is required to perform the switch operation (i.e. exchanging information in subspace 1 with information in subspace 2). Understanding how transient dynamics in extra dimensions contribute to complex context-dependent computation warrants further systematic investigation.

In summary, through low-rank neural network modeling, our work provided a parsimonious mechanistic account for how information can be selected along different pathways, making significant contributions towards understanding the intriguing selection mechanisms in context-dependent computation.

We investigated networks of $N$ neurons with $S$ input channels, described by the following temporal evolution equation

In this equation, $x_i\left(t\right)$ represents the activation of neuron $i$ at time $t$, $\tau$ denotes the characteristic time constant of a single neuron, and $\varphi$ is a nonlinear activation function. Unless otherwise specified, we use the tanh function as the activation function. The coefficient $J_{ij}$ represents the connectivity weight from neuron $j$ to neuron $i$. The input $u_s\left(t\right)$ corresponds to the $s$-th input channel at time $t$, with feedforward weight $I_{si}$ to neuron $i$, and ${ϵ}_i\left(t\right)$ represents white noise at time $t$. The network’s output is obtained from the neuron’s activity $\varphi \left(x\right)$ through a linear projection:

The connectivity matrix $J$, specified as $J=\left\{J_{ij}\right\},i=1,\dots ,N,j=1,\dots ,N$, can be a low-rank or a full-rank matrix. In the low-rank case, $J$ is restricted to a low-rank matrix $\frac{1}{N}{\sum }_{r=1}^R{m}_r{n}_r^T$, in which $m_r$ is the $r$-th output vector, and $n_r$ is the $r$-th input-selection vector, with each element of $m_r$ and $n_r$ considered an independent parameter (Dubreuil et al., 2022). In this paper, we set the time constant $\tau$ to be 100 ms, and use Euler’s method to discretize the evolution equation with a time step $\Delta t=20ms$.

We modeled the click-version CDM task recently investigated by Pagan et al., 2022. The task involves four input channels $u_1\left(t\right)$, $u_2\left(t\right)$, $u_1^{ctx}\left(t\right)$, and $u_2^{ctx}\left(t\right)$, where $u_1\left(t\right)$ and $u_2\left(t\right)$ are stimulus inputs and $u_1^{ctx}\left(t\right)$ and $u_2^{ctx}\left(t\right)$ are context inputs. Initially, there is a fixation period lasting for $T_{fix}=200ms$. This is followed by a stimulus period of $T_{sti}=800ms$ and then a decision period of $T_{decision}=20ms$.

For trial $k$, at each time step, the total number of pulse inputs (#pulse) is sampled from a Poisson distribution with a mean value of $40\Delta t=0.8$. Each pulse has two properties: location and frequency. Location can be either left or right, and frequency can be either high or low. We randomly sample a pulse to be right with probability ${\displaystyle p^k}$ (hence left with probability ${\displaystyle 1-p^k}$) and to be high with probability $p_{high}^k$ (hence low with probability $1-p_{high}^k$). The values of ${\displaystyle p^k}$ and $p_{high}^k$ are independently chosen from the set {39/40, 35/40, 25/40, 15/40, 5/40, 1/40}. The stimulus strength for the location input at trial $k$ is defined as ${\displaystyle 2\times p^k-1}$, and for the frequency input, it is defined as $2\times p_{high}^k-1$. The input $u_1\left(t\right)$ represents location evidence, calculated as 0.1 times the difference between the number of right pulses and the number of left pulses (#right-#left) at time step $t$. The input $u_2\left(t\right)$ represents frequency evidence, calculated as 0.1 times the difference between the number of high-frequency pulses and the number of low-frequency pulses (#high-#low) at time step $t$.

The context is randomly chosen to be either the location context or the frequency context. In the location context, $u_1^{ctx}=1$ and $u_2^{ctx}=0$ throughout the entire period. The target output $z_k$ is defined as the sign of the location stimulus strength (1 if ${\displaystyle p^k>0.5}$, otherwise –1). Thus, in this context, the target output $z_k$ is independent of the frequency stimulus input. Conversely, in the frequency context, $u_1^{ctx}=0$ and $u_2^{ctx}=1$. The target output $z_k$ is defined as the sign of the frequency stimulus strength (1 if ${\displaystyle p_{high}^k>0.5}$, otherwise –1). Thus, in this context, the target output is independent of the location stimulus input.

To uncover the computational mechanism enabling each RNN to perform context-dependent evidence accumulation, we utilize linearized dynamical system analysis to ‘open the black box’ (Mante et al., 2013). The dynamical evolution of an RNN in context $c$ is given by:

where $r=\varphi \left(x\right)$ is the neuron activity. First, we identify the slow point of each RNN in each context using an optimization method (Mante et al., 2013). Let ${\displaystyle x}$ be the discovered slow point, i.e., ${\displaystyle x\approx J{r}^{\ast }+I_c^{ctx}}$ where ${\displaystyle r=\varphi (x^{\ast })}$. We define a diagonal matrix ${\displaystyle G=\text{diag}\left({\varphi }^{\prime }\left(x^{\ast }\right)\right)}$, with the $i$-th diagonal element representing the sensitivity of the $i$-th neuron at the slow point. Near the slow point, we have ${\displaystyle \mathrm{\Delta }r=r-r^{\ast }\approx G{\varphi }^{\prime }\left(\mathrm{\Delta }x\right)}$. Then, we can derive:

Thus, the dynamics of neuron activity around the slow point can be approximated by a linear dynamical system:

where $E$ is the identity matrix, $M$ denotes the state transition matrix of the dynamical system, and ${\tilde{I}}_r=G{I}_r,r=1,2$ are stimulus input representation directions. Similar to previous work (Mante et al., 2013), we find that for every network, the linear dynamical system near the slow point in each context is roughly a linear attractor. Specifically, the transition matrix $M$ has a single eigenvalue close to zero, while all other eigenvalues have negative real parts. The right eigenvector of $M$ associated with the eigenvalue close to zero defines the stable direction of the dynamical system, forming the line attractor $\rho$(unit norm). The left eigenvector of $M$ associated with that eigenvalue defines the direction of the selection vector $s$. The norm of selection vector $s$ is chosen such that $s\cdot \rho =1$. Previous work (Mante et al., 2013) has shown that a perturbation $\Delta r_0$ from the line attractor will eventually converge to the line attractor, with the distance from the starting point being $s\cdot {\Delta r}_0$.

Based on linearized dynamical systems analysis, Pagan et al. recently defined input modulation and selection vector modulation (Pagan et al., 2022) as:

Specifically, the input modulation and selection vector modulation for stimulus input 1 are defined as ${\displaystyle (s^{ct{x}_1}+s^{ct{x}_2})/2\cdot \left({\stackrel{\sim }{I}}_1^{ct{x}_1}-{\stackrel{\sim }{I}}_1^{ct{x}_2}\right)}$ and $(s^{ct{x}_1}-s^{ct{x}_2})\cdot ({\tilde{I}}_1^{ct{x}_1}+{\tilde{I}}_1^{ct{x}_2})/2$, respectively. Similarly, the input modulation and selection vector modulation for stimulus input 2 are defined as ${\displaystyle (s^{ct{x}_1}+s^{ct{x}_2})/2\cdot \left({\stackrel{\sim }{I}}_2^{ct{x}_2}-{\stackrel{\sim }{I}}_2^{ct{x}_1}\right)}$ and $(s^{ct{x}_2}-s^{ct{x}_1})\cdot ({\tilde{I}}_2^{ct{x}_1}+{\tilde{I}}_2^{ct{x}_2})/2$, respectively. Proportion for selection vector modulation (Figure 7D and G) is defined as $\frac{mo{d}_{sel}}{mo{d}_{inp}+mo{d}_{sel}}$.

The rank-1 RNNs in Figure 2 are trained using backpropagation-through-time with the PyTorch framework. These networks are trained to minimize a loss function defined as:

Here, $z_k$ is the target output, ${\widehat{z}}_{k,t}$ is the network output, and the indices $t$ and $k$ represent time and trial, respectively. $M_t$ is a temporal mask with value {0, 1}, where $M_t$ is 1 only during the decision period. $L_{reg}$ is the L2 regularization loss. For full-rank RNNs (Figure 7), $L_{reg}=w_{reg}{\sum }_{ij}{J}_{ij}^2$ and for low-rank RNNs, $L_{reg}=w_{reg}{\sum }_r({‖m_r‖}_2^2+{‖n_r‖}_2^2)$. The loss function is minimized by computing gradients with respect to all trainable parameters. We use the Adam optimizer in PyTorch with the decay rates for the first and second moments set to 0.9 and 0.99, respectively, and a learning rate of ${10}^{-3}$. Each RNN is trained for 5000 steps with a batch size of 256 trials.

For Figure 2, we trained 100 RNNs of $N=512$ neurons with the same training hyperparameters. The rank of the connectivity matrix is constrained to be rank-1, represented as $J=\frac{1}{N}{m}_{dv}{n}_{dv}^T$. We trained the elements of the input vectors $I_1,I_2,I_1^{ctx},I_2^{ctx}$, connectivity vector $m_{dv},n_{dv}$, and the readout vector $w$. All trainable parameters were initialized with random independent Gaussian weights with a mean of 0 and a variance of $1/N^2$. The regularization coefficient $w_{reg}$ is set to ${10}^{-4}$.

For the trial-averaged analyses shown in Figures 2f and 3f, we followed a procedure similar to Mante et al., 2013. Specifically, at each time point $t$ and for each neuron $i$, we fit the following linear regression model to characterize how different task variables contribute to that neuron’s activity:

Here, $choice\left(k\right)$, $u_{1,k}$, $u_{2,k}$, and $context\left(k\right)$ represent the values of the corresponding variables on trial k. Next, for each regressor (choice, input 1, input 2, context), we pooled the resulting regression coefficients across neurons to get the time-vary regression vectors ${\displaystyle {\mathit{\beta }}_{choice;t},{\mathit{\beta }}_{in{p}_1;t},{\mathit{\beta }}_{in{p}_2;t,}{\mathit{\beta }}_{context;t}}$. For regressor $v$, we identified the time $t_v^{max}$ when the regression vector ${\beta }_{v;t}$ has maximum norm, and got the time-independent regression vectors:

Next, we assembled these vectors into a matrix ${\displaystyle B=\left[{\mathit{\beta }}_{choice}{\mathit{\beta }}_{in{p}_1}{\mathit{\beta }}_{in{p}_2}{\mathit{\beta }}_{context}\right]}$ and used QR decomposition to get the orthogonalized regression basis ${\displaystyle \left[{\mathit{\beta }}_{choice}^{\mathrm{\perp }},{\mathit{\beta }}_{in{p}_1}^{\mathrm{\perp }},{\mathit{\beta }}_{in{p}_2}^{\mathrm{\perp }},{\mathit{\beta }}_{context}^{\mathrm{\perp }}\right]}$. Finally, we averaged neuronal activity across trials sharing the same condition (choice, context, input 1), and then projected this average activity onto the choice and input 1 axes. This process resulted in the trial-averaged population dynamics illustrated in Figures 2F and 3F.

The transition matrix of neuron activity in the rank-1 RNN is given by:

Multiplying $M^T$ on the right by ${\displaystyle {\mathit{n}}_{dv}}$, we obtain

where the $⟨\cdot ,\cdot ⟩$ symbol is defined as ${\displaystyle ⟨\mathit{a},\mathit{b}⟩=\frac{1}{N}\sum _{i=1}^N{a}_i{b}_i}$ for two vectors of length-$N$. The requirement for the linear attractor approximation, ${\displaystyle ⟨{\stackrel{\sim }{\mathit{m}}}_{dv},{\mathit{n}}_{dv}⟩\approx 1}$, is met by all our trained and handcrafted RNNs (Figure 4—figure supplement 1). This demonstrates that $n_{dv}$ is the left eigenvector of the transition matrix, and hence $\frac{1}{N}{‖{\tilde{m}}_{dv}‖}_2{n}_{dv}$ is the selection vector in each context. Therefore, the direction of the selection vector is invariant across different contexts for the rank-1 model, indicating no selection vector modulation and this is consistent with the training results shown in Figure 2.

First, we provide the implementation details of the rank-3 RNNs used in Figure 3. The network consists of 30,000 neurons divided into three populations, each with 10,000 neurons. The first population (neurons 1–10,000) receives stimulus input, accounting for both the information flow from the stimulus input to the intermediate variables and the recurrent connection of the decision variable. The second (neurons 10,001–20,000) and third populations (neurons 20,001–30,000) handle the information flow from the intermediate variables to the decision variable and are modulated by the context signal. To achieve this, we generate three Gaussian random matrices ${\displaystyle (M^{\left(1\right)},M^{\left(2\right)},M^{\left(3\right)})}$ of shape 10,000×3. Let $M_{:,r}^{\left(p\right)}$ denote the $r$-th column of matrix $M^{\left(p\right)}$. The stimulus input $I_1$ is given by the concatenation of three length-10,000 vectors $\left[M_{:,1}^{\left(1\right)};0;0\right]$, 0 where denotes a length-10,000 zero vector. The stimulus input $I_2$ is given by $\left[M_{:,2}^{\left(1\right)};0;0\right]$. The context input $I_1^{ctx}$ is given by $\left[0;M_{:,1}^{\left(2\right)};0\right]$. The context input $I_2^{ctx}$ is given by $\left[0;0;M_{:,1}^{\left(3\right)}\right]$. The connectivity vectors $m_{iv1},m_{iv2}$ and $m_{dv}$ are given by $\left[0;M_{:,2}^{\left(2\right)};M_{:,2}^{\left(3\right)}\right]$, $\left[0;M_{:,3}^{\left(2\right)};M_{:,3}^{\left(3\right)}\right]$ and $\left[M_{:,3}^{\left(1\right)};0;0\right]$, respectively. The input-selection vectors $n_{iv1},n_{iv2},$ and $n_{dv}$ are given by $\left[10M_{:,1}^{\left(1\right)};0;0\right]$, $\left[10M_{:,2}^{\left(1\right)};0;0\right]$, and $\left[3M_{:,3}^{\left(1\right)};-g{M}_{:,2}^{\left(2\right)}+M_{:,3}^{\left(2\right)};M_{:,2}^{\left(3\right)}-g{M}_{:,3}^{\left(3\right)}\right]$, respectively, where $g={\int }_{-\infty }^{+\infty }\frac{1}{\sqrt{2\pi }}{e}^{-x^2/2}{\varphi }^{`}\left(x\right)dx$ represents the average gain of the second population in context 1 or the third population in context 2. The readout vector $w$ is given by $\left[{4M}_{:,3}^{\left(1\right)};0;0\right]$. We generate 100 RNNs based on this method to ensure that the conclusions in Figure 3 do not depend on the specific realization of random matrices. Linearized dynamical system analysis reveals that all these RNNs perform flexible computation through pure selection vector modulation. Please see the section ‘The construction of rank-3 RNN models’ for a mean-field-theory-based understanding.

The dynamical evolution equation for low-rank RNN with $R$ ranks and $S$ input channels in context c is given by

Assuming $x\left(0\right)=I_c^{ctx}$ at $t=0$, the dynamics of $x\left(t\right)-x\left(0\right)$ are always constrained in the subspace spanned by $\left\{m_r,r=1,\dots ,R\right\}$ and {$I_s,s=1,\dots ,S\}$. Therefore, $x\left(t\right)$ can be expressed as a linear combination of these vectors: $x\left(t\right)=I_c^{ctx}+{\sum }_{r=1}^R{k}_r\left(t\right)m_r+{\sum }_{s=1}^S{k}_{p_s}\left(t\right)I_s,$ leading to the following evolving dynamics of task variables:

where $k_{p_s}\left(t\right)$ denotes activation along the $s$-th input vector, termed the input task variable and $k_r\left(t\right)$ denotes activation along the $j$-th output vector, termed the internal task variable.

Therefore, for rank-1 RNNs, the latent dynamics of decision variable (internal task variable associated with $m_{dv}$) in context $c$ is given by

Similar to the linearized dynamical systems analysis introduced earlier, we linearized Equation 15 around $I_c^{ctx}$, obtaining the following linearized equation:

where ${\tilde{m}}_{dv}^{ct{x}_c}=G_c{m}_{dv}$ (termed as the decision variable representation direction) and ${\tilde{I}}_s^{ct{x}_c}=G_c{I}_s$ (termed as the input representation direction), with $G_c$ equal to $\text{diag}\left({\varphi }^{`}\left(I_c^{ctx}\right)\right)$. By denoting $\frac{1}{N}{n}_{dv}^T{\tilde{m}}_{dv}^{ct{x}_c}$ and $\frac{1}{N}{n}_{dv}^T{\tilde{I}}_s^{ct{x}_c}$ as $E_{dv\to dv}^{ct{x}_c}$ and $E_{p_s\to dv}^{ct{x}_c}$, respectively, together with the fact that $\frac{1}{N}{n}_{dv}^T\varphi \left(I_c\right)$ is close to zero for all trained rank-1 RNNs, we obtain

which is Equation 3 in the main text.

Using a similar method, we can uncover the latent dynamics of rank-3 RNNs, as shown in Figure 5. Note that the rank-3 RNN in Figure 3 is a special case of this more general form. The latent dynamics for internal task variables in context $c$ can be written as:

By applying the same first-order Taylor expansion, we obtain the following equations:

We consider the case in which the intermediate variables ($k_{i{v}_s},s=1,2$, internal task variables associated with ${\displaystyle {\mathit{m}}_{i{v}_s},s=1,2}$, respectively) only receive information from the corresponding stimulus input, and the effective coupling of the recurrent connection is 1 in both contexts. Specifically, we assume:

Our construction methods for rank-3 RNNs in Figures 3 and 5 guarantee these conditions when the network is large enough (for example, $N=\mathrm{30,000}$ in our setting). Under these conditions, Equations 20 and 21 can be simplified to:

Suppose at $t=0,$ the network receives a pulse from input 1 with size $A_1$ and a pulse from input 2 with size $A_2$, which correspond to $u_1\left(t\right)=A_1\tau \delta \left(t\right)$, $u_2\left(t\right)=A_2\tau \delta \left(t\right)$. Under this condition, the expression for $k_{dv}$ is given by

From Equation 27, as $t\to \infty$, $k_{dv}$ will converge to ${\displaystyle \sum _{s=1}^2{A}_s\left(E_{in{p}_s\to dv}^{ct{x}_c}+E_{in{p}_s\to i{v}_s}^{ct{x}_c}{E}_{i{v}_s\to dv}^{ct{x}_c}\right)}$, providing a theoretical basis for the pathway-based information flow formula presented in Figures 4 and 5.