Neural circuit models for evidence accumulation through choice-selective sequences

Domain: biology
Published: 2026-03-17
Editorial outcome: minor_revision

Abstract
Decision making is traditionally thought to be mediated by neurons that accumulate evidence through persistent activity. However, recent decision-making experiments in rodents have observed neurons across the brain that fire sequentially, rather than persistently, with the subset of neurons in the sequence depending on the animal’s choice. We developed two candidate circuit models in which neurons are active sequentially and transfer evidence faithfully to the next active population. One model encodes evidence in the relative firing of two competing chains of neurons, and the other in the network location of a stereotyped, bump-like pattern of neural activity. Neural recordings from four brain regions during an evidence accumulation task revealed that different regions displayed evidence tuning consistent with different candidate models. This work provides a mechanistic explanation for how graded information may be precisely accumulated within and transferred between neural populations, and suggests that different brain regions may accumulate evidence through different circuit mechanisms.

Introduction
The accumulation of evidence for different alternatives is thought to be a fundamental cognitive operation. As such, neural recordings have been performed in many tasks designed to probe how evidence is accumulated to make a decision, including tactile1,2, visual3–11, olfactory12, auditory13,14, motor15, and value discrimination tasks16–20. Experiments using these tasks have been conducted with the goal of understanding how the neural circuitry within and across different brain regions contributes to the accumulation of evidence, subsequent maintenance of this information in working memory, and ultimate commitment to a decision based on this evidence21.

A major class of models proposed to describe this decision making process is the drift diffusion model22–24, in which evidence is accumulated with a noisy drift process, and decisions are made based on when the resulting accumulation reaches some threshold. Many circuit-based models have been shown to be equivalent to the drift-diffusion model25. Such circuit-based models include mutual inhibition models26, feedforward inhibition models27, and pooled inhibition models28. These basic circuit models have been extended to explain a number of different decision-making experiments with recordings in different regions where neurons show persistent activity29–34.

The mechanistic essence of the circuit-based models is the accumulation of evidence along a low-dimensional attractor representing the decision variable. Thus, these models predict that the activity of single neurons ramps up or down with evidence throughout a decision-making task35,36. Such ramping activity during decision-making tasks has been observed in a number of different brain regions in different paradigms, including in the lateral intraparietal area (LIP), medial temporal area (MT), and ventral intraparietal area (VIP) of monkeys in a random dot motion discrimination task37–40, in the dorsal premotor cortex of monkeys during discrimination tasks41,42, and in the posterior parietal cortex (PPC) and frontal orienting fields (FOF) of rats in an auditory evidence accumulation task14,43. These classes of neural circuit models that accumulate evidence through persistent, ramping activity successfully capture much observed experimental data.

However, a mounting body of evidence from large scale recordings during a range of decision-making tasks suggests that, in many decision-making contexts, neurons across cortical and subcortical regions are not persistently active. Instead, neurons fire transiently and sequentially in a choice-specific manner. This occurs in decision-making tasks both with navigation44–49 and without navigation50–58, including in evidence accumulation tasks44,49. These data challenge existing models of evidence accumulation by contradicting a core tenet of previous models: the presence of persistent, rather than sequential, neural activity.

For a model to accumulate evidence through sequential neural firing, it must perform two key computations. First, as in traditional models based on persistent activity, it must accumulate evidence received at each position. Second, unlike traditional models, it must transfer information between neurons at sequential positions.

Here, we develop two classes of neural circuit models that perform these computations. In each model class, the transfer of accumulated evidence between neurons encoding sequential positions is accomplished via a position-modulated gating signal. The key difference in the models is in their encoding of evidence, which is either in the relative amplitude of activity in two competing chains of neurons or in the network location of a stereotypically shaped, unimodal pattern (bump) of neuronal firing. We tested these differentiating evidence tuning predictions in a dataset of 14,247 imaged neurons, collected across 4 cortical and subcortical regions in mice performing an evidence accumulation task in virtual reality59. This analysis revealed that while seemingly similar choice-selective sequences are present in this task across these regions, neurons in different regions differ in their evidence tuning properties. In the neocortex, neuronal tuning curves mainly exhibit a monotonically increasing or decreasing encoding of evidence, characteristic of the competing chains model class, whereas in the hippocampus neurons tend to have narrow, non-monotonic tuning curves, consistent with the bump model. Thus, our results reveal different encodings of accumulated evidence across brain regions and suggest that different regions could use distinct circuit mechanisms to form their evidence representations.

Results
Traditional Integrator Models Fail to Explain Choice-Selective Sequences in the Accumulating Towers Task

The models presented below are motivated by a large dataset (n = 14,247 neurons; N=26 mice) of previously unpublished calcium imaging data (see Methods) from anterior cingulate cortex (ACC) and dorsomedial striatum (DMS), as well as recently published imaging data from hippocampus (HPC)49 and retrosplenial cortex (RSC)44. These data were all obtained during the same navigation-based, accumulation of evidence task in virtual reality (accumulating towers task, Fig. 1A). In this task59 (see Methods), mice navigate a T-maze with Poisson counts of visual cues in the form of towers appearing to each side. After a delay region with no cues, the mouse reaches the arms of the maze and is rewarded for turning to the side with more cues (Fig. 1A). To make the correct decision, the mouse must integrate the number of cues and hold this information in working memory during the delay region. Mice are able to perform this task with good accuracy that increases with the magnitude of the difference between the number of cues to each side59.

Across brain regions, neurons fired sequentially rather than persistently (Fig. 1B-E; Supplementary Fig. 1). Such sequences were present throughout the task, beginning in the cue region and continuing throughout the delay region. Moreover, many of the neurons in these sequences were choice-selective, such that neurons preferentially fired depending upon whether the animal will ultimately choose to turn left or right (Fig. 1B-E).

Previous models of evidence accumulation do not produce the observed choice-selective sequences. This includes traditional circuit models that neurally instantiate the drift diffusion model (Fig. 2A-B), as well as bump-based integrators (Fig. 2C-D) commonly used to accumulate head velocity signals into a representation of head direction. Further, canonical sequence models such as synfire chains60–62 and more recent models of replay sequences63 do not accumulate evidence, while other sequence models accumulate external signals but are not choice-selective64,65. This motivates our construction of sequence-based evidence-accumulation models.

Competing Chains Models

We first present a class of models consisting of two chains of neurons. In these models, the current position of the animal is represented by the location of elevated neural activity along the chains and the accumulated evidence is represented by the relative amplitude of firing of the two neurons corresponding to this location. Here, we consider a model with mutual inhibition between neurons at the same position in the chains (mutually inhibiting chains model, Fig. 3A), but other competing chains models are possible (see below, Supplementary Fig. 2, and Supplementary Text). The temporal dynamics of the firing rate r of the ith neuron in the left (L) side chain are given by the equation

dri,L = −ar +[br +cr −er +f1 (t)+P(t)−T]+ (1),
i,L i, L i−1, L i,R left i
dt

and an analogous equation describes the firing rates r of neurons in the right (R) side chain. Here, a is the exponential decay rate of activity in the absence of input, b is the weight of self-excitation (Fig. 3A, green self-loops), c is the strength of feedforward synaptic connections from the previous neuron (green feedforward connections), and e is the strength of the inhibitory connection from the neuron at the same position in the opposite chain (brown connections). The external inputs to be integrated by the network consist of brief pulses of the form 1 (t), which is defined to be 1 if there is a left cue within left
0.5 cm of the current position and 0 otherwise, and project to the chain with synaptic strength f. P(t) is the position-gating signal, and T is the firing threshold. The position-gating signal P(t) controls the i
progression of activity along the chain by raising the neuron’s input above the firing threshold T, and is zero for all neurons except the pair of neurons (one from each chain) corresponding to the position of the animal in the environment. The firing rates of neurons not at the active position decay exponentially with time constant 1/a. Firing rates and inputs during a subset of an example trial are shown in Fig. 3B. See Methods for complete details and parameterization of the model.

At a given position, this model has the same structure as the traditional mutual inhibition model26 (Fig. 2A), and this motif is repeated across different positions to form two chains. This motif accumulates and maintains evidence through positive feedback within and between the two neural populations receiving the active position gating signal. Competition is mediated by mutually inhibitory positive feedback between the chains, such that increases of activity in one chain results in corresponding decreases of activity in the other chain. We note that the accumulation process saturates at an upper bound when neurons in the non-dominant chain reach zero firing. This is because, when one chain is at zero firing rate, the disinhibitory component of positive feedback between the chains is lost, so that further increases in the magnitude of evidence cannot be maintained. This saturation is not a core feature of the model, however, and the non-saturating range of evidence levels can be modified by the amount of external input to each chain or the synaptic strength of the visual input connections (Supplementary Fig. 2J). As in the case of traditional, non-sequential accumulator models, perfect accumulation up to this upper bound requires fine-tuning of the parameters to prevent leak or uncontrollable growth in the integrator35 (see Supplementary Text). Previous studies of such networks have also suggested that such accumulator networks may be made more robust to these fine-tuning conditions through incorporation of bistable components66,67 or corrective feedback mechanisms68,69.

Only the neurons receiving the position gating signal actively represent and accumulate evidence. When the animal moves to a new position, the position gating signal activates the next pair of neurons in the chain, and evidence is transferred to the newly activated pair of neurons by the feedforward synaptic connections64 (see Supplementary Text for conditions on the feedforward connectivity for faithful transfer of evidence between positions). Repeating this process as the animal navigates down the track results in a sequence of neurons firing, with each successive pair of activated neurons representing the accumulated evidence in the difference of their firing rates.

We simulated the model on the accumulating towers task (see Methods), demonstrating that this model performs the task perfectly in the absence of noise for appropriately tuned parameters (Fig. 3C, see Supplementary Text). Previous analysis of the accumulating towers task and other evidence accumulation tasks suggest that noise in an animal’s performance is dominated by sensory noise, so that only a fraction of input cues are integrated13,59. When we correspondingly introduce input noise into our model (see Methods), the performance decreases in a manner similar to that observed in animals, with a more gradual increase in performance with increasing magnitude of the difference between cues (Fig. 3C). The neurons in the model exhibit choice-selective (or, more precisely, accumulated-evidence-selective) sequential activity, with neurons in each chain showing greater activity for the side of the evidence that it receives (Fig. 3D). Thus, accumulated evidence is encoded monotonically (up to saturation) in the relative firing rates of the two chains in this class of models (Fig. 3E).

Variant architectures of competing chains models are possible, including an uncoupled competing chains model where competition between opposing inputs occurs without inhibition between the chains (Supplementary Fig. 2A, see Supplementary Text). At a given position, this model has a push-pull feedforward arrangement of evidence inputs, with cues favoring one side exciting one chain and inhibiting the opposing chain, as in traditional decision-making models based on uncoupled competing accumulators27. This motif is then repeated to form two chains. As with the model based on mutually inhibiting chains of neurons, this uncoupled competing chains model is able to perform the task accurately (Supplementary Fig. 2B) and produce choice-selective sequences (Supplementary Fig. 2C) in which individual neurons have monotonic tuning curves to evidence (Supplementary Fig. 2D). As in the mutually inhibitory chains model, this model can either saturate or not depending on the strength of the synaptic input connections and assumptions about the dynamic range of the neurons (Supplementary Fig. 2K).

