Distributed evidence accumulation across macaque large-scale neocortical networks during perceptual decision making

Despite the traditional view of parietal cortex as an important region for perceptual decision-making, recent evidence suggests that sensory accumulation occurs simultaneously across many cortical regions. We explored this hypothesis by integrating connectivity, cellular and receptor density datasets and building a large-scale macaque brain model able to integrate conflicting sensory signals and perform a decision-making task. Our results reveal sensory evidence accumulation supported by a distributed network of temporal, parietal and frontal regions, with flexible sequential decision pathways which depend on task difficulty. The model replicates experimental lesioning effects and reveals that the causal irrelevance of parietal areas like LIP for decision performance is explained by compensatory mechanisms within a distributed integration process. The model also reproduces observed temporal gating effects of distractor timing during and after the integration process. Overall, our work hints at perceptual integration during decision-making as a broad distributed phenomenon and provides multiple testable predictions.


Introduction
Decision-making is a fundamental cognitive function in the brain, involving the processing of sensory information towards a categorical decision and a subsequent motion response (Hanks and Summerfield, 2017).Electrophysiological recordings have characterized DM as an evidence accumulation process, where a decision is made when the firing rate of specific neuronal populations, typically in association parietal areas such as LIP, exceeds a threshold (Roitman and Shadlen, 2002;Shadlen and Newsome, 2001).This principle of sensory integration within local (parietal) circuits has been used to understand a vast range of decision-making paradigms, from simple perceptual decision-making tasks (Brunton et al., 2013;Roach et al., 2023;Wang, 2002;Wong and Wang, 2006) to context-dependent choices (Mante et al., 2013;Siegel et al., 2015;Song et al., 2016), multisensory integration (Chandrasekaran, 2017;Meijer et al., 2020;Raposo et al., 2014;Wierda et al., 2023) or even perceptual choices in the context of consciousness or arousal modulation (Beerendonk et al., 2023;Joglekar et al., 2018;Klatzmann et al., 2023;Mashour et al., 2020).
In recent years, however, compelling evidence has emerged against the local (parietal) sensory integration theory.Notably, lesioning and optogenetic studies in rodents have revealed that the inactivation of parietal cortex has little effect on the performance of decision tasks, even though the activity of those areas is highly correlated with decision speed and accuracy (Erlich et al., 2015;Lohuis et al., 2022).Similar experiments in macaques led to similar conclusions, suggesting that sensory areas like MT have a greater influence on achieving good performance than LIP (Katz et al., 2016).It seems implausible, therefore, that sensory integration necessary to guide decisions would occur solely on parietal circuits.Recent studies have explored whether the involvement of multiple brain regions could support perceptual decision-making (Pinto et al., 2022;Siegel et al., 2015;Steinmetz et al., 2019;Wilming et al., 2020), revealing a wide range of differentiated roles across cortical regions (Akrami et al., 2018;Hanks et al., 2015) and emergent interactions across different areas (Siegel et al., 2015).Promising conceptual advances have been facilitated by computational modeling work involving two or three interacting areas (Jaramillo et al., 2019;Murray et al., 2017;Wimmer et al., 2015).However, the lack of a data-constrained, brain-wide framework to compare and evaluate the contributions of individual areas involved in decision-making tasks has hindered further progress.
In this work, we integrated multiple brain-wide data sets and used them to develop a computational model of the large-scale network of the macaque neocortex with the capacity to simulate decision-making tasks.The data integrated into the model included cortex-wide variations in neuron density, autoradiography-derived NMDA and GABAA receptor densities, and average dendritic spine count per cell (Elston et al., 2011;Elston, 2007;Froudist-Walsh et al., 2023b), which contribute to capture part of the existing regional heterogeneity in cortex.By replicating a classical perceptual decision-making task (random dot motion discrimination), our model showed the involvement of multiple cortical areas in sensory accumulation, in contrast with the traditionally assumed local integration in parietal cortex.We also observed different roles, such as sensory, accumulators and categorizers, emerging across cortical areas.In addition, the model revealed the existence of two distinct dynamical regimes governing the flow of information during sensory integration, with the level of stimulus coherence determining which of them drove the decision process.By simulating lesions in different cortical regions, we reproduced the results of experimental inactivation studies, supporting the hypothesis of the causal irrelevance of area LIP for decision making and proposing a more important role for temporal lobe regions.Lastly, we propose that temporal gating during decision-making, a phenomenon observed in mice so far, may also be present in macaques due to the robustness of the distributed network supporting integration and maintenance of choice memory (Mejias and Wang, 2022).