Although these variant competing chains models cannot be reliably distinguished by their tuning since both models can be modified to have either linear or saturated tuning curves, the models instantiate different hypotheses that could be distinguished by functionally-targeted optogenetic perturbations70 (Fig. 3F, Supplementary Fig. 2E). We simulated the delivery of an excitatory input to a neuron (or similarly tuned neuronal population) while the animal is at the neuron’s active position (Fig. 3F, Supplementary Fig. 2E) for a trial with no external cues so that the accumulated evidence prior to perturbation is zero. Due to the gating of neural activity by the forward-traveling position signal, both models predict that only neurons tuned to later positions will be affected by the optogenetic perturbation. For the mutually inhibiting chains model, subsequent neurons in the same chain will have increased activity due to this feedforward excitation while subsequent neurons in the opposing chain will have decreased activity due to the mutual inhibition between the chains (Fig. 3G). By contrast, for the models based on uncoupled chains in which the anti-correlated activity of the two chains is solely due to opposing external cue inputs (Supplementary Fig. 2A), an excitatory perturbation in one chain increases firing only within the same chain but has no impact on the other chain (Supplementary Fig. 2F).

Another difference between the mutually inhibiting and uncoupled competing chains models is that the model based on mutual inhibition is more robust to the tuning of the position-gating. Specifically, the mutual inhibition model requires only that the position-gating signal P be greater than or equal to the threshold T at the active position (Supplementary Fig. 2G,H), whereas the uncoupled chains model requires that P be exactly equal to T to avoid continual growth of neural activity from temporally integrating the position-gating signal71 (see Supplementary Text).

Taken together, these competing chains models represent a class of models in which position-gating gives rise to choice-selective sequences that encode graded evidence signals monotonically in the neuronal firing rates at each position.

Position-Gated Bump Attractors

We next consider a class of models in which evidence is encoded in the location of a stereotypical, unimodal pattern (bump) of neural activity in the population. Such bump attractor models have been used to describe the neural circuitry that computes head direction in the rodent72–77 and the fly78–84, where the location of the bump corresponds to the animal’s heading. Similarly, two-dimensional bump attractor models have been proposed to describe path integration in the hippocampal place cell system85–87. In these models, motion is temporally integrated to determine head direction or to perform path integration, and analogously our models temporally integrate evidence from the visual cues. A traditional bump model that only accumulates evidence does not produce sequences like those observed in the neural data because the bump would recur at the same location whenever the accumulated evidence is at the same level and would remain at a fixed location throughout the delay period when the evidence level does not change (Fig. 2D). Thus, we modify these traditional bump attractor models to have separate dimensions for encoding position and evidence (Fig. 4 and Supplementary Fig. 3).

The model consists of layers of recurrently connected neurons, with each layer corresponding to a given position and neurons within the layers representing the accumulated evidence at this position. The activity r of a neuron at position i and evidence level j evolves according to

dri,j = −ar +F(∑W r +br +c(I (t)+I (t))+P(t)−T)(2)
dt i,j j,k i,k i−1,j i,j+1,L i,j−1,R i
k

where a is the exponential decay rate of activity in the absence of input, P(t) is the position-gating i
signal that selects the active layer of neurons, and T sets the soft threshold of the (1 + tanh(⋅)) neuronal output nonlinearity F. Within a layer, W is the strength of the connection from the neuron representing j,k
evidence level k onto the neuron encoding evidence level j. W is excitatory for connections of a j,k
neuron onto itself and its immediate neighbors (Fig. 4A, green recurrent connections) and inhibitory onto all other neurons in that layer (Fig. 4A, brown recurrent connections). The coefficient b is the synaptic strength of feedforward excitatory connections between layers of neurons encoding neighboring positions, and we assume for simplicity that these connections are only between neurons encoding the same evidence level (Fig. 4A, green feedforward connections). The terms I and I
i,j+1,L i,j-1,R
are the firing rates of the shifter neurons (analogous to the rotation cells in Skaggs et al. (1995)72) (Fig. 4A, orange and purple circles) that shift the bump left or right in response to external cue inputs (Fig. 4A, left and right tower inputs; Fig. 4B, shift in location of the bump following an input). The left shifter (orange) neuron is activated when jointly receiving input from evidence accumulating neuron j+1 and a left-side external cue. In this manner, if the peak of the bump of neuronal activity was previously at j+1, this peak activity then (through a synaptic connection of strength c from the shifter neuron) drives the neuron at evidence level j. Mathematically, the firing rates of the left, and corresponding right, shifter neurons are given by

I = r 1 (t) (3)
i,j,L i,j left
I = r 1 (t) (4),
i,j,R i,j right

where 1 (t) equals 1 if there is a left or right cue, respectively, within 0.5 cm of the current position. left/right
For the complete model specification, see Methods.

Within the active position, accumulated evidence is represented and maintained in the location of a bump of activity along the evidence axis. As in traditional bump attractor models88–90 (Fig. 2C), the shape of the bump is maintained in the absence of input by local symmetric excitation and long range inhibition. Asymmetric input connections cause the bump to shift left or right in response to a left or right cue, respectively73,77, as described above (see Supplementary Text). This shift of the localized bump of activity occurs despite the external cue signal being global (Fig. 4A, tower inputs). Mechanistically, the ability to shift a localized bump of activity with a global external signal is achieved via the multiplicative gating of the external input by the activity of the evidence-encoding neurons (Eqs. (3) and (4)) so that the global external input signal only activates the shifter neurons close to the location of the current bump. Alternatively, this gating could be accomplished with additive inputs plus thresholds84,88.

As the animal progresses forward in the maze, the active position changes so that information must be transferred to the next layer of neurons, requiring both a change in the position gate and feedforward connections between layers (Fig. 4B). Here, position-gating is controlled by the position signal P(t) and the soft threshold T, so that the output of F is only significantly greater than zero when i
the position signal is active (i.e., P(t) ≥ T). Away from the active position (i.e., where P(t) = 0), activity i i
decays exponentially with time constant 1/a (Fig. 4B). Feedforward excitation between the layers then transfers information to neurons at the next position. We have modeled these feedforward connections as only projecting between neurons representing the same evidence level, but this projection could be broader as long as the feedforward input is symmetric and centered around the current evidence level (see Supplementary Text). As in the chains-based models, we have presented the position signal as an external input. However, position could simultaneously be generated within the network through an integration of velocity inputs, by including a similar arrangement of excitatory and inhibitory connectivity along the position axis. In this case, the network forms a planar bump attractor that jointly encodes both evidence and position (Supplementary Fig. 3, see Supplementary Text).

We simulated the position-gated bump attractor model for the same set of trials as the competing chains models. As in the competing chains models, for appropriately tuned parameters, the position-gated bump attractor also performs the task perfectly in the absence of noise, with performance worsening when input noise is introduced (Fig. 4EC). The simulated neural activity also produces choice-selective sequences (Fig. 4D). As seen by the downward curvature of the sequences in Figure 4D, the sequences are marked by the feature that more neurons respond at late positions than respond at early positions in the sequence. This is because the random walk of the bump location as evidence is accumulated leads to a greater range of evidence levels being observed at later positions in the trial; thus, on average across trials, there are more different neurons active at later positions, with each neuron corresponding to a different possible evidence level. However, we note that this is not a fundamental prediction of the model because, experimentally, the exact shape of the choice-selective sequence depends on the number of neurons residing within each uniformly spaced evidence-position location in the model (Fig. 4A). We also note that our assumption of perfect symmetry simplifies the analysis, but is not required to produce bump attractor models91,92 and that the robustness of bump attractor models to diffusion and drift of the bump in the presence of noise has been extensively studied91,93–95 and may be reduced through error-correcting codes96, extending representations to multiple bumps97, or relaxing the continuous attractor assumption to have only a discrete set of memory states98.

The most straightforward difference between the position-gated bump attractor and the competing chains models is in the shape of their evidence tuning curves. Due to the structure of the network connectivity, the neurons in the position-gated bump attractor are non-monotonically tuned to specific evidence levels at their active position (Fig. 4E), in contrast to the monotonic tuning curves of the competing chains models (Fig. 3E).

To further differentiate the position-gated bump attractor model from the competing chains models, we consider predictions for the effects of optogenetic excitation of a single neuron during a trial with no cues presented (Fig. 4F). Here, we consider the regime where the stimulus is sufficiently large to cause the bump to move to the stimulated location. Due to the similar feedforward connectivity patterns between positions in both model classes, stimulation of a neuron active at one position only affects the activity of neurons tuned to later positions, as in the competing chains models. However, unlike the competing chains models, we see that rather than impacting every neuron at a later position equally, the effects are sparse and differ in magnitude across neurons (Fig. 4G). These differences are because the stimulation has caused the bump to move to a new evidence level, which is then sustained at subsequent positions. Thus, neurons with a similar evidence tuning to the neuron that was stimulated increase in activity, while those tuned to where the bump of activity would be located in the absence of stimulation (at 0 evidence) decrease their activity since the bump has moved to the new location (Fig. 4H). We note that for smaller excitatory stimuli that do not shift the bump center to the site of stimulation, and dependent upon the exact details of the excitatory and inhibitory connections, smaller shifts in the bump either toward or away from the stimulated location are possible95,99–101, but the perturbations still maintain the properties that changes are sparse and that changes within the same choice-preferring population as the stimulated neuron may be bidirectional.

Neural Activity in Different Brain Regions Correlates with Different Classes of Models

Above, we presented two circuit models for accumulating evidence through sequences (competing chains models: Fig. 3, Supplementary Fig. 2; position-gated bump attractor model: Fig. 4). While both models produce neurons that are sequentially tuned to position, the models predict different evidence coding schemes. The competing chains models make two core predictions: first, evidence should be monotonically encoded in neuronal activity and, second the monotonic encoding within the population should be graded across a range of evidence levels. The bump attractor model likewise makes two core predictions: first, evidence should be encoded with relatively narrow, unimodal tuning and, second, the peaks of these tuning curves should be broadly distributed across the entire evidence range. We looked for signatures of these different evidence encoding schemes in recordings from ACC, RSC, DMS and HPC during the accumulating towers task (Fig. 1A).

Towards this end, we first identified cells that were evidence-tuned, showing any consistent response pattern at different accumulated evidence values. For these evidence-tuned cells, we then identified whether neurons were monotonic-like or unimodal-like in their evidence tuning and also, for unimodal tuning, we characterized whether the peaks of the tuning curves were broadly distributed across the range of evidence levels (Supplementary Fig. 4). Within the class of monotonically evidence-tuned neurons, we next assessed whether the evidence tuning curves had broadly graded tuning to evidence, versus more binary tuning that would be consistent with a choice-like representation, and also how this distribution varied with position in the maze. We further examined whether graded evidence could be decoded from the population as a whole and whether the encoding of evidence across the population was consistent with the prolonged, temporal response kernels required for graded accumulation of evidence.

ACC and RSC Neurons Primarily Have Broad Monotonic Evidence Tuning While HPC Neurons Have Narrow Unimodal Evidence Tuning