A large-scale macaque cortex model incorporating data from multiple macroscopic gradients
We combined multiple neuroanatomical data sets with local dynamical models to build a largescale model of the macaque cortex (Figure 1A; see Materials and Methods for further details).We first considered a local circuit of two stimulus-selective excitatory populations and one nonselective inhibitory population.This local circuit was our template to describe each of the 40 cortical areas considered in our large-scale model, with properties of each individual circuit varying across areas following available area-specific data.This ensures that each cortical area was unique in its configuration and dynamics, and reflected the high level of neural and circuit heterogeneity found in real brains.The area-specific properties revealed macroscopic gradients across cortex (Wang, 2020), and included variations in the effective strength of synapses inferred from local properties such as the number of dendritic spines per pyramidal cell (Elston, 2007;Mejias and Wang, 2022) and the area-specific density of NMDA and GABAA receptors per neuron (Froudist-Walsh et al., 2023b, 2023a)(Figure 1B).
Before assembling these areas into a large-scale brain model, we analyzed the organization of the receptor density per neuron data when compared to the position of each area in the anatomical hierarchy (Markov et al., 2013;Nikola T. Markov et al., 2014) and area-specific dendritic spine numbers (Elston, 2007;Mejias and Wang, 2022).A strong positive correlation was found between NMDA receptor density and the cortical hierarchical position inferred from laminar connectivity data (r=0.731)(Figure 1C), and NMDA receptor density had an almost linear relationship with GABA receptor density (Suppl.Fig. S1).As previously shown (Chaudhuri et al., 2015;Mejias and Wang, 2022), dendritic spine count positively correlated with cortical hierarchical (Figure 1D).We identified here some parietal areas (such as LIP, 7m, and 7a) as outliers based on their low pyramidal spine count but high cortical hierarchy value and NMDA receptor density.This suggested that NMDA-mediated synapses in those areas preferentially target inhibitory neurons.
The data-constrained local circuits were then connected between them to form an extended brain network, using tract-tracing connectivity data of the macaque neocortex (Froudist-Walsh et al., 2021;Markov et al., 2013;Nikola T. Markov et al., 2014;N. T. Markov et al., 2014), leading to a brain model constrained at both the local and large-scale levels by neuroanatomical data.Using the layer-specificity and projection-directionality provided by this dataset, and following previous work (Feng et al., 2023;Mejias and Wang, 2022), we implemented a counterstream inhibitory bias in our network, which assumes that feedforward and feedback projections along the cortical hierarchy slightly but preferentially target excitatory and inhibitory neurons, respectively.Background input to cortical areas was varied across areas to mimic the differentiated thalamocortical projections across cortex and to facilitate that all cortical areas displayed a similar spontaneous activity level.