Identifying neurons that were tuned to evidence requires jointly, rather than independently, fitting position and evidence to account for the correlated, nonuniform sampling of these two variables (Supplementary Fig. 4A-B). For example, neurons tuned to early positions in the maze are only active when accumulated evidence is small, given that little evidence has been presented at early positions. Thus, if such neurons were directly fit to evidence tuning alone, rather than jointly fit to evidence and position, they would appear to be tuned to small absolute evidence levels even if they had no tuning to evidence at all. Likewise, because the evidence levels sampled at later parts of the maze were most concentrated around moderate-to-large absolute values (Supplementary Fig. 4A), neurons untuned to evidence but tuned to later positions in the maze would artifactually appear to have a bimodal tuning to evidence.

To capture many different tuning curve shapes across the population, we first modeled the joint tuning for position and evidence as a product of a Gaussian function of position and a Gaussian function of evidence (Supplementary Fig. 4B, Fig. 5C,D,G,H,K,L, Supplementary Fig. 5C,D; see Methods). The Gaussian in position intuitively captures the non-monotonic, sequential position tuning. The Gaussian in evidence, although non-monotonic in its mathematical form, allowed us to identify not only non-monotonic but also monotonic evidence tuning, because in the latter case the best fit Gaussian has a peak near the extreme of the evidence range (i.e., normalized peaks near +/- 1; Supplementary Fig. 4D). Broad monotonic ramps like those of the competing chains models also have a large standard deviation relative to the observed evidence range.

Jointly fitting position and evidence tuning across the 14,247 neurons revealed that 56% of all neurons in ACC, 48% in RSC, 25% in HPC, and 41% in DMS had significant evidence tuning (showing any pattern of response to evidence with statistical comparison to pseudosessions; see Methods). In analyzing neurons with significant evidence tuning in each region, we saw distinguishing patterns emerge in the nature of the evidence tuning (Fig. 5, Supplementary Fig. 5).

In ACC, many individual neurons exhibited primarily monotonic and relatively broad evidence tuning (Fig. 5A-B, Supplementary Figs. 6A,7). To accurately visualize evidence tuning, in addition to plotting raw position and evidence fields (Fig. 5A), we also fit both a Gaussian and a logistic function to the peak position tuning (Supplementary Fig. 4C) for a number of example neurons, and plotted the best fit function (Fig. 5B,F,J). Similar to these example neurons, among the population of significantly evidence-tuned neurons, peak evidence tuning was towards the extreme of the observed evidence range, consistent with monotonic tuning (Fig. 5C-D). The presence of monotonic tuning is most consistent with predictions of the competing chains models (Fig. 3E, Supplementary Fig. 2D,J,K, Supplementary Fig. 4E,F), although the logistic fits revealed many neurons had steeper sigmoidal tuning than produced by the competing chains models (Fig. 3E, Supplementary Fig. 2D), especially at later positions in the maze (Fig. 5B, Supplementary Fig. 6A). The steepness of such sigmoidal tuning, and its variation across different positions of the maze, is treated further in the next section.

In contrast, HPC neurons showed primarily non-monotonic and narrow evidence tuning (Fig. 5E,F, Supplementary Figs. 6B,8), which tiled the range of observed evidence (Fig. 5F-H). These narrow tuning curves occurred throughout the cue and delay region (dark blue neurons in Fig. 5G). The presence of non-monotonic, narrow tuning curves that tile evidence space is consistent with predictions of the bump model (Fig. 4E, Supplementary Fig. 4G,H). In addition to narrow tuning curves, the bump model also predicts that the population should track the accumulated evidence on individual trials. Such narrow tuning curves make decoding challenging because, to perfectly decode evidence, recordings would need to include cells tuned to each combination of position, evidence, and other variables previously identified in this task49, and the HPC recordings are sparse. However, previous analysis of these data suggests that evidence is represented in the population on a nonlinear neural manifold in HPC49, and we similarly find that evidence can be decoded above chance using another nonlinear embedding technique102 (Supplementary Fig. 11).

In RSC, the majority of neurons showed monotonic patterns (Fig. 5I-J, Supplementary Figs. 6C,9) like those observed in ACC, while a smaller population showed more narrow tuning curves like, but typically somewhat wider than, those observed in HPC (example 2 in Fig. 5I, examples 3, 9, 16, 22, 28, and 34 in Supplementary Fig. 9). These observations are seen in the population summary of evidence and position tuning (Fig. 5K-L), which shows that many neurons have wide evidence tuning with peaks at the extreme (deep red) while many others have narrow intermediate evidence tuning (dark blue). RSC also has a particularly large population of neurons that, when fit to a Gaussian, have intermediate widths and evidence peaks not fully at the extremes (Fig. 5K, green); such neurons typically exhibit a steep sigmoidal tuning curve such as that shown in example 3.

DMS was dominated by monotonic evidence tuning, but with some narrowly tuned neurons (Supplementary Fig. 5, Supplementary Figs. 6D,10). However, DMS contained relatively few evidence-modulated neurons with peak position tuning in the cue region. Instead, most evidence-modulated neurons had their peak position tuning in the delay region and had choice-like steep sigmoidal evidence tuning (Supplementary Fig. 5C,E).

Shift from Evidence Tuning in the Cue Period to Choice Tuning in the Delay Period in ACC and RSC

While the monotonicity of the evidence tuning of individual cells observed in ACC and RSC is consistent with the competing chains models, a candidate integrator must represent graded evidence beyond monotonic, step-like evidence tuning driven by choice. Since the units of our models should not be thought of as individual cells, but rather populations of neurons with similar position tuning, the graded representation predicted within our model units could emerge from either individual cells having graded representations or individual cells stepping on at different evidence thresholds such that population averages are graded103. Therefore, we examined whether the regions dominated by monotonic tuning for evidence contained representations of evidence in the population beyond binary choice tuning, using both single-cell and population-level analyses.

We first asked whether there were signatures of evidence accumulation within the population by deriving cue response kernels from a linear cue encoding model fit to the difference in right- and left-preferring population activity (Fig. 6A-C, see Methods). As required for integration, we found cue-locked, sustained, and relatively constant steplike responses to the cues, with left and right cues having a similar pattern but with opposite signs (Fig. 6D,E). When accounting for the amount of evidence accumulated so far, the magnitude of the cue responses decreased, consistent with a saturation at more extreme evidence levels an d a more choice-like representation as the maze progresses (Supplementary Fig. 13A,B). Interestingly, the cue response kernels decreased less as a function of position in the maze than as a function of accumulated evidence (Supplementary Figure 13C,D; note nonzero kernels even for later positions), suggesting that the transition observed as a function of previously accumulated evidence may be due to animals committing to a decision with larger levels of evidence104,105, rather than solely being due to being at a later position in the maze. Most importantly, this analysis suggests that accumulated evidence is represented beyond choice at the population level (at least for intermediate values of evidence).

We next sought to quantify the mixture of graded evidence and binary choice activity in individual cells. To do so, we fit single-neuron evidence vs. choice encoding models at each position bin, using evidence and behavioral choice as regressors for both correct and incorrect trials (Fig. 6F, see Methods). At each position, we calculated the fraction of active neurons (see Methods) whose firing rate was significantly modulated by the linear evidence or the binary choice variable, while accounting for the contribution to neural activity of the other variable. In both ACC and RSC, we observed neurons with significant encoding of evidence beyond that explained by choice during the cue period, with decreasing proportions of cells encoding evidence during the delay period (crimson traces, Fig. 6G,H). In contrast, the fraction of significant choice encoding neurons (beyond that explained by the evidence) increases throughout the maze, with large fractions of active neurons encoding choice in the delay period (lime traces, Fig. 6G,H). The above analysis primarily differentiates behavioral choice-tuned neurons from evidence-tuned neurons by their change in tuning to evidence between correct and incorrect trials. However, a subtlety is that a neuron that has binary tuning to the current best choice (i.e., the current sign of the evidence), rather than graded tuning to evidence, may be identified as evidence tuned despite not being graded in its response. Therefore, to test whether there was graded response in the population of cells identified as evidence-tuned, we plotted the trial-averaged firing of the active cells in this subpopulation and confirmed that this subpopulation exhibited graded patterns of activity with evidence during the cue period, consistent with our model predictions for accumulation (Supplementary Fig. 13E,F).

To complement this single-neuron encoding analysis, we also asked whether evidence could be linearly decoded from the population (Fig. 6I). Similar to the findings from the single-neuron analysis, evidence can be linearly decoded from the population of active cells above chance early in the maze, falling to chance levels in the delay region (Fig. 6J,K). This is true even when controlling for the sign-of-evidence and considering only nonzero evidence levels to eliminate the possibility of a binary encoding of just the current sign-of-evidence (Supplementary Fig. 13G,H). Together, these results imply that graded evidence representations in the cue region transition to a more binary choice representation in the delay period.

We sought to summarize these results with a population level visualization. The population-averaged position-evidence responses in both ACC and RSC show graded increases in activity with preferred-evidence (i.e., evidence calculated in reference to each neuron’s preferred direction of evidence) during the early portion of the cue period (Fig. 6L,M, dark blue traces; population responses based on averaging all evidence-modulated neurons, see Methods). Consistent with the other analyses (Fig. 6G,H,J,K), with progression along the maze, this evidence tuning sharpens, transitioning to a step function by the start of the delay period (magenta and purple traces) that resembles tuning to choice rather than graded evidence tuning. In ACC (and more subtly in RSC), population responses also increase later in the maze, due to the fact that there are more evidence-tuned neurons in the delay region that contribute to these population averages (Supplementary Fig. 12A,B, left panels).

Augmenting the Circuit Models to Recapitulate the Transition from Graded Evidence to Binary Choice Tuning

Neither of the competing chains models produce this transition from relatively more evidence-encoding to relatively more choice-encoding, since the models were designed to perfectly accumulate evidence. This transition could be realized in the models in multiple manners. First, a population of choice readout cells at each position could compare the accumulation of the left and right chains at that position and respond in a binary fashion when either the left or right chain has greater activity, where the number of readout cells at each position increases over the course of the maze (Fig. 7 A-C). Second, this pattern could be realized if the parameters of the mutually inhibiting chains model are slightly mistuned such that the integrator becomes increasingly unstable later in the maze (Fig. 7 D-F). Such slight instability has been suggested in previous integrator models, where the instability has been functionally interpreted to provide an urgency signal that accelerates decision commitment106,107. These mechanisms are not mutually exclusive. Thus, the prevalence of choice cells in the data complements the integration processes focused upon in our models by suggesting that there are additional processes at play driving the transition from evidence accumulation to choice representation.

Methods
The Accumulating Towers Task

The accumulating towers task is a behavioral paradigm, performed in virtual reality, for studying accumulation-of-evidence based perceptual decision making. The complete details of the task and experimental setup can be found in Pinto et al. (2018)59. Here, we briefly summarize the key features of the task. Head-fixed mice run on an 8-inch styrofoam ball within a virtual environment, generated by ViRMEn143, projected onto a surrounding screen. The virtual environment has a T-maze structure. In the first 30 cm (the precue region), there are no visual cues. In the next 200 cm (cue region), visual cues appear as towers on both sides, with the total count of towers on each side determined by a Poisson distribution of different means. Towers are positioned randomly with an enforced minimal distance between them. A tower appears when the mouse is within 10 cm of its position and disappears 200 ms later. Following the cue region, there is a 100 cm delay region during which the mouse must remember which side was presented with more towers. At the end of the delay region, the mouse receives a sweetened liquid reward for turning to the side with more cues. For water restriction and shaping protocols, see Pinto et al. (2018)59. Recording sessions were approximately 1.0-1.5 h in duration, consisting of approximately 150-250 trials per session.