Distributed and heterogeneous evidence accumulation process in decision making
We simulated our large-scale macaque cortical model performing a classical two-choice randomdot motion discrimination task for decision making (Roitman and Shadlen, 2002;Shadlen and Newsome, 2001), where subjects have to choose which of two possible movement directions has a larger number of dots involved (Figure 2A).In our model, this visual input is modelled as two 700-ms-duration external currents entering both selective excitatory populations in V1 (Figure 2A and B), with dot movement coherence, or motion contrast, reflected in the difference between these currents (Figure 2B).This sensory input drives a competition between both selective excitatory V1 populations, and subsequently between the equivalent populations across other cortical areas.
Simulations of the large-scale brain model (Figure 2C) revealed that the competing visual stimuli induced the accumulation of sensory evidence and the resulting competition and winner-take-all dynamics in area LIP, as shown in experiments and classical models (Wang, 2002;Wong and Wang, 2006).However, such dynamics were also present in other brain areas, such as temporal lobe regions (TEO, TEpd), prefrontal cortex (9/46v), and many others (Suppl.Fig. S2).This suggests that the decision process, including sensory accumulation and winner-take-all dynamics, occurs across a distributed cortical network rather than in localized parietal regions.The dynamics presented important differences across the cortex: areas like LIP, TEO and 7M showed a slow and gradual ramping process, frontal areas like F5 and 9/46v underwent a sharp transition reflecting a categorical choice dynamic, and temporal areas like MT and TEpd displayed mutual ramping after stimulus onset and then bifurcated to encode the final decision.As in classical models and experimental evidence on single areas, the accumulation process was dependent on the level of stimulus coherence (Figure 2D): as illustrated by this example in LIP, ramping activities bifurcated earlier in high-coherence trials compared to low-coherence ones, and experienced distinct trajectories for two situations.In both cases, the winning population displayed sustained activity after the integration process, due to the bistability of the system, and related to working memory capacity of the network (Mejias and Wang, 2022).When adopting a majority rule to infer the behavioral decision from the winning populations across all cortical areas, the large-scale brain model was also able to replicate psychometric curves in agreement with macaque behavioral output data (Figure 2E), as well as plausible chronometric curves (Figure 2F).The reaction time in chronometric curves was defined as the time at which populations displaying gradual ramping activity (like LIP) reached a threshold of 15 Hz, signaling a joint contribution of association areas towards a decision.(Roitman and Shadlen, 2002).
Winning onset times uncover motion information flow across the cortex While evidence accumulation was observed across many cortical areas, integration times greatly varied across regions.We analyzed the temporal integration across multiple areas by comparing the normalized firing rate traces of the populations encoding the winning choice.For low coherence levels, temporal and prefrontal areas (like TEO and 9/46v, respectively) led the integration and were the first ones to reach their half-maximum level (Fig. 3A, left).However, the situation changed for trials with high coherence levels, as visual areas like V4 and MT were the first ones reaching their half-maximum level (Fig. 3A, right) and temporal and frontal areas evolved at a slower pace.This suggests that the information flow during decision-making is stimulus-dependent, and that the functional relationships between cortical areas (i.e. which area leads the integration), rather than being a fixed structure, flexibly depends on task difficulty.
To estimate the moment in which a winner emerges in the local competition between selective populations in each area, we defined the 'winning onset' in a cortical area as the moment in which the difference between their firing rates becomes larger than a certain threshold (zero, for simplicity) and continues evolving towards the decision threshold without returning back to zero (black dot in Fig. 3B).Intuitively, the winning onset signals the time at which a winner population gains an edge over the other one, and a differentiated integration emerges.The winning onset provides a simple way to evaluate which areas lead the integration process, which aligns reasonably well with decision sequences observed in categorical decision tasks in macaques (Siegel et al., 2015)(see Suppl.Fig. S3 for a direct comparison).We simulated, for a given coherence level, 2000 decision-making trials and computed the winning onset for each area.By sorting these winning onset times and ranking all areas from first/fastest to last/slower, we obtained a trial-averaged distribution of winning-onset ranks for each area (Fig. 3C).
For low coherence values, frontal areas like 9/46v and 45A, as well as temporal areas like TEO, STPc and STPi were usually ranked in the first places.Areas linking sensory and association areas, like MT and LIP, as well as vision-unrelated areas such as areas 2 and 5 (primary and secondary somatosensory cortex), fell behind in the ranking.This suggests that, under low coherence input, higher association areas led the sensory integration and decision process, which was later reflected in the activity of areas like LIP and MT.For high coherence values, the pathway for sensory integration drastically changed (Fig. 3C, right).In this condition, LIP, MT and similar areas connecting sensory and association cortex become fast integrators and reclaim the first places in the winning onset ranking.Association areas, on the other hand, became slower and decreased in the ranking (although some of them, like TEO, remain in moderately high ranks).Overall, high-ranked areas seemed to be hierarchically close to early sensory areas V1 and V2, indicating that sensory integration follows the anatomical hierarchy for high coherence input.
Two distinct mechanisms accounted for decision making in our model, depending on the coherence level.We focused on the rank order distribution of area MT to further illustrate the intrinsic mechanisms.Area MT experienced a remarkable transition from low to high winningonset ranks when coherence increased from 0% to 100% (Fig. 3D).For low coherence, the differentiated sensory drive received by MT was weak, and due to the low self-coupling strength of this area, it took some time for its firing rate to reach the winning onset.Association areas, which typically have stronger self-coupling, reached their winning onset points faster in this condition, and hence ranking higher.For higher coherence levels, the differentiated sensory drive was strong enough to trigger a winning onset (and therefore a choice) in MT and related areas, and this information was then passed to higher association areas which followed MT.
The rank switch phenomenon could be further illustrated with a toy model (Fig. 3E) including only V1, MT and 9/46v.Here we rescaled the anatomical connections while keeping their relative ratios, and kept other parameters constant (see Methods).As in Fig. 2C, we observed activity ramping in both MT and 9/46v (Suppl.Fig. S4).Winning onset times decreased with coherence level for both areas (as easier trials lead to faster decisions), and the intersection between both curves signaled a switch in the winning onset ranking.Inactivation reveals the functional irrelevance of LIP during decision making Inactivation of specific brain regions, either by lesioning, pharmacology or optogenetics, has provided valuable insight for decision-making circuits.In macaques, inactivation via the GABAA receptor agonist muscimol revealed a profound impact of MT inactivation in behavioral performance, while LIP inactivation had no measurable impact on performance (Katz et al., 2016).We used our model to simulate inactivation of MT and LIP areas (Fig. 4A), and we compared our modeling results with experimental evidence.We found a strong agreement between our results and the experimental findings of Katz and colleagues (Katz et al., 2016), with LIP inactivation barely affecting the psychometric curve for our model and MT inactivation leading to a flattening of the sigmoidal, reducing the accuracy of the decision process (Fig. 4B).Our results therefore support the evidence of LIP playing a relatively minor role in decision making, despite accurately reflecting the sensory accumulation process occurring across other distributed cortical networks.The effects of inactivation were significantly smaller for higher values of the global synaptic coupling strength (G=0.52), as in this case the inactivation of any area could be easily covered by the joint contribution of other areas (Suppl.Fig. S5).
To better explain this result, we evaluated how the percentage of valid trials, defined as those for which the model was able to make a decision, depended on the coherence level for different global coupling values.For small global coupling, lesioning MT prevents the sensory evidence to reach higher cortical areas, leading to noise-driven trials and poor performance.Larger global coupling values ensure instead that trials are evidence-driven, leading to a good performance for MT inactivation due to compensatory effects of other regions.The inactivation of LIP, however, can be compensated by sensory accumulation in other areas for both small and large global coupling, leading to a causal irrelevance of LIP for the decision process (Fig. 4C).
Our distributed computational model was also able to infer how lesioning one area would affect the dynamics of other areas, especially their ramping speed which is a key feature of the evidence accumulation process.To precisely define ramping speed, we introduced three areaspecific metrics besides the winning onset: (i) the winning rate, defined as the firing rate of the winning population at the winning onset, (ii) the reaching rate, defined as 75% of the attractor firing rate minus 25% of the winning rate, and (iii) the reaching onset, defined as the time at which the winning population firing rate arrives at the reaching rate (Figure 4D).Ramping speed was defined as the change in firing rate over the ramping time window, and the lesioning effect of area  on area  was the percentage change between pre-and post-inactivation ramping speed.To better capture the ramping dynamics, we use a winning onset threshold of 0.5 instead of zero.
We then observed that lesioning MT had a stronger effect on other non-sensory areas than lesioning LIP (Figure 4E).Though lesioning LIP has a stronger impact in certain areas like 7B or STPi, the ramping speed of most areas (especially LIP, 7m, 8l and 10) changed more drastically when lesioning MT.Notably, lesioning LIP led to ramping speed increases in areas like F3 and TEpd; the reason was that long-range projections from LIP on those areas targeted inhibitory populations more strongly than excitatory ones, so removing LIP had a positive effect on those areas.This effect was not observed when lesioning MT, which follows from its early position on the feedforward cortical hierarchy and its role in the excitatory feedforward pathway.

Effects of lesioning areas across the whole cortex
We extended the study of the previous section by lesioning all non-sensory areas one by one (Fig. 5A) and analyzed how any particular inactivation affects ramping speeds across cortex.The full inactivation matrix is displayed in Fig. 5B.We observed, for example, that lesioning prefrontal area 46d had a strong effect on the ramping speed of parietal areas LIP and 7m, indicating a supportive role of PFC on parietal cortex.Furthermore, lesioning TEO greatly decreased the ramping speed of TEpd and increases that of area 45A and the OPro.Lesioning temporal areas like STPc, STPi and MT affected virtually all association areas.
When we took mean values of the lesioning effect matrix, effectively applying the majority rule, we could infer how lesioning specific cortical areas globally affected the macaque's decisionmaking performance (Fig. 5C) and the mean ramping speed across areas (Fig. 5D).For example, frontal areas like OPro, PRoM, 32 and 25 had little effect on the psychometric curve (as indicated by changes in its slope), since their interareal connectivity with other regions was relatively weak.This link between structural connectivity and lesioning effects was reflected in a negative correlation of  = −0.39 between both factors (Fig. 5E).In contrast, inactivation of temporal areas led to strong effects, suggesting a more salient role of the temporal lobe in decision making.

Temporal gating in the macaque brain
In our large-scale brain model, information about the decision taken is stored in working memory as distributed sustained activity patterns, as in previous work (Mejias and Wang, 2022) and experimental evidence (Christophel et al., 2017;Leavitt et al., 2017;Voitov and Mrsic-Flogel, 2022).For categorical decision-making tasks, the robustness of such memory representations is important for translating decision outcomes into motor output.This has been experimentally tested in rodents and macaques, by presenting distractors at different time points across the delay period, revealing that the robustness of decision outcomes greatly depends on the distractor timing, a phenomenon known as temporal gating (Finkelstein et al., 2021;Seidemann et al., 1998;Suzuki and Gottlieb, 2013).We tested the robustness of our decisionrelated working memory in a decision-making task with distractors (Fig. 6A).Two versions of the task were considered: one in which distractors have a fixed onset but variable duration, and one in which they have a fixed duration but variable onset: at the sampling, early delay, or late delay phases (Fig. 6B).
We first studied the system in the fixed duration version of the task, and calculated the fraction of trials in which the decision was changed as a consequence of a distractor at different onsets (Fig. 6C).We observe that early distractors, especially those presented during the sampling phase, have stronger effects on changing the decision output than late ones, in agreement with experimental observations (Finkelstein et al., 2021;Seidemann et al., 1998).The effects were similar for weak and strong distractors, with the latter simply scaling up the fraction of trials with reversed decisions.We observed a positive relationship between the distractor onset and the minimal strength needed for the distractor to be effective, i.e. to change the decision outcome (Fig. 6D, left panel).For the fixed onset version of the task, we found an inverse relationship between the distractor strength and the minimal duration required for the distractor to be effective (Fig. 6D, right panel).
The results were similar for networks with a stronger global coupling (Fig. 6E), although with more abrupt limits in distractor effectiveness.For example, in the fixed duration version of the task (left panel), the decision was nearly impossible to be reversed with distractors of very late onset, irrespectively of the distractor strength.Likewise, decisions could not be reversed with very weak distractors (right panel), irrespectively of the distractor duration.Changes in distractor duration (for the fixed onset version) or distractor strength (for the fixed duration version) were determinant factors for the success of the distractor in overturning the decision (Suppl.Fig. S6).

Discussion
In this work we investigated the distributed computations underlying decision-making tasks, by simulation of large-scale cerebral network models.Traditionally, decision making has been characterized by a gradual evidence accumulation process reflected in ramp-up activity in local cortical circuits (Roxin and Ledberg, 2008;Wang, 2002;Wong and Wang, 2006;Ye and Li, 2021).While a few computational studies have approached this topic considering two interconnected areas (Jaramillo et al., 2019;Murray et al., 2017;Wimmer et al., 2015), our work provides a formal leap as a novel and detailed computational study of decision making encompassing a large-scale, data-constrained cortical network.Additionally, our study opens the door to explore cognitive functions within large cortical networks in an integrative manner, as its architecture and properties align with a recent proposal regarding the distributed nature of working memory (Mejias and Wang, 2022).
A salient feature of our work is the integration of multiple data sets of gradients of biophysical properties across the cortical network: (i) a pyramidal cell dendritic spine gradient informing on local synaptic density, (ii) a tract-tracing connectome embedding the model with directed and weighted long-range connections, (iii) a neuroanatomical cortical hierarchy used to infer targeted populations in long-range projections, and (iv) gradients of receptor density per neuron for the NMDA and GABAA receptors which allowed to constrain the excitatory/inhibitory balance and strength.This level of detail allowed us to identify outliers like LIP or 7A, which had a high hierarchical rank and NMDA receptor density per neuron but a low level of pyramidal dendritic spines.This pointed at an unexpected abundance of NMDA receptors on inhibitory neurons, which contributed to the observed slow ramping activity in parietal cortex compared to other areas such as prefrontal cortex.The combination of data sets also allowed our model to explain sustained activity as reported in areas V4 and MT (Bisley et al., 2004;Hayden and Gallant, 2013), which previous models had missed.Combining these data sets with additional refinements, such as causality-reconstructed connectivity maps and hierarchical interactions (Chaudhuri et al., 2015;Li and Wang, 2022;Mejias et al., 2016;Nunes et al., 2021;Zhou et al., 2014) or additional receptor per neuron densities (Froudist-Walsh et al., 2023b, 2021;Hansen et al., 2022;Zilles and Palomero-Gallagher, 2017) may further increase the model's predictive capacity.