In addition to the accumulation of evidence trials described above, sessions also included some warm-up maze trials (cues on one side only, no distractors, tall visual guides at the end of the maze) at the beginning of the session and interspersed in the session if the animal’s performance dropped below some threshold, as described in Pinto et al. (2018)59. Since these trials do not require evidence accumulation, we did not include them in our analyses.

Neural Recordings

Animal procedures were conducted in accordance with standards from the National Institutes of Health and under the approval of the Princeton University Institutional Animal Care and Use Committee.

The head-fixed setup of the accumulating towers task allows for neural imaging throughout the task. Our neural datasets consist of calcium fluorescence data previously recorded from retrosplenial cortex (RSC, 8 animals, 41 sessions (29 sessions from layers 2/3 and 12 sessions from layer 5), 8579 total neurons)44 and hippocampus (HPC, 7 animals, 7 sessions, 3144 total neurons)49, as well as newly acquired recordings from anterior cingulate cortex (ACC, 3 animals, 7 sessions, 1720 total neurons) and dorsomedial striatum (DMS, 8 animals, 11 sessions, 804 total neurons). See below for details of ACC and DMS recordings.

For DMS and ACC data, mice were co-housed with same-sex littermates prior to and following surgical procedures. Mice were maintained on a 12-hr/12-hr light/dark cycle. All surgical procedures and behavioral training occurred in the dark cycle. Ambient temperature was maintained at 21–26 °C and humidity at 30–70%.

Animals and Surgery for ACC Recordings

Three female mice of CaMKIIa-tTA (JAX 007004) crossed with tetO-GCaMP6s144 (JAX 024742) were used for two-photon imaging of ACC. These mice, which were heterozygous for each transgene, showed robust and stable expression of GCaMP6s in a subset of CaMKIIa-positive, glutamatergic neurons in the cortex. All mice were between 12-16 weeks when they underwent sterile stereotaxic surgery under isoflurane anesthesia (4% for induction, 1-1.5% for maintenance). The skull was exposed by making a vertical cut with a scalpel. Two craniotomy holes were made, one over ACC (anteroposterior axis: 0.8 mm, mediolateral axis: -0.5 mm, relative to bregma, diameter of ~1.2 mm) and another over DMS (anteroposterior axis: 0.75 mm, mediolateral axis: -1.1 mm, relative to bregma, diameter of ~0.3 mm). Adeno-associated virus (AAV; AAVretro2 serotype) encoding mCherry under CaMKIIa promoter (diluted to 2.5 x 1012 genome copies/mL, packaged at Princeton Vector Core) was injected to the DMS (anteroposterior axis: 0.75 mm, mediolateral axis: -1.1 mm, dorsoventral axis: -3.1 and -2.8 mm, relative to bregma) to provide information about a projection target and to be used for better motion correction. A total of 300 nL of AAV was infused at each site along the dorsoventral axis, at a rate of 30 nL/min. After injection, the needle was held in the same place for an additional 10 minutes. The needle was slowly withdrawn over 10 minutes to prevent backflow. After injections, the dura below the craniotomy was carefully removed within the ACC craniotomy site, and the cortical tissue above ACC was carefully aspirated up to 300 𝜇m below the surface. Then we implanted a 1.0-mm diameter GRIN lens (4.3 mm length, Inscopix, 1050-002184) into the ACC, down to 1.5 mm in the dorsoventral axis (relative to bregma), using a 3D-printed custom lens holder. After implant, a small amount of diluted metabond (Parkell) was applied to affix the lens to the skull using a 1-mL syringe and 25-gauge needle. After 10 minutes, the lens holder grip was loosened while the lens was observed through the stereoscope to make certain that there was no movement of the lens. After applying a more diluted metabond around the lens, a titanium headplate was positioned over the skull using a custom tool and aligned parallel to the stereotaxic frame using an angle meter. The headplate was affixed to the skull using a metabond. Finally, a titanium ring to hold the immersion medium was placed on the headplate, with lens in the middle, and then affixed to the headplate using dental cement blackened with carbon. After the dental cement was fully cured, mice were unmounted from the stereotaxic frame and protective head caps were placed. Recovery of the mice was monitored for ~2 hours. After one or two weeks of recovery, mice started water-restriction, and after 5 days, behavioral training for the accumulating towers task began, as described in Pinto et al. (2018)59. Histology showed that GRIN lenses sampled mostly deep layers (layers 5 and 6) in two animals, while in one mouse a GRIN lens was medially implanted to cover both superficial (layers 2 and 3) and deep layers. For the data analysis, we did not make distinctions between neurons in different layers but instead pooled neurons from all layers together.

Animals and Surgery for DMS Recordings

D1R-Cre (n = 3, female, EY262Gsat, MMRRC-UCD), A2a-Cre (n = 3, male, KG139Gsat, MMRRC-UCD), and D2R-Cre (n = 2, male, ER44Gsat, MMRRC-UCD) mice were backcrossed to a C57BL/6J background (Jackson Laboratory, 000664) and bred in-house. All transgenic mice were heterozygous for the transgene. These mice were used for DMS two-photon imaging. Surgical procedures were identical to those described above for the ACC except for the following details. Approximately one week prior to GRIN lens implantation, a Cre-dependent AAV expressing GCaMP6f under either the CAG promoter (AAV5-CAG-DIO-GCaMP6f, Addgene, diluted to 3.5 x 1012 genome copies/mL) or synapsin promoter (AAV1-Syn-DIO-GCaMP6f, diluted to 3.8 x 1012 genome copies/mL) was delivered (500ul per injection site; 200 ul/min) to the DMS (anteroposterior axis: 0.7 mm, mediolateral axis: 1.15 and 1.85 mm, dorsoventral axis: -3 mm, relative to bregma) and the surgical incision was closed using monofilament non-absorbable sutures (Johnson and Johnson). This allowed for expression of GCaMP6f in D1R- or D2R-expressing medium spiny neurons according to the transgenic mouse line used, and facilitated the settling of brain tissue prior to GRIN lens implantation. Following recovery (5-10 days), sutures were removed and the viral injection craniotomy was expanded and dura removed to accommodate a 1.0-mm diameter GRIN lens with a prism (4.3 mm length, Inscopix, 1050-004601). The lens was implanted with prism window facing posterior and just anterior to viral injection sites (anteroposterior axis: 0.75 mm, mediolateral axis: 1.35 mm, dorsoventral axis: -3.25 mm, relative to bregma). Fixation of the lens to the skull, as well as the custom-printed headplate, water immersion ring, and protective covering were identical to that described above for ACC. Across all DMS targeted GCaMP6f imaging mice we obtained three distinct fields of view from one mouse, two distinct fields of view from a second mouse, and the remaining mice contributed single fields of view (8 animals, 11 sessions, 804 total neurons).

Optical Imaging and Data Acquisition for ACC and DMS Recordings

Imaging was performed using a custom-built two-photon microscope compatible with virtual reality (VR). The microscope was equipped with two high-power femtosecond lasers (Alcor, Spark Laser), with fixed wavelengths for 920 and 1064 nm. The scanning units used a 5-mm galvanometer and a 8kHz resonant scanning mirror (Novanta). The collected photons were passed to a dichroic mirror (XYZ) that separates near-infrared light from the visible light and then split into two channels by another dichroic mirror (FF562-Di03, Semrock). The light for respective green and red channels was filtered using a bandpass filter (FF01-520/60 and FF02-607/70, Semrock), and then detected using GaAsP photomultiplier tubes (1077PA-40, Hamamatsu). The signal from the photomultiplier tubes was amplified using a high-speed current amplifier (59-179, Edmund). Output of the amplifier was sent to a data acquisition system (National Instruments). Black rubber tubing was attached around the objective (Zeiss, W N-Achroplan 20x/0.5 M27, 20x magnification, 0.5 NA) to shield the light from the virtual reality screen. The objective was aligned to the implanted GRIN lens within the implanted titanium headring. Double-distilled water was filled inside the headring and used as the immersion medium. The amount of laser power at the objective was 30~40 mW for ACC and 40~50 mW for DMS. Control of the microscope and imaging acquisition was performed with the ScanImage software (MBF Bioscience), on a separate computer from the VR-running one. Images were acquired at 30 Hz with a size of 512 x 512 pixels. Synchronization between behavioral logs and acquired images was achieved by sending behavioral information from the VTR controlling computer to the scanning computer via an I2C serial bus each time the ViRMEn environment was refreshed. Specifically, numeric information about the current session, trial, and frame were written in the header of TIFF imaging files.

Data Processing for ACC and DMS Recordings

Both ACC and DMS imaging data were processed with Suite2p. For ACC recordings in which mCherry signals are present, rigid motion correction was performed with that channel and then applied to the functional green channel with GCaMP6s, with the following parameter settings: nimg_init = 1000, smooth_sigma = 1.15, maxregshift = 0.1. For DMS recordings, the same motion correction was performed, but with the green channel. Regions-of-interest (ROIs) corresponding to individual neurons were detected, and calcium signals were extracted and neuropil-corrected (neucoeff = 0.7). Subsequently, an experimenter-mediated curation was performed on these ROIs to discriminate putative neurons from noise. Neuropil-corrected fluorescence from these curated neurons was exported to MATLAB, and the ∆F/F was calculated, taking baseline fluorescence as the mode of the fluorescence within 3 minute long windows (the same baseline as applied in the RSC and HPC datasets).

Spike Inference

In order to compare neural recordings across datasets with different calcium indicators, we performed spike inference on our data based on a recently developed and validated algorithm that has only two free parameters, the average estimated firing rate and the calcium decay rate145,146. For each ∆F/F dataset, we assumed a baseline firing rate of 2 Hz (but our results were not sensitively dependent on this specific choice) and set the calcium decay rate according to the values fit by Jewell et al. (2020)145 for each calcium indicator (0.9885 for GCaMP6s (ACC), 0.9704588 for GCaMP6f (all other regions)). From these inferred spikes, we inferred a continuous firing rate by smoothing with a gaussian window function (window length = 1 s, standard deviation = 0.25 s).

Position Binning of Neural and Evidence Data

To compare neural data along maze positions and evidence levels, we resampled the neural (continuous firing rate as calculated above) and behavioral (cumulative evidence up to that point) data in position bins. Specifically, firing rates were averaged based on the position of the animal within 5 cm bins from the start of the trial (-30 cm) to the end of the delay region (300 cm) for each trial. Cumulative evidence for each bin was defined by the number of right minus the number of left cues observed before the start of the position bin.

Trial Selection Criteria

We discarded a small fraction of trials from analysis using the following criteria. As mentioned above, any trials from the warm-up maze were excluded, as these trials do not require evidence accumulation. Additionally, trials with excessive travel were discarded (3.2% of trials in ACC sessions, 8.0% of trials in DMS sessions), which were identified as trials where mice showed a large view angle deviation (more than 90 degrees, either left or right) before the end of the delay region.

Determining Choice Selectivity of Neurons

To identify choice-selective neurons (Fig. 1B-E), we first jointly fit the position and evidence tuning of each neuron in our dataset as described in Data Analysis below. We then considered only neurons that had significant evidence tuning. To produce the sequence plots, neurons were first sorted as right-preferring if the mean of the neuron’s evidence tuning was greater than 0 and left-preferring otherwise and then sorted by their position mean. Using this sorting, we then took the average activity of each neuron at different positions for correct left and correct right trials. The resulting traces for each neuron were then peak normalized.