Winning onsets uncover multiple information pathways
The sequential order of winning onsets observed in our model was a reflection of how motion information may be processed and broadcasted in the brain.Interestingly, the information pathways revealed by the model were not fixed, and rather dependent on the input coherence and therefore the task difficulty.We observed that higher cognitive and association areas like TEO and 9/46d led the decision process for hard trials, while early sensory areas led for easier trials.The alteration in winning onset order could indicate a change in the role of certain areas when the macaque encountered a tough vs simple task.This suggests that an efficient performance on difficult decision-making tasks might benefit from top-down modulation, a prediction aligned with conscious perception theories (Joglekar et al., 2018;Klatzmann et al., 2023;Mashour et al., 2020;Vugt et al., 2018).There is also evidence of monkeys taking different strategies depending on task difficulty (Brunet and Jagadeesh, 2019;Hyafil and Moreno-Bote, 2017).
In addition, previous experimental attempts to uncover the information flow during categorical perceptual decision tasks (Siegel et al., 2015) revealed a decision-related cortical pathway given by MTà LIPà V4à ITà FEFà IPFC.By setting the input to full coherence, therefore approximating a categorical perceptual decision task, our model led to a similar cortical pathway for information flow attending to winning onset order: V4à MTà IT(TEO)à LIPà FEF (8I)à IPFC (9/46v) (Suppl.Fig. S3), with only minor misalignments with the data for the positions of V4 and TEO/IT.Furthermore, the activation of this sequence occurs within 0.2 seconds for both experimental data and our model predictions.This indicated that dataconstrained large-scale models are well-suited to explore decision-related information pathways in cortical networks, and to uncover the mechanisms behind sequential decision stages (Finkelstein et al., 2021;Guo et al., 2014).

Causal irrelevance of LIP is replicated by simulated inactivation
Despite the traditional focus on parietal cortex as a circuit involved in sensory accumulation, recent evidence on rodents and primates has disputed this idea (Akrami et al., 2018;Erlich et al., 2015;Hanks et al., 2015;Katz et al., 2016;Pinto et al., 2022).However, a mechanistic explanation as of why parietal cortex would play a relatively minor role in decision making despite its activity's high correlation with behavior was crucially lacking.By simulating lesioning effects across different areas, our model demonstrates that a distributed network of redundant evidence accumulation is sufficient to explain the experimentally observed causal irrelevance of LIP, as its inactivation may be compensated by the integration process occurring in parallel in other areas.Inactivation of other areas such as MT led by contrast to a much more substantial impairment in performance, as also observed experimentally.This is mostly due to the relatively low hierarchical position of MT, acting as an important bridge in the sensory pathway.When further exploring currently unknown effects of lesioning other areas, the model predicted an important role of temporal lobe areas such as TEO and TEpd, a prediction that requires experimental validation.

Temporal gating in decision-related memory
In decision making tasks, monkeys are sometimes expected to maintain their choice in memory for a short time, which posits the robustness of choice-related memories as fundamental for performance.We varied global coupling strength to individually simulate both strongly and weakly distributed circuit models.As in (Finkelstein et al., 2021), distractors with different duration and strength were fed to the model, revealing distinct reactions to distractors and suggesting that early distractors are most effective than late ones.In weakly distributed models, it was easier for distractors to disturb the ramping activity, perhaps reflecting the performance of poorly trained monkeys.Future experiments could address how monkeys at different trained stages react to distractors, to test the existence of these two effects behind temporal gating (Treue and Maunsell, 1999).

Materials and Methods
The large-scale brain model presented here uses previous modeling work on distributed working memory as a basis (Mejias and Wang, 2022), and substantially improves upon it with the introduction of (i) a larger number of cortical areas interconnected by neuroanatomical data, (ii) area-specific data on the density of NMDA and GABA receptors across the brain, obtained via autoradiography experiments, and (iii) novel theoretical assumptions on the relative impact on the effects of the density of dendritic spines and postsynaptic receptors on intra-and inter-areal interactions.
We describe our model in three steps: the local circuit employed, the differentiation of such circuit across the whole neocortex via synaptic gradients, and the inter-areal projections.We then briefly define several metrics for behavioral performance and inactivation effects, ending with a short summary on our simplified 'toy model' for winning onset results.

Computational model: local neural circuit
We employed the Wong-Wang model (Wong and Wang, 2006) to characterize the neural dynamics of a local microcircuit representing a cortical area.This model, in its three-variable version, captures the temporal evolution of the firing rates of two input-selective excitatory populations, as well as the firing rate dynamics of an inhibitory population.The populations are interconnected with each other, as depicted in Figure 1A.The model is governed by the following equations: Here, SA and SB are the NMDA conductances of selective excitatory populations A and B respectively, and SC is the GABAergic conductance of the inhibitory population.Values for the constants are τN=60 ms, τG=5 ms, γ=1.282 and γI=2.The variables rA, rB and rC are the mean firing rates of the two excitatory and one inhibitory population, respectively.
They are obtained by solving, at each time step, the transcendental equation  ) =  ) ( ) ) (where  is the transfer function of the population, detailed below), with Ii being the input to population 'i', given by ( =  ',  % +  ',  & +  ''  ( +  -( +  ./# (+  ( () (Eq.6) In these expressions, Js, Jc are the self-and cross-coupling between excitatory populations, respectively, JEI is the coupling from the inhibitory populations to any of the excitatory ones, JIE is the coupling from any of the excitatory populations to the inhibitory one, and JII is the selfcoupling strength of the inhibitory population.The parameters I0i with i=A, B, C are background inputs to each population.Fixed parameters across the cortex are Jc=0.0107nA, JEI=-0.31nA, JII=-0.20 nA, and I0C=0.26nA (background currents for populations A and B are specified later).The term I i net denotes the long-range input coming from other areas in the network and will be detailed later.Sensory stimulation can be introduced here as extra pulse currents of strength IstimA =0.3(1+c) and IstimB=0.3(1-c)to the V1 excitatory populations A and B respectively, with a duration of Tpulse =0.7 sec, unless specified otherwise.
The last term xi(t) with i=A, B, C is an Ornstein-Uhlenbeck process, which introduces some level of stochasticity in the system.It is given by Here, ξi(t) is a Gaussian white noise, the time constant is τnoise=2 ms and the noise strength is σA,B=0.01nA for excitatory populations and σC=0 for the inhibitory one.
The transfer function ϕi(t) which transform the input into firing rates takes the following form for the excitatory populations (Abbott and Chance, 2005): (Eq.8) The values for the parameters are a=135 Hz/nA, b=54 Hz and d=0.308 s.For the inhibitory population a similar function can be used, but for convenience we choose a threshold-linear function: The notation [] A denotes rectification.The values for the parameters are gI=4, cb=615 Hz/nA, ca=177 Hz and r0=5.5 Hz.Finally, it is sometimes useful for simulations (although not a requirement) to replace the transcendental equation  ) =  ) ( ) ) by its analogous differential equation, of the form: The time constant can take a typical value of τr=2 ms.