Traditional Integrator Models Simulations

To show the activity of traditional integrator models (Fig. 2), we simulated equations corresponding to previously developed models of integration: a competing accumulator model and a bump attractor model.

Competing Accumulator Model

We modeled two competing populations of neurons, one left-preferring population and one right-preferring population. The activity r of the neuron representing the left population is given by

drL = −ar +[br −er +f1 (t)+X(t)]+
L L R left
dt

and a corresponding equation governs the right population, where

1 min|p(t)− ℓ| < 0.5 cm
1 (t) = { ℓ∈L
left
0 otherwise

and

1 min|p(t)− R| < 0.5 cm
1 (t) = { R∈ℛ
right
0 otherwise

where p(t) is the position of the animal at time t, L is the set of positions of left cues, and ℛ is the set of positions of right cues.

In our simulations, we used parameters a = 50 s-1, b = 10 s-1, e = 40 s-1, f = 100 s-1, and, common external input to both populations, X(t) = 500 Hz/s. The model was simulated as described in Simulation of Differential Equations.

Bump Attractor Model

For the bump attractor model, we assumed each neuron has its peak response at a different evidence level. The activity of the neuron with peak response at evidence j is given by

drj
= −ar +F(∑ W r +c(I (t)+I (t))+X(t)),
j k j,k k j+1,L j−1,R
dt

where

q(1+tanh(x))
F(x) = ,
2

I (t) = r1 (t),
j,L j left

and

I (t) = r1 (t)
j,R j right

where

1 min|p( t)− ℓ| < 0.5 cm
1 (t) = { ℓ∈L
left
0 otherwise

and

1 min|p(t)− R| < 0.5 cm
1 (t) = { R∈ℛ
right
0 otherwise

where p(t) is the position of the animal at time t, L is the set of positions of left cues, and ℛ is the set of positions of right cues.

We assigned the neuron at each evidence level an angle, evenly spaced between 0 and 2π, and set the synaptic connection weights as

W = ω (cos(θ −θ )+ω ).
j,k 0 j k 1

We note that, since evidence is not periodic, a sufficiently large number of evidence levels were used such that substantial interactions do not occur between the nominally periodically connected endpoints. For the shown simulation, we used 35 neurons, corresponding to evidence levels of -17 to 17. We used the parameters a = 55 s-1, X(t)=0.04 Hz/s, q = 1250Hz/s, ω = 0.12 s-1, ω = -1.0, and c = 0.2 s-1. Neurons r , r , and r were initialized to 15 Hz, 17.5 Hz, and 15 Hz activity, respectively, and all other neurons were initialized to 0 Hz activity. The model was simulated as described in Simulation of Differential Equations.

Parameterization of the Mutually Inhibiting Chains Model

As in Eq. (1), the activity of the neuron at position i in the left chain can be described as

dri,L = −ar +[br +cr −er +f1 (t)+P(t)−T]+.
i,L i, L i−1, L i,R left i
dt

A neuron’s active position is defined by a square pulse over a range of positions, where p(t) is the position of the animal at time t,

p(t)
T+X(t) i ≤ < (i+1)
P i (t) = { P0 .
0 otherwise

The external input is of the form

1 min|p(t)− ℓ| < 0.5 cm
1 (t) = { ℓ∈L
left
0 otherwise

and

1 min|p(t)− R| < 0.5 cm
1 (t) = { R∈ℛ
right
0 otherwise

where L is the set of positions of left cues, and ℛ is the set of positions of right cues.

In our simulations, we used parameters a = 50 s-1, b = 10 s-1, c = 50 s-1, e = 40 s-1, f = 100 s-1, T = 15000 Hz/s, and X(t) = 500 Hz/s. (See Supplementary Text for necessary conditions on the parameters.) We take P = 20 cm, which divides the total length of the maze (330 cm) into 17 different sections of active positions. Therefore, in these models we simulate 17 neurons in each chain for a total of 34 neurons. This choice of parameters was selected to recapitulate typical firing rates of the neurons that we recorded as well as approximate the typical position width observed in the experimental choice-selective sequence plots (Fig. 1).

The parameters above produce saturation in the evidence tuning at the extreme evidence levels (Fig. 3E). As demonstrated in Supplementary Fig. 2J, for different parameters, this model will instead show linear tuning curves over the observed range of evidence (same parameters as above, but X(t) = 1000 Hz/s). We note that a similar effect of linear evidence tuning throughout the observed evidence range could have been achieved by decreasing the weight f on the cues. Both manipulations prevent the non-dominant chain from reaching 0 Hz, which results in saturation of the dominant chain.

In Supplementary Fig. 2G, we changed the form of the position gate to be a smooth gaussian, rather than a square pulse. To do so, we used the original parameters but took T = 1100 Hz/s and

−(p(t)−μ )2
i
P(t) = αexp( )
i β

where α = 60 Hz/s, β = 1250 cm2, and μ = 20i - 10 cm.
i

In Supplementary Fig. 2HT, we change the form of the position gate to have heterogeneous widths. We use the original parameters but only 12 neurons in each chain and rather than a fixed width P for each position, we draw a set of transition points ρ uniformly between -30 cm and 300 cm. If (ρ - ρ) < 5 cm, we take ρ = ρ + 5 cm, to enforce a minimum position width of 5 cm. We then define

T+X(t) ρ ≤ p(t) < ρ
P(t) = { i i+1.
i
0 otherwise

The Mutually Inhibiting Chains Model for the Transition to Choice Representations

In Figure 7, we propose two models that could explain the transition from more graded evidence representations early in the maze to more choice-like representations later in the maze. The first model (Fig. 7A-C) incorporates a population of choice-readout cells. In this model, the evidence accumulator cells of the chains are modeled exactly as described in the previous section, but we introduce a population of choice readout cells at each position. The firing rate of a left choice readout cell at position i, r , is given by
i,L,readout

r = α(sgn(r −r +ξ)+1)
i,L,readout i,L i,R

where α = 5 Hz, sgn is the signum function which returns the sign of its input, and ξ ~ 𝒩(0, 2) Hz with the analogous equation for a right choice readout cell. We assume that there are i choice readout cells for each choice at position i.

In the second model (Fig. 7D-F), we set the parameters of the mutually inhibiting chain model such that the integrator becomes increasingly more unstable with advancing position down the chain. Specifically, we incrementally increase the weight of self-excitation b in Eq. 1 with progression down the chain by making it a function of the position i:

b = b +Δi
i 0

where Δ = 0.5 s-1, and b = 10 s-1 is the value of self-excitation that gives a perfect integrator.
0

Parameterization of the Position-Gated Bump Attractor

In Eq. (2), we defined the evolution of the firing rate for a neuron at position i and evidence j as

dri,j
= −ar +F(∑ W r +br +c(I (t)+I (t))+P(t)−T).
dt i,j k j,k i,k i−1,j i,j+1,L i,j−1,R i

A neuron’s active position is defined by a square pulse over a range of positions, where p(t) is the position of the animal at time t,

p(t)
T+X(t) i ≤ < (i+1)
P i (t) = { P0 .
0 otherwise

The neuronal nonlinearity is

q(1+tanh(x))
F(x) = .
2

The strengths of the synaptic connections between neurons were chosen to be symmetric with excitatory connections between neurons nearby in their evidence tuning and inhibitory connections between neurons farther apart in their evidence tuning . This was done by assigning each evidence level an angle, evenly spaced between 0 and 2π, and setting

W = ω (cos(θ −θ )+ω ).
j,k 0 j k 1

We note that, since evidence is not periodic, a sufficiently large number of evidence levels were used such that substantial interactions do not occur between the nominally periodically connected endpoints.

When an input is present, asymmetric connections cause the active evidence neurons to shift towards the input. This is modeled through shifter neurons (Fig. 4A, purple and orange small circles) that are activated by the coincidence of external cue signals and neural activity from the current location of the bump by

I (t) = r 1 (t),
i,j,L i,j left

and

I (t) = r 1 (t)
i,j,R i,j right

where

1 min|p(t)− ℓ| < 0.5 cm
1 (t) = { ℓ∈L
left
0 otherwise

and

1 min|p(t)− R| < 0.5 cm
1 (t) = { R∈ℛ
right
0 otherwise

where p(t) is the position of the animal at time t, L is the set of positions of left cues and ℛ is the set of positions of right cues.

As in the competing chains models, the position-gating of each layer of neurons is activated over a length 20 cm (P = 20), again giving 17 different active positions. For the neurons to span the range of evidence of the experiments, we used 35 neurons in the evidence dimension, corresponding to evidence levels of -17 to 17. Thus we simulated a total of 595 neurons.

The parameters used in the simulation are: a = 55 s-1, b = 0.088 s-1, T = 300 Hz/s, X(t)=0.04 Hz/s, q = 1250 Hz/s, ω = 0.12 s-1, ω = -1.0, and c = 0.2 s-1. Neurons r , r , and r are initialized to 15 Hz, 17.5 Hz, and 15 Hz activity, respectively, and all other neurons are initialized to 0 Hz activity.

Simulation of Differential Equations

Differential equations for each model (Eqs. (1), (2), and (S2)) were simulated in Python using scipy.odeint, which uses the lsoda algorithm, with a maximum stepsize of 0.01 s (or equivalently 0.5 cm when assuming a constant velocity of 50 cm/s) to find the activity levels of the neurons in each model for points between -30 and 300 cm in 0.1 cm increments. In all simulations, we assume the animal travels at a constant velocity of 50 cm/s through the maze so that time and position are interchangeable. We simulated the neural activity for 1000 trials with cue positions derived from 1000 experimental trials from the behavioral and neural recording data. For the competing chains models, we say that the animal makes a left choice if the activity of the last neuron in the left chain is greater than the activity of the last neuron in the right chain. For the bump model, we say that the animal makes a left choice if the index k of the neuron with maximal firing in the last layer is less than 0.5J, where J is the total number of neurons in the layer. If the maximal firing occurs exactly at 0.5J, we determine the animal’s decision by randomly choosing between the two alternatives with equal probability.

Simulations with Input Noise

For each model, we simulated 25 sessions with input noise (Figs. 3C, 4C, Supplementary Fig. 2B). For each session, we selected 150 trials drawn randomly without replacement from the 1000 experimental trials that we used for simulation of the model. Let L and R be the set of left cues and right cues, respectively, in the original trial. For each cue x in L, we drew z, a random uniform number between 0 and 1, and if z<0.33, we added x to the set L’. We repeated the same process for each cue in R, yielding two reduced sets of cues when input noise is present, L’ and R’, corresponding to the animal ignoring or not perceiving 67% of the cues. We then simulated the model with the input terms 1 (t) and 1 (t) determined according to the reduced cue sets. We then calculated the psychometric left right
curves (as described in Psychometric Curves below) with the difference in cues based on the original cue sets, but whether the animal turns left or right based on the simulation of the trial with the reduced cue sets.

Evidence Tuning Curves for Simulated Model Data

In order to find the evidence tuning curves (Figs. 3E, 4E, Supplementary Fig. 2D,J,K), we averaged the activity of each neuron binned by position and accumulated evidence for each of the 1000 simulated trials. Position bins were of size 0.1 cm, ranging from -30 to 300 cm. Evidence bins were of size 1 tower, ranging from -15 to 15 towers. Bins that had fewer than 2 samples were not included in further analysis and not plotted. For a given neuron, we averaged across evidence levels to find the position with maximum firing. We then plotted the average firing rate versus evidence level at this position.

Simulated Single-Neuron Perturbations

Single-neuron perturbations were simulated in the same manner as no-perturbation trials but in the absence of any cues. A subthreshold excitatory term O(t) was added within the thresholded i
dynamics of each model (i.e. within the square brackets of Eq. (1) or within the function F of Eq. (2)) to simulate the effects of stimulating an individual neuron, where O(t) = 25 Hz/s for the mutually inhibiting i
competing chains model, O(t) = 12.5 Hz/s for the uncoupled competing changes model, and O(t) = 0.2 i i
Hz/s for the position-gated bump attractor, if neuron i is the neuron being stimulated, and 0 otherwise.

Psychometric Curves

For each session, we sorted trials based on the difference in the final number of cues at the end of the maze. We binned these trial differences with bin size of 3 towers, ranging from -14 to 14 towers. Within each bin, we then found the fraction of trials in each bin for which a left choice was made (Fig. 3C, Fig. 4C, Supplementary Fig. 2B).

Data Analysis

Determining Position and Evidence Selectivity

In our neural recording data, we sought to identify the set of neurons whose firing encoded accumulated evidence. As described in the main text, due to the nonuniform sampling of position and evidence and fact that neurons are only active at some positions, we found that robustly identifying evidence selectivity required that we jointly identify the position range over which a neuron is active and the evidence selectivity within this set of positions. To do so, we assumed the firing rate of each neuron had the form,

FR(p,e) = aP(p)E(e)+b

where FR(p,e) is the firing rate of the neuron at position p and evidence e, and P and E were described by Gaussians,

2
−(pP−μ )
p
P(p) = exp( )
2σ2
p

and

−(e−μe )2
E(e) = exp( ).
2σ2
e

The above functional form, which assumes evidence and position tuning can be approximated by independent Gaussians, was not used as a detailed descriptor of neuronal firing (for less constrained fits around a given position, see Neural Evidence Tuning Curves Around the Most Active Position). Rather, we found that it was effective for the purpose of identifying evidence-modulated neurons and determining if the mean evidence selectivity favors left or right evidence (as used to sort neurons by evidence-selectivity in Fig. 1B-E). Here, the Gaussian was allowed to extend beyond the range of observed evidence values, so that it can accommodate both monotonic and non-monotonic tuning curves (i.e., if the mean evidence occurs at the extremes of the evidence range, the evidence tuning will be monotonic throughout the range, Supplementary Fig. 4D).

For each neuron, we normalized the firing rate observations (binned as in Position Binning of Neural and Evidence Data) to the range [0, 1]. We then fit the parameters of the above model using an iterative fitting procedure, alternating fitting μ , σ , a, and b with fitting μ , σ , a, and b with a nonlinear p p e e
fitting method (scipy curvefit, which employs the trust region reflective algorithm, with bounds -50 ≤ μ ≤ p
350, 0 ≤ σ ≤ 200, e - 1 ≤ μ ≤ e + 1, 0 ≤ σ ≤ 30, 0 ≤ a ≤ 10, and 0 ≤ b ≤ 1). The bounds on μ are p low e high e e
defined by e = min(E’) and e = max(E’), where E’ is the set of all evidence levels observed by the low high
neuron within 5 cm of μ . These bounds are updated on each iteration based on the current fit of μ . p p
Due to the iterative nature of our fitting procedure, after the final position update, μ may no longer e
satisfy these bounds, and we do not include such neurons when analyzing the distributions of the normalized evidence parameters.

We then tested for significance of the evidence tuning fit using the pseudosession method147. For each session of firing rate and position data, we generated new evidence for each trial in the session by resampling trials from actual sessions, and applied the iterative fitting procedure to the new pseudosession evidence and position data. We repeated this process for 50 pseudosessions and calculated the mean squared error between the true firing rates and the model predictions for each session. We then used a one sample t-test to compare the fit of the model, as measured by the mean squared error, in the true session to the fit of the null distribution from the pseudosession, and only considered neurons for which the fit was significantly better (p<0.05) than the null distribution for further analysis.

In order to compare the parameterization of the evidence tuning across positions (Fig. 5D,H,L, Supplementary Fig. 5D), we normalized μ and σ to account for the different observed evidence e e
ranges at different positions. Let Ɛ be the set of all evidence values observed by the neuron at positions within 5 cm (one position bin) of μ . Define Ɛ+ = {e ∈ Ɛ | e > 0}, i.e. the set of positive p
observed evidence levels, and Ɛ- = {e ∈ Ɛ | e < 0}. For μ < 0, the normalized mean is given by e

μ̅̅̅ = μe ,
e e max(|e|)
e∈ℰ−

and for μ > 0, the normalized mean is given by e

𝜇̅̅ 𝑒 ̅ = ma 𝜇 x 𝑒 (e) .
e∈ℰ+

The normalized evidence standard deviation is given by

σ̅̅̅ = σe .
e e max(ℰ)−min(ℰ)

For visualization purposes (Fig. 5, Supplementary Fig. 5), we plot only the 80% of significantly evidence-selective neurons in each region with the best fit. For robustness to single large observations making a neuron appear narrowly tuned, for neurons with σ < 3, we identified the nearest position-e
evidence bin to (μ , μ ). If there were fewer than 4 observations in this bin or if there was a single outlier p e
in this bin (defined by one observation greater than 1.5 times the interquartile range), the neuron was excluded from plotting. For visualization of all significantly evidence-selective neurons without any further criteria applied, see Supplementary Fig. 12. In all density plots, kernel density estimates are generated using the python package seaborn, which uses scipy.stats.gaussian_kde to estimate the density with Gaussian kernels (bandwidth determined by Scott’s Rule).

Neural Evidence Tuning Curves Around the Most Active Position

To generate the 1-D evidence tuning curves around the most active position (Fig. 5B,F,J, Supplementary Fig. 5B) for each neuron with significant evidence tuning, we identified the subset of data within 0.5σ of μ . (Note that if this range extended to less than 0 cm or greater than 300 cm, only p p
the data within 0-300 cm were used.) Using the data from within this range, the tuning curves were computed by computing the average of the firing rates at each evidence level. We then fit the averages to two functional forms, a Gaussian and a logistic function, and plotted the curve that better fit the average tuning curve. By fitting to the averages, rather than the individual data points, we found we could better account for the uneven distribution of sampling across evidence levels at different positions (Supplementary Fig. 4A). The Gaussian evidence tuning curve fit took the form

−(e−μ)2
FR(e) = aexp( )+c
σ2

where FR(e) denotes the firing rate as a function of the evidence e. The fit was performed with a nonlinear fitting method (scipy curvefit, which employs the trust region reflective algorithm, with bounds m < μ < M, 0 < σ < 30, 0 < a < 10, 0 < c < 1, where m is the minimum and M is the maximum observed evidence within the active position). The logistic evidence tuning curve took the form

a
FR(e) = +c
(1+exp(−k(e−x )))
0

and was fit with a nonlinear fitting method (scipy curvefit, which employs the trust region reflective algorithm, with bounds -1 < k < 1, -15 < x < 15, 0 < a < 10, 0 < c < 1). To determine the curve with the 0
better fit, we calculated the mean squared error between the predicted firing rates from each fit curve and the true firing rates. We note that each form of the tuning curve has the same number of fit parameters, so there was no need to correct for the numbers of parameters when comparing model predictions.

For comparison to our fit tuning curves, we also plotted the raw evidence tuning curves around the most active position (Supplementary Fig. 6). We note that such raw evidence tuning curves around a given position require a reasonable assessment of the peak position tuning, as was obtained in our joint position-evidence fits described above. However, once the most active positions were obtained, the fits no longer depended on requiring the evidence to have the same functional form at every position or, in the case of the raw fits, on a particular functional form. Thus, we found that the fits could provide a more quantitatively precise visualization of the evidence tuning around the most active position and, for the neurons best fit by logistic functions, a better quantitative assessment of the steepness of the sigmoid. However, beyond this small quantitative difference, the fits typically did not differ substantially from the joint Gaussian fits, and required fitting to two separate functional forms, so we used the joint Gaussian fits for the population-wide assessment of tuning in Figure 5C,D,G,H,K,L.

Cue Encoding Model

To find the average response to a cue (Fig. 6A-E), we fit a linear cue encoding model to the difference between the average right-preferring activity and average left-preferring activity on each session, using the subset of neurons that are significantly evidence tuned:

f (p) = β +∑(K *ℓ)(p)+∑(K *r)(p)
pop 0 L R
Lℓ∈L r∈R

where f (p) is the difference between the average activity of the right evidence-preferring cells and the pop
average activity of the left evidence-preferring cells at position p, L is the set of left cues, R is the set of right cues, and K and K are the cue response kernels to left and right cues respectively. As above, L R
neurons that were significantly evidence tuned were categorized as right preferring if they had μ > 0 e
and left preferring if μ < 0.
e

For this analysis, we filtered the inferred spikes from the Ca++ imaging data using a causal half-Gaussian (window length = 1 s, standard deviation = 0.25 s), so that responses to cues would be purely causal. In order to fit the shapes of the cue response kernels of width 300 cm, each cue was convolved with a cubic spline basis set with 7 degrees of freedom148. We fit a ridge regression model with regularization strength 1 to predict the difference in population activity. The regression used the convolution of the cues with each of these splines as regressors, so that the cue encoding model was given by

Ns Ns
f (p) = β +∑β (s *ℓ)(p)+∑β (s *r)(p)
pop 0 k,L k k,R k
k k

where s is the kth basis function, N = 7 is the number of basis functions, and ℓ and r are binary vectors k s
with value 1 at the onset (defined as the time at which the cue becomes visible, 10 cm ahead of its position) of each left and right cue respectively.

Left cues and right cues were convolved separately to derive separate cue response kernels for each cue type. In Supplementary Figure 13A,B, separate cue kernels were fit for left and right cues when the magnitude of the current evidence prior to cue appearance, e, was low (|e| ≤ 1), medium (2 ≤ |e| ≤ 4), or high (|e| ≥ 5). For the medium and high evidence cues, the width of the cue kernel was 200 cm. In Supplementary Figure 13C,D, separate cue kernels were fit for left and right cues when the cue occurred in the early (before 70 cm), middle (between 70 and 140 cm), or late (after 140 cm) part of the cue period. For early cues, the width of the cue kernel was 300 cm, for middle cues, the width of the cue kernel was 230 cm, and for late cues, the width of the cue kernel was 160 cm.

Single-Neuron Evidence vs. Choice Encoding Analysis

For each neuron at each position bin p, we fit a single-neuron evidence vs. choice encoding model (Fig. 6F-H) to describe the firing rate of the ith neuron

FR (e,c,p) = β +β e(p)+β c
i 0 1 2

as a function of the current evidence e(p) at that position and the upcoming behavioral choice on that trial (c), using ordinary least squares regression. For this analysis, we fit to data from both correct and incorrect trials. For each coefficient, we calculated the F-statistic, and used the F-test (one-sided) to determine whether a coefficient was significant at a 0.01 significance level. At each position, for each session, we found the fraction of active neurons at that position that had a significant coefficient. For this analysis, a neuron was defined as active at position p if |p - μ | < σ , based on the results of the p p
joint position-evidence fitting procedure described above.

Average Evidence Population Responses for Trials with Different Evidence Levels

In Supplementary Figure 13E,F, we examined the average population responses across positions for correct trials with different final evidence levels. We binned trials by final evidence for left and right trials for low (1 ≤ e ≤ 3), medium (4 ≤ e ≤ 6), and high evidence levels (e ≥ 7). At each position bin, we only analyzed neurons that met two criteria: (1) they were active at this position bin, where active at position p is defined as above as having |p - μ | < σ , and (2) had a significant evidence p p
coefficient β (at a 0.01 significance level) at any position bin, based on the results of the single-neuron 1 L
encoding model described above. For these neurons, we calculated the difference between the average firing rates (extracted from the Ca++ imaging using the causal half Gaussian filter, as described above) of the right-preferring and left-preferring populations. As above, a neuron is defined as right preferring if μ > 0 and left preferring if μ < 0.
e e

Linear Evidence Decoding Model

To evaluate the ability to linearly decode evidence from the population on correct trials (Fig. 6I-K), we fit two decoders at each position bin, one for trials in which the evidence was greater than or equal to 0 and one for trials for which the evidence was less than or equal to zero. Separate decoders were calculated so that positive correlations between actual and decoded evidence do not result just from correctly decoding the animal’s choice on correct trials. To further eliminate correlations resulting from predicting just the sign-of-evidence, we also applied the same analysis but with one decoder for evidence strictly greater than zero and one for evidence strictly less than zero (Supplementary Fig. 13G,H).

For each decoder, we used nested five-fold cross-validation149 to fit a ridge regression model to predict evidence from the activity of the active neurons at that position. As above, neurons were defined as active at position p if |p - μ | < σ . In the inner loop of cross-validation, the regularization p p
strength was chosen from [.0001, .001, .01, .1, 1, 10, 100, 1000], with the best performing regularization strength used for evaluation on the test set in the outer loop. Decoders were evaluated by the average correlation between the predicted and actual evidence in the test set.

To compare to shuffled data, we repeated the same fitting procedure described above 5 times, but randomly shuffled the evidence values across trials of the same sign-of-evidence for each position.

Nonlinear Evidence Decoding in HPC

In Supplementary Figure 11, we use a nonlinear manifold method102 to decode evidence from the neural population across timepoints, unlike the linear evidence decoding models described above. On each fold of a five-fold cross validation, trials were assigned to training or test sets. To maximize the amount of data used to build the embedding, both correct and incorrect trials were used. We then fit a single CEBRA manifold102 (https://cebra.ai/) embedding of dimension d to the full trial trajectories in the training set for values of d from 2 to 7, using the following parameter settings: model_architecture = 'offset10-model', batch_size = 512, learning_rate = 3e-4, temperature_mode = 'auto', output_dimension = d, max_iterations = 10000, distance = 'cosine', conditional = 'time_delta', device = 'cuda_if_available', time_offsets = 10. At each point, we provided the current evidence value as the label for the neural data. We trained a k nearest neighbors decoder (k = 3) to predict the evidence level from the embedded training set. We then embedded the test set into the CEBRA manifold and used the decoder to predict evidence from the embedded test set. We then found the Pearson correlation between the predicted and actual evidence levels on the test set.

To compare performance to chance when controlling for choice, we repeated this process for 5 shuffles of the data in which the neural data was preserved but the evidence labels were permuted across trials within trials of the same behavioral choice.

Population-Averaged Tuning Curves

In Figure 6L,M, we calculated the population-averaged response across all significantly evidence-selective neurons as follows. For each neuron, we calculated the average firing rate in each position-evidence bin. If a bin had fewer than 3 observations for a particular neuron, that neuron was not included in the population average for that bin. To combine data for neurons that preferred positive evidence with data from neurons that preferred negative evidence, for each neuron, we defined its preferred-evidence as the sign of the neuron’s μ value and remapped the evidence axis to be in terms e
of preferred evidence rather than absolute evidence (by multiplying the evidence values by -1 if the preferred tuning was for leftward values of evidence instead of rightward). We then took the average across all neurons in each bin of position by preferred-evidence. Bins which were not sampled three or more times on at least 10% of the sessions for a particular brain region were not plotted.

Data Availability

The ACC, DMS, and HPC data have been uploaded to Figshare (https://doi.org/10.6084/m9.figshare.30921038)150. Previously published data from RSC is available at (https://datadryad.org/dataset/doi:10.5061/dryad.cvdncjt53)151. Source data are provided with this paper.

Code Availability

Code for the simulation of all models and the data analysis is available at https://github.com/lindseysbrown/evidence_accumulation_through_sequences152.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this manuscript, the authors developed two candidate circuit model classes to explain sequential (rather than persistent) firing patterns of neurons in an evidence-accumulation decision-making task. Such responses have been observed in many recent experiments. The models combined components from previously developed sequence models and evidence accumulation models to achieve sequential firing and evidence accumulation at the same time. The study highlights differences in neural tuning for evidence accumulation between the two models and describes how perturbation might additionally tell them apart. The study then compares tuning predictions to recorded calcium imaging data. Different forms of tuning are found in different areas, potentially suggesting that these areas might be using the different model classes. Overall, the structure and logic of the manuscript were both clear and the supplementary material contains useful workups of key features of the models, such as the stability of differences in amplitudes of the two chains across positions in the absence of external cue inputs. That being said, the models are not particularly novel, in the sense that they are first in essence a composition of previously known models in the literature, and second more of a specific instantiation of two already existing conceptual ideas for how evidence might be accumulated in sequences, competing chains (directly inspired from the data) and a progression in a space of evidence (as suggested by the co-authors in previous work and referred to in the discussion). It is interesting that different brain regions appear to have different tuning, but the non-monotonic tuning is mostly in hippocampus where previous work by one of the authors has shown place-cell like encoding for progression in non-spatial spaces (auditory) and thus may be less surprising. Nonetheless, having these models as concrete anchor points for these ideas is a useful contribution to the field.\nQuestions:\n1. The results from simulation and data analysis were always shown as an average across trials (and in places as averages across neurons as well). I understand that the task structure is more varied than some decision-making tasks in the way that visual stimuli are presented as Poisson sequences, but given that responses that are shifted in time average poorly I was surprised that there was almost no single-trial population level analysis. Both in general given previous literature, and since in classical attractor models, the nature of drift can be very informative of the underlying dynamics. Moreover, single trial analysis could perhaps show directly what are the changes in neural activity that happen when the same evidence is presented in two different positions, or the same position is traversed with or without evidence. In addition, it might give a better sense of how the models overall match the data.\n2. The fitting of firing rate as a function of position and evidence and the fitting of evidence tuning curves around the most active position were both done by solving a least squares problem, i.e., by minimizing the sum of squares of residuals. However, the metric used by the authors to evaluate the fit of the models was the correlation coefficient between the predicted firing rates and true firing rates. The two metrics are clearly related and likely highly correlated, but was there a reason for this switch between the optimized quantity during model fitting (minimized residuals) and the actual evaluation of models (correlation coefficient)? How will the results change (e.g. Fig 5) if residual was used as the evaluation metric of models?\n3. In the last results section “Progression in Graded Evidence Versus Binary Choice Tuning…” I may have misunderstood the analyses and my suggestion may have actually been performed, but I thought a helpful analysis would be comparing the classification (evidence-tuned vs choice-tuned) of a neuron with analysis in trials where the accumulated evidence is close to zero (i.e. the numbers of left and right stimuli are similar). It seems to me that neurons that are choice-tuned would still be able to predict the animal’s choice reliably in these trials, but neurons whose choice tuning is inherited through the association between evidence and choice would perform poorly. Was this the analysis described in the last paragraph? I think not but am unsure.\n4. Many example evidence-tuned neurons visually seem to be very choice selective, especially when the neurons are selective at the delay regions (e.g. Extended Data Fig 6, examples 27, 28, 29, 30, 31, 33, 34, 36). The large difference between negative and positive evidence (i.e., choice) are clearly visible, but it is hard to see graded firing rates under different evidence levels. Is this just a visualization issue or are they only weakly selective to evidence?\nMinor comments:\n1. Are the evidence tuning curves from data (Fig. 5B) informative of the difference between the coupled version of the competing chains model (Fig. 3D) and the uncoupled version (Extended Data Fig. 2D)?\n2. Would it be possible to directly fit the parameters of the competing chains model or the position-gated bump attractors using the experimental data?\n3. Reference 119: I believe the authors accidentally referenced a change in a manuscript and not the actual manuscript.\n4. In figure 1D, H, L I found the y-label “density” for what is described as a histogram in the legends a bit confusing, in particular as the scale changes five-fold from the left panel to the right panel.", "reviewer_2": "The current study presents a modelling effort to describe evidence accumulation in neural networks where neurons are activated sequentially as a rodent traverses an environment. Two types of models are considered (monotonic and ‘bump-shaped’ tuning). The authors analyze data from new and previously published datasets collected from mice performing the ‘towers task’ from several brain areas — in which such sequential activity is present — and study which model type corresponds more closely to the neural activity in each brain area.\nThe data analysis is nice, although the idea of joint evidence and spatial encoding is not novel. The modelling effort feels a bit raw.\nMajor Comments\n1. Weak link between model and neural activity. The authors construct models satisfying a series of desiderata for linear evidence accumulation in sequential models. However, the evidence that this phenomenon is taking place in the recordings is fairly limited.\nIn the areas showing ‘monotonic tuning’ (RSC and ACC) neurons seem to largely represent the best current categorical guess (Figs. 5 and 6). Although the authors argue that linear evidence accumulation transitions to binary choice encoding (lines 532-533) their Fig. 6D,F shows that, at best, the small fraction of neurons that have non-zero regression coefficients for either variable is similar. Then only binary choice becomes significantly represented later in the corridor. Furthermore, even at these early positions we are talking about ~ 2% of neurons, and also, at these early positions distinguishing between binary choice and graded evidence becomes more difficult (in the limiting case there only one tower has been seen they cannot be distinguished).\nIn HPC, convincing evidence is shown that a significant proportion of neurons have ‘bump shaped’ tuning curves for evidence, but this is also fairly indirect evidence that the population as a whole is tracking the running accumulated evidence in the task. There are certainly more direct ways to probe if the HPC is representing instantaneous accumulated evidence or not.\nThus, the bulk of the constraints that are used to develop the models derive from functional considerations that are not experimentally demonstrated. If, for instance, the data in Figs 1,5-6 was presented first, it would be more salient that the modelling in Figs 3,4 only loosely maps to the data. The data convincingly demonstrates is a dichotomy in terms of the shape of tuning curves to evidence, but this fact, without the requirements for graded evidence accumulation, provides a much weaker constraint on the mathematical models. If the authors wish to make a tight link between data and modelling, less constrained models (presumably a larger model class) should be considered. If the authors wish to explore possibilities for modelling sequential graded accumulation of evidence with a mostly theoretical motivation, then the theoretical and experimental aspects of the manuscript become less tightly associated. It would need to be made explicit that the data relates to these models fairly indirectly, and instead bears on a much larger model class, which seems to partly defeat the purpose of having the data and models together as it is done in the present version.\n2. Schematic and somewhat unrealistic modelling approach. While I appreciate the effort of translating the desiderata into specific mathematical implementations, the models presented seem somewhat artificial. A large number of tuning conditions need to be satisfied. The treatment of the all-or-none position signal is also quite artificial. The data in Fig1 shows graded and heterogeneous spatial ‘tuning curves’ and connected neurons will interact throughout the corridor, not just through the ‘position gate’.\nSome specific issues I had with the modelling:\n2.1 I might be mistaken, but it seems to me like the formula between lines 113 and 114 in the Supp. Text is incorrect? The upper limit of the integral is P_0, but it should be quantity with dimensions of time. Presumably P_0/speed, i.e., the time during which the position gate is open? Given the values of these parameters, a*P_0/speed = 1 Hz* 20 cm / 50 cm/s = 0.4, so that each position is actually integrated quite partially. As the authors point out, full integration would require this quantity to by >>1.\n2.2 Why does the neural activity profiles in figs 3C, 4C look so discontinuous along space? If the decay rate of activity is 1Hz and the mice runs at 50cm/sec (see comment below on speed), then graded decay should be evident on a 300 cm region.\n2.3 What is the role of speed on the activity of the network? If the position gate is fixed integration of evidence will proceed more slowly for faster mice? Is this reasonable?\n2.4 What is the shape of the dynamics of the sum mode in the competing chains model? Are there any qualitative features that could be compared with the data? The global activity of the network (and/or evidence-tuned neurons is readily accessible).\nOverall, readers will benefit from a more thorough characterisation of the models. More generally, while the models seem like a good first step, they seem too schematic and unrealistic for a finished product, specially given the fact that they are structurally unstable and the adding both functional heterogeneity (present in the data) and parameter heterogeneity (present in the brain) is not simply a matter of cosmetics, but rather represents an important challenge to the desired model function.\nOther comments:\n1. Fig 2 seems expected. Maybe it could go to the supplementary info.\n2. The optogenetic perturbation simulations in Figs 3 and 4, while interesting, take a lot of space and are not pursued experimentally. The authors might consider again moving them to the supplementary information and using that space to describe more precisely how the models work (in a sense deconstructing panel C by showing temporal profiles of synaptic currents from different sources etc), effect of different parameters, etc\n3. In lines 571-2 the authors speculate that the position signal (which in the model is totally external), might come place cells in the hippocampus, but the authors are actually modelling the hippocampus. This seems strange."}, "author_response": "We thank the reviewers for their feedback. Reviewer 1 praised the modeling as “a useful contribution to the field” and found that “the structure and logic of the manuscript were both clear,” the data analysis results were “interesting,” and the supplementary material was “useful.” Reviewer 2 also found that the “data analysis is nice” and “convincingly demonstrates … a dichotomy in terms of the shape of tuning curves to evidence.”\nIn addition, both reviewers provided important suggestions to improve both the modeling and analysis aspects of the paper. The two main criticisms of the work were, first, that our data analysis did not clearly show a representation of graded evidence, as assumed in the models and, second, that the models did not recapitulate a key feature of the data, in particular the presence of a binary choice population that increased in size as the animal traversed the maze.\nWe have addressed this constructive feedback by:\n1. Providing a more complete analysis of evidence representations in the brain regions that we analyze. In the broadest sense, we consider evidence tuning to be any predictable response in firing as a function of evidence, which can incorporate the (i) graded, monotonic evidence tuning of the chains-based models, (ii) the non-monotonic tuning of the position-gated bump attractor, and (iii) choice tuning, since choice will be correlated with evidence on correct trials. However, for a brain region to be a candidate integrator, it must have decodable representations of graded evidence. In the original version of the manuscript, we characterized the evidence tuning of individual neurons in each region and presented correlational analyses (single-neuron encoding models) that demonstrated that some neurons in ACC and RSC had graded encoding of evidence.\nWe now both improve upon and add to these analyses in the following ways:\n○ Revising our individual neuron evidence vs. choice encoding analysis models to explicitly separate behavioral choice from evidence. We have now updated our single-neuron evidence vs. choice encoding analysis to more directly account for behavioral choice by fitting our model to evidence and choice (rather than sign-of-evidence) on both correct and incorrect trials (rather than just correct trials). Incorporating incorrect trials allows us to better distinguish cells that respond to behavioral choice from those that respond to evidence, since the two are highly correlated on correct trials. This modification revealed a much larger portion of cells that responded to linear evidence beyond behavioral choice.\n○ Building evidence decoders to show that the regions we analyze carry graded evidence information. We now demonstrate that evidence can be linearly decoded, even after controlling for the current best choice, from the populations in ACC and RSC, lending further support that these regions have evidence representations. In addition, our decoding analysis shows that there is a transition from evidence coding to an increasingly choice-like representation later in the maze, where decoding of evidence falls to chance levels in the delay. This observed transition is consistent with recent work suggesting that neural dynamics change at the time of decision commitment (Luo et al. 2023).\nSeparately, we have added an analysis showing that evidence could be decoded from the HPC by utilizing a nonlinear manifold.\n○ Using a cue encoding model to uncover signatures of evidence integration. We fit a cue encoding model to the population activity that models neural activity as a linear sum of a temporal kernel associated with each cue (“cue kernel”), and found that cue response kernels in both ACC and RSC had the causal, sustained, approximately step-like responses to the appearance of cues. This signature is a necessary feature of an integrator.\n○ Incorporating single-trial analyses to further probe evidence responses in these regions. In addition to single-trial decoding analyses and population cue encoding model fits, we performed additional analyses of activity on individual trials. We showed that average trial trajectories showed graded patterns based on the final evidence during the cue region, which then collapsed to choice-like representations in the delay region, again consistent with a transition from evidence coding to an increasingly choice-like representation later in the maze.\n2. Making the connection between the data and the models more explicit. Our original models were built to integrate evidence, but our data analysis revealed that in addition to cells tuned to evidence, there is a large population of cells that respond to behavioral choice, especially toward the end of the maze. We now include two expanded circuit models that instantiate hypotheses for how such choice-tuning emerges as the animal traverses the maze, as well as a model for how position encoding may emerge from an integration of velocity.\n○ Explicitly modeling a population of choice cells that can read out the representations of both models. One of the major differences between the predictions of the chains-based models and the patterns observed in the neural data in regions with monotonic evidence tuning was the transition to a more choice-like representation later in the maze in the neural data. We now explicitly model two ways this could occur (Fig. 7): first, we model a population of choice-cells that respond to the difference in the two chains at different positions and show that, when there are increasingly more of these cells at later positions, we can reproduce the experimental, population-data results of Figure 6L,M. Alternatively, we show that the shift towards a more binary choice tuning can occur within the integrator if it is slightly unstable, as suggested for other accumulation of evidence tasks (e.g., Wong & Wang, 2006).\n○ Incorporating a mechanistic model of velocity integration into our bump-based model to allow position and evidence to be jointly integrated. In our previous version of the manuscript, we assumed that the position gating signal that controls which neurons are active derives from the observed spatially selective cells that do not encode evidence (these cells are now shown in Extended Data Figure 1). However, it is also possible that position and evidence are integrated within the same population. We now include in the manuscript an alternative, planar bump attractor model that integrates position along one axis and evidence along the other (Extended Data Fig. 3). This velocity-to-position integration mechanism could also be used to derive a position gating signal, separate from the evidence accumulation mechanism.\nIn addition to the major improvements outlined above, we have also made smaller updates based on reviewer feedback, including extending our discussion of the robustness of our models and demonstrating robustness to smoothness and heterogeneity in the position gate. The reviewers’ comments (black italics) and our individual responses (blue) to each point are included below.\nOverall, our work contributes two models that provide potential explanations for the choice-selective sequences that are present across regions during evidence accumulation and the different shapes of tuning curves in the hippocampus and in cortical regions. The additional analyses suggested by the reviewers helped to further characterize the differences across regions. We now show how these models can be modified to account for cells that respond to binary choice in the data. That being said, we note that one of the biggest experimental questions in the field is where the integration in evidence accumulation based decision making tasks is occurring, as it is difficult to distinguish the locus of integration from a readout of an upstream integrator – we hope that the predictions of our mechanistic models lay the groundwork for perturbation experiments that can further probe the coordination of evidence accumulation across the brain in cases where sequential, rather than persistent, neural activity is observed."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The reviewers responded in detail to my concerns and have performed a substantial amount of work in the revision. In my opinion, some of the responses, such as the addition of a comparison of correct and incorrect trials were useful and added to the paper. Others, like the clarification of regarding the use of \"evidence-tuned\" were important clarifications. However, my main concern, adding more single trial analyses, was less well answered. My comment was that there should be aspects of single-trial dynamics that are quite different between the two types of models. Just as the authors describe perturbation experiments that would better differentiate between the models, there should be single trial analyses that do the same. For instance, the authors write that (starting line 384): \"This is because the random walk of the bump location as evidence is accumulated leads to a greater range of evidence levels being observed at later positions in the trial; thus, on average across trials, there are more different neurons active at later positions, with each neuron corresponding to a different possible evidence level\". This an example property that could have been better probed with a single trial analysis. I don't believe this is the best example, I just mention it as an example of some single-trial analysis that makes expected model properties (sharper bumps) more apparent. My suggestion was for one or more analyses along such lines that clarify and then contrast properties of the two suggested models. The population analyses performed instead were more focused on population-level versions of overall properties. Analysis (i) shows that evidence can be decoded (which is not surprising given the single neuron results) analysis (ii) looks for sustained responses, which unless I am mistaken have to exist in both models and analysis (iii) is more of a visualization that better reveals graded encoding. These are all interesting, but do not address my specific comment which is on single trial dynamics that differentiate the two models. I appreciate the authors' response that drift is tricky and I fully appreciate that analysis and some others can't be done given the recordings available. My comment referenced that as one specific example of a broader set of analyses which I believe would have strengthened the paper.\nThat being said, I believe the paper is suitable for publication.\nTwo more minor comments: 1. I agree with reviewer 2 that Figure 2 is not particularly useful and can be sent to supplementary. 2. I personally find the model schematics very difficult to parse. To me the equations and verbal descriptions were more useful. While I imagine the authors have already dedicated some thought to this, it would be nice if they could find some better schematics (which can be for instance tested by asking people not directly involved with the study what they infer from the schematics).", "reviewer_2": "The authors have largely satisfactorily addressed my concerns.", "reviewer_3": "I have been brought in as an additional reviewer following missing responses from R2 in this review process. After reading the responses provided to R2's concerns, the authors have thoroughly and carefully answered all outstanding questions.\nR2's general critique was slightly vague at times, yet the authors have gone above and beyond, congratulations."}, "author_response": "We thank the reviewers for their positive feedback and recommendation that the manuscript is now “suitable for publication”. All reviewers found that we “responded in detail,” (Reviewer 1) “performed a substantial amount of work” (Reviewer 1), “largely satisfactorily addressed [their] concerns” (Reviewer 2), and “thoroughly and carefully answered all outstanding questions” (Reviewer 3), with Reviewer 3 even praising that our responses “have gone above and beyond.”\nThe reviewers’ comments (black italics) and our individual responses (blue) to each point are included below."}]

Decision letter