Computational model: Gradient of synaptic strengths
Before considering the large-scale network and the inter-areal connections, we look into the areato-area heterogeneity to be included in the model.
Our large-scale cortical system consists of N=40 local cortical areas, for which inter-areal connectivity data is available.Each cortical area is described as a Wong-Wang model of three populations like the ones described in the previous section.Instead of assuming areas to be identical to each other, here we will consider some of the natural area-to-area heterogeneity that has been found in anatomical studies.For example, work from Elston (Elston, 2007) has identified a gradient of pyramidal cell dendritic spine density, from low spine numbers (~600) found in early sensory areas to large spine counts (~9000) found in higher cognitive areas.From an electrophysiological point of view, excitatory postsynaptic potentials (EPSP) have similar values both in early sensory (~1.7+1.3 mV) and higher cognitive areas (~0.55+0.43mV).The combination of these findings suggests an increase of local recurrent strength as we move from sensory to association areas.In addition, cortical areas are distributed along an anatomical hierarchy (Felleman and Van Essen, 1991;Nikola T. Markov et al., 2014).The position of a given area 'i' within this hierarchy, namely hi, can be computed with a generalized linear model using data on the fraction of supragranular layer neurons (SLN) projecting to and from that area.
In particular, we assigned hierarchical values to each area such that the difference in values predicts the SLN of a projection.Concretely, we assign a value Hi to each area Ai so that SLN(Aj à Ai) ~ f (Hi-Hj), with 'f' being a logistic regression.The final hierarchical values are then obtained by normalizing hi=Hi/Hmax.Further details on the regression are provided elsewhere (Chaudhuri et al., 2015;Nikola T. Markov et al., 2014).
In the following, we will assign the incoming synaptic strength (both local and long-range) of a given area as a linear function of the pyramidal dendritic spine count values observed in anatomical studies, with age-related corrections when necessary.Alternatively, when spine count data is not available for a given area, we will use its position in the anatomical hierarchy, which displays a high correlation with the spine count data, as a proxy for the latter.After this process, the large-scale network will display a gradient of local and long-range recurrent strengths, with sensory/association areas showing weak/strong local connectivity, respectively.We denote the local and long-range strength value of a given area i in this gradient as hi, and this value normalized between zero (bottom of the gradient, area V1) and one.We assume therefore a linear gradient of values of Js, with its value going from Jmin to Jmax: The above equations consider the impact of the number of dendritic spines per neuron in the synaptic strength of each brain area, but the density of NMDA and GABAA receptors per neuron must also be taken into account.Since excitatory signals target not only other pyramidal neurons but also interneurons, we assume that the local NMDA receptor density has a positive correlation with the sum of all excitation signals in a given brain region.We use a saturating gradient, or logistic function, to model such area-specific parameters: 3A7 ) 89: (6D(EFG% ' 6! ) )) +  3 (Eq.12) Here, the sum of the three excitatory projections in the left-hand side define the strength of excitatory transmission, which depends on  ) , the experimental NMDA receptor density of brain area 'i', via a logistic function.The constants  3 = 0.5241,  3 = 1.01,  3 = 0.01846,  3 = 0.12,  = 9.35 are used constrain the slope and roof of the curve.Since receptor density ratio between NMDA and GABA measures the relative excitatory-inhibitory synaptic strength balance, GABAergic receptor density can be introduced via the following balance equation: H%&% " ' EFG% ' (Eq. 13) where  ) is the experimental GABAA receptor density of brain area 'i'.
If association areas have large values of  " , it can influence their spontaneous activity, even without considering inter-areal coupling.To ensure that the spontaneous firing rate of these areas remains within physiologically realistic limits, a viable approach is to enforce that the fixed point of spontaneous activity is the same for all areas, which is a reasonable approximation.This can be done by adjusting the background currents  -(for both excitatory populations A and B) on an area-specific basis, which aligns with the differentiated thalamocortical input observed in real brains.We therefore use adaptive background current to balance spontaneous firing rate: ) −  ,' K'L T (Eq.16) Here  - K'L is the background current in parietal area LIP and its value has been set  - K'L = 0.3294 as in previous work (Wong and Wang, 2006) to fit temporal ramping excitatory neuron firing rates in LIP.The parameter  *M = 0.03566 is the constant value of spontaneous NMDA gating variable in a local disconnected circuit.
Unless specified otherwise, we assume a range of  C). =0.225  and  C51 = 0.42  (i.e., below the critical value), so that the model displays distributed attractors as in the case of working memory (Mejias and Wang, 2022) Computational model: Inter-areal projections In the present study, the inter-areal projections that connect isolated areas contribute to the formation of the expansive cortical network.We assume that these inter-areal projections originate solely from excitatory neurons, as inhibitory projections typically exhibit a more localized pattern within real circuits and are selective for excitatory neurons.Thus, the network or long-range input term within a specific area "x" from all other cortical areas, can be expressed as follows: Here,  is the global coupling strength, which controls the overall long-range projection strength in the network (G=0.52 unless specified otherwise). = 1.2 is a factor that considers the relative balance between long-range excitatory and inhibitory projections.
Aside from global scaling factors, the effect of long-range projections from population y to population x is influenced by two factors.The first one,  1N , is the anatomical projection strength as revealed by tract-tracing data (Markov et al., 2013).We use the fraction of labelled neurons (FLN) from population y to x to constrain our projections values to anatomical data.We rescale these strengths to translate the broad range of FLN values (over five orders of magnitude) to a range more suitable for our firing rate models.We use a rescaling that maintains the proportions between projection strengths, and therefore the anatomical information, that reads Here, the values of the rescaling are k1 =1.2 and k2 =0.et al., 2011;Roxin, 2011), this was done to have a better control of the heterogeneity levels of each area, and to minimize confounding factors such as the uncertainty on volume injections of tract tracing experiments and the influence of potential homeostatic mechanisms.In addition, and as done for the local connections, we introduce a gradient of long-range projection strengths using the spine count data:  1N → ( * ()/ PQ9 )  1N , so that long-range projections display the same gradient as the local connectivity presented above.
The second factor that needs to be taken into account is the directionality of signal propagation across the hierarchy.Feedforward (FF) projections that are preferentially excitatory constitute a reasonable assumption which facilitate signal transmission from sensory to higher areas.On the other hand, having feedback (FB) projections with a preferential inhibitory nature contributes to the emergence of realistic distributed WM patterns (Figure 4) (see also previous work (Nikola T. Markov et al., 2014;Tsushima et al., 2006)).This feature can be introduced, in a gradual manner, by linking the different inter-areal projections with the SLN data, which provides a proxy for the FF/FB nature of a projection (SLN=1 means purely FF, and SLN=0 means purely FB).In the model, we assume a linear dependence with SLN for projections to excitatory populations and with (1-SLN) for projections to inhibitory populations, as shown above.
Following recent evidence of frontal networks having primarily strong excitatory loops (Markowitz et al., 2015), it is convenient to ensure that the SLN-driven modulation of FB projections between frontal areas is not too large, so that interactions between these areas are never strongly inhibitory.In practice, such constraint is only necessary for projections from frontal areas to 8l and 8m (which are part of the frontal eye fields) and has little effect on the behavior of our model otherwise.The introduction of this limitation has two minor consequences: (i) it allows area 8l and 8m to exhibit a higher level of persistent activity during distributed WM -as their hierarchical position and recurrent strength are not strong enough to sustain activity otherwise, as previously suggested in anatomical studies (Markov et al., 2013;Nikola T. Markov et al., 2014), and (ii) it slightly shifts the transition point in cortical space.
Unless specified otherwise, we consider that the SLN-driven modulation of FB projections to 8l and 8m is never larger than 0.4.

Behavioral performance: psychometric curves
Psychometric curves were fitted with a binomial generalized linear regression model (GLM): where  was the probability to choose population A as the final decision and  was the contrast strength. = ( 3 ,  -) were the parameters of the fit, with  3 representing the slope of the psychometric curve and  -the value at zero contrast.

Evaluation of lesioning effects
To replicate the results of the lesioning protocol, we employed a manual approach to induce lesioning effects in the target area's three populations.Specifically, we set the firing rates of all three populations to zero.We evaluated the impact of such lesioning by assessing the percentage change in ramping speed, which measures the speed at which a specific area encodes the final decision: Here, RSi is the ramping speed for area  in the decision-making task with the intact recurrent network.As mentioned in the main text, the winning onset (WOi) is the time at which the difference between the firing rates of both excitatory populations becomes larger than a certain threshold (zero for Fig. 3, 0.5 for Figs. 4 and 5) and increases from there without returning to zero (see Fig. 3B).The winning rate (WRi) is defined as the firing rate of the winning population at the winning onset, the reaching rate (RRi) is defined as (0.75 x attractor firing rate -0.25 x winning rate), and the reaching onset (ROi) is defined as the time at which the winning population firing rate arrives at the reaching rate.We compare how lesioning one area affects the other area under same realization of Gaussian white noise with the following expression: where LEi,j is the degree by which lesioning area 'j' affects the ramping speed of area 'i', and RSi,j is the ramping speed in area 'i' when area 'j' is lesioned.
A toy model accounting for the shift of winning onset order The toy model of Figure 3E consists of three cortical areas V1, MT and 9/46v, with V1 transmitting motion information, MT as an intermediate area connecting association areas with sensory areas, and 9/46v as an association area to sustain delay activity.The computational model structure is the same except for the long-range projections, which are now as follows: Here,  * 1N ,  6 1N are parameters assigning individual roles to the three areas.All parameter values of the model can be found in Supplementary Table S1.

Figure 1 :
Figure 1: Scheme and anatomical basis of the large-scale macaque brain model.(A) Lateral view of the macaque cortical surface with modelled areas in color.In each area, cortical dynamics follow a local winner-take-all model, and areas are connected using anatomical connectivity data.(B) Area-specific values for the number of dendritic spines per neuron (top), the densities of NMDA receptors per neuron (middle) and GABAA receptors per neuron (bottom).(C) Correlation between NMDA receptor density and anatomical hierarchy as defined by layerdependent connections.(D) Correlation between spine count data and anatomical hierarchy as defined by layerdependent connections.

Figure 2 :
Figure 2: Large-scale brain model reveals a distributed sensory accumulation process for decision making.(A) Random-dot motion discrimination task simulated in our work (left panel), where monkeys must choose which is the preferred direction of movement of the visual stimuli.The stimuli can display any coherence level from 100% right to 100% left movement (right panel).(B) Circuit schematic of the simulation.Two external stimuli representing motion information enter the excitatory populations of V1. (C) Activity of selected cortical areas during the decision-making task, with a 5%-coherent selective visual input duration of 700 ms, displaying a sensory accumulation process distributed across areas.(D) Sensory accumulation activity in LIP, for two different coherence

Figure 3 :
Figure 3: Winning onset analysis reveals information flow during decision making.(A) Normalized firing rate difference between both excitatory populations for selected cortical areas under zero (left) and full coherence (right).

Figure 4 :
Figure 4: Comparison between MT lesion and LIP lesion in decision making.(A) Schematic of the inactivation protocol.(B) Psychometric curves for baseline, MT lesion and LIP lesion in macaques (left) and in our computational brain model with global coupling strength of G= 0.42 and   = . (right).(C) Relationship between coherence and fraction of valid trials, for different global coupling strength values.(D) Illustration of winning rate, reaching onset, reaching rate, and ramping period.(E) Relative slope change (percentage of ramping speed change) across individual areas when lesioning MT or LIP, averaged over 300 trials.

Figure 5 :
Figure 5: Model prediction of unrecorded cortical areas.(A) Schematic of the inactivation protocol for the neocortex.(B) Slopes of psychometric curve for baseline (0.45) and single area inactivation.(C) Inactivation matrix, with rows indicating the lesioned area, columns the affected areas, and color the relative slope change.(D) Mean of relative ramping speed changes when lesioning individual areas.(E) Relationship between lesioning effect (mean relative slope change) and structural connectivity.The distribution of corresponding variables is also displayed (left and bottom).

Figure 6 :
Figure 6: Temporal gating effects during decision making.(A) Decision making task with distractor.The model was presented with a 100%-coherence stimulus of strength 0.1 and duration 250ms, and then a distractor was presented in one of three time points (sample, early delay, or late delay).(B) Distractors in this task can be considered as either fixed onset (here at 125 ms) and variable duration, or as fixed duration (here, 250 ms) and variable onset.(C) Behavioral impact was quantified as the fraction of trials in which the network changed its decision due to the distractor, for the case of a weak (I=0.17,top panel) or strong distractor (I=0.80,bottom panel) and G=0.42.(D) Minimal effective distractors in terms of duration and strength, for the fixed duration (left) and fixed onset (right) versions of the task, global coupling G=0.42 and input noise of   = ..(E) Same as panel D, but for a network with a stronger global coupling G=0.52 and   = ..

Figure S2 :
Figure S2: Evidence accumulation is distributed across multiple cortical areas.Integration of sensory evidence across all 40 cortical areas in the model, for a visual input of 5% coherence entering V1.We can distinguish between sensory-driven areas (such as V1 and V2), accumulators (like TEO, LIP and 9/46D) and classificators (like 45A and 24c) depending on their response profile and speed of integration.

Figure S4 :
Figure S4: Evidence accumulation at prefrontal cortex (9/46v) during the decision-making task for the simplified three-area model (V1, MT, 9/46v).As in the full-network model, trials with low coherence lead to longer decision times.

Figure S5 :
Figure S5: Lesioning effect on behavioral performance for strongly connected networks.Psychometric curve of the full-network model for the control case (baseline), inactivation of MT (green) and inactivation of LIP (blue)for global coupling of G=0.52.Due to the strong interaction between cortical areas, lesioning either MT or LIP does not lead to substantial performance drops in the task, in contrast to the case of weakly connected networks described in the main text.

Figure S6 :
Figure S6: Temporal gating effects.Impact of temporal gating of distractor signals (purple) on the maintenance of cue information (blue) for V1 and TEpd.Top and middle-top: effects of changing the duration from 280 ms to 270 ms in a fixed onset task.Middle-bottom and bottom: effects of changing the distractor strength in a fixed duration task.