arXiv CSopen access

SeekBrain: An Autonomous Multi-Agent System for Accelerating Neuroscience Discovery

SeekBrain: An Autonomous Multi-Agent System for Accelerating Neuroscience Discovery Jiamin Wu1, 3 * , Peishan Xiang1 * , Jingyang Chen1 † , Yuqing Zhu1 † , Yuxi Li2 † , Ling Luo1 † , Qihao Zheng1 † , Jialiang Zu2 , Yongchao Wu4 , Mindong Liu1 , Haitao Wu1 , Chaofan Hu1 , Yijie Sun1 , Yuqi Hang6 , Yu Zhu1 , Shuo Li1 , Yue Fan1 , Shiyang Feng1 , Wanghan Xu1 , Tianlei Zhang2 , Jie Zhang5 , Wenlong Zhang1 , Bo Zhang1 , Kai Wang2 , Lei Bai1 ♠ , Mianxin Liu1 B♠ , Wanli Ouyang1, 3 B , Jiulin Du1, 2 B , Chunfeng Song1 B♠ 1 Shanghai Artificial Intelligence Laboratory, 2 State Key Laboratory of Brain Cognition and Brain-Inspired Intelligence

Technology, Center for Excellence in Brain Science and Intelligence Technology, Chinese Academy of Sciences, 3 The Chinese University of Hong Kong, 4 Southeast University, 5 Fudan University, 6 New York University

arXiv:2607.29347v1 [cs.MA] 31 Jul 2026

∗ Co-first Author

† Core Contributor

♠ Project Lead

B Corresponding Author

Modern neuroscience relies on integrating multi-scale, multimodal datasets to uncover the neural principles underlying intelligence. However, analytical challenges posed by highly heterogeneous data and fragmented workflows increasingly constrain discoveries. Here we introduce SeekBrain, an autonomous multi-agent framework designed to accelerate neuroscience discovery through domain-grounded hierarchical planning and cross-modal data analysis. SeekBrain dynamically constructs a repertoire of analysis recipes extracted from code-paper pairs. By coupling this codified expertise with agentic planning and execution engines, the framework scalably generates hypotheses and analytical pipelines on demand. Systematic evaluation on the expert-annotated BrainArena benchmark demonstrates that SeekBrain substantially outperforms state-of-the-art agent baselines across various analysis tasks. Crucially, when deployed in real-world research, SeekBrain integrated behavioral, neural, and anatomical data to reveal structured, distributed neural representations of larval zebrafish behavior and a shared axis of regional decoding strength across the brain in a mouse decision-making task. These results establish SeekBrain as a scalable and practical tool for accelerating data-driven discoveries in neuroscience. Project Page: https://ai4neurolab.github.io/SeekBrain-page/ Code: https://github.com/AI4NeuroLab/SeekBrain

1. Introduction Decoding biological intelligence is a central challenge in modern science. To achieve this, neuroscience research must link multiple biological scales to uncover the mechanisms underlying neural computation. Recent experimental technologies now allow us to observe the brain across molecular profiles, structural connectivity, neural dynamics and behavioral [1, 2, 3, 4, 5, 6, 7, 8]. Yet, despite the increasing availability of these high-dimensional datasets [1, 3, 4, 5, 6], progress is constrained by the complexity of analyzing them. Extracting neuroscientific insight from such diverse data requires specialized expertise, cross-disciplinary reasoning, and extensive trial and error. Consequently, routine data processing and coding consume much of the research cycle, leaving many cross-scale hypotheses untested and valuable data underused [2, 7, 8]. Accelerating discovery in neuroscience therefore demands automated frameworks that scale domain expertise to drive an iterative cycle of hypothesis exploration and rigorous data analysis. To address similar bottlenecks in other scientific domains, Large Language Model (LLM)-driven agents have recently emerged to automate key research stages [9, 10, 11, 12, 13, 14, 15, 16, 17]. They have demonstrated the capability to generate testable biomedical hypothesis [11, 18], identify drug repurposing candidates [9], analyze single-cell RNA (scRNA) data [12], and optimize empirical # [email protected]

scientific software [13]. AI-assisted discovery has primarily advanced in fields supported by robust computational tools and unified data representations. In bioinformatics, for instance, frameworks such as CellVoyager automates repetitive scRNA-seq analysis tasks by leveraging established software packages (e.g., scanpy [19]) that encode decades of methodological expertise. Similarly, generalpurpose scientific agents like Biomni [20] rely on clean, structured input data ready for direct computational analysis. However, transferring these agent systems to neuroscience remains severely bottlenecked. The underlying challenge is a mismatch between current agent architectures and the fragmented nature of neuroscience research. Unlike the standardized pipelines in bioinformatics, neuroscience presents a highly heterogeneous landscape where analytical procedures and experimental datasets vary drastically across molecular, structural, functional and behavioral scales. Analyzing such complex data requires deep domain knowledge to design appropriate analytical workflows and interpret findings across different scales. Crucially, the expertise necessary to ensure analytical validity [21, 22] is not centralized in unified software, but instead scattered across disparate publications and isolated code repositories. Without an architecture that explicitly embeds this distributed domain knowledge into research workflows, general-purpose agents remain prone to methodological hallucinations [23, 24], preventing AI systems from reliably automating neuroscience research. Here, we introduce SeekBrain, an autonomous multi-agent framework designed to accelerate neuroscience discovery through domain-grounded research planning and cross-modal data analysis. To scale domain expertise across heterogeneous data modalities, SeekBrain dynamically constructs a Neuroscience Analysis Repertoire. Instead of relying on static toolkits, the framework crystallizes methodological knowledge on demand from relevant literature and open-source codebases, encoding it into adaptable analysis recipes ranging from preprocessing to advanced modeling. A dual-engine architecture, comprising a Research Planning Engine and an Analysis Execution Engine, couples this structured recipe with generalized LLM reasoning to produce hierarchical research plans and orchestrate analytical pipelines across diverse neuroscience tasks. A researcher-in-the-loop paradigm enables experts to provide feedback that could be distilled into refined recipes and reintegrated into the repertoire, allowing the system to generalize to previously unseen challenges. We evaluated SeekBrain using BrainArena, an expert-annotated benchmark comprising diverse data analysis tasks across species and modalities. Our framework outperformed the state-of-the-art agent frameworks Claude Code and Codex by 11.7% and 17.7%, respectively. Furthermore, we demonstrated the system’s capacity for scientific discovery in two real-world research scenarios. In researcher-steering mode, it analyzed an multimodal dataset from free-moving zebrafish, revealing low-dimensional neural representations of behavior with structured manifold geometry and distributed brain-wide tuning. In a SeekBrain-led rediscovery task using a mouse decision-making dataset [6], it formulated and tested a region-level decoding hypothesis, identifying a shared low-rank axis of regional decoding strength across five decoder maps, with the stimulus, choice, and feedback maps retaining spatial variation after regressing out speed and velocity. Together, these results show that SeekBrain can propose hypotheses, perform cross-modal analyses, and facilitate neuroscience discovery, providing a promising path toward AI systems that accumulate and apply domain expertise to uncover the neural principles underlying intelligence.

2. Results 2.1. SeekBrain methodology The SeekBrain is a domain-specific multi-agent system that accelerates neuroscience discovery by formulating research plans and performing cross-modal data analyses in neuroscience (Fig. 1). To

2

a Input

Output

SeekBrain Workflow

Research Query Research Plan

Do mice use trial history to form subjective priors and integrate them with sensory evidence during visual decisions?

Research Planning Engine

Multimodal Data

Analysis Task Sub-question Question

Crystallization Agent Paper Code

Structure

Analysis Execution Engine

…

Neural recordings Transcriptomics

Route Extract

Behavior

Analysis Recipe

Diverse species

b

Figures

Report Agent

Neuroscience Analysis Repertoire Registration & Alignment

Preprocess Structured Synthesis

Report

Analytical Evidence

Feature Extraction

Cross-modal Analysis

Computational Modeling

Visualization

c Analysis Execution Engine

Research Planning Engine

Research Query

Survey Agent

Orchestrator Agent

Reviewer Agent

Analysis Task GLM priors correlate with ideal observer but slope is sub-optimal.

Formulation Agent

d

e

Code Agent

Validation Agent

Researcher-in-the-loop Paradigm

Neuroscience Analysis Repertoire

Customize

Crystallization Agent

Molecule

Gene-expression profiling Cell-type profiling Developmental trajectory

Initialize recipes

Multimodal Analysis Recipes

Molecule-Structure Coupling

Neuroscience Analysis Repertoire

Cell localization in atlas Molecular-to-atlas mapping Cross-modal interpretation

Structure

Connectome analysis Neuron morphology Anatomical quantification

Structure-Function Coupling

SC-FC association testing Network-based S-F coupling Gradient-based S-F comparison

Neural Dynamics

Neural activity profiling Population dynamics FC & network dynamics

Stimuli-Dynamics Coupling

Sensory classification Stimulus decoding Sensorimotor activity coupling

Behavior

Behavioral motif discovery Latent state inference Behavioral strategy modeling

Dynamics-Behavior Coupling

Dynamics-behavior coupling Continuous behavioral decoding Decision-variable decoding

Inject recipes

Incremental update

Unimodal Analysis Recipes

Interpretation Agent

SeekBrain Workflow

Feedback

Researcher

Review/approve Validated User Traces

Distillation Agent

Figure 1 | Overview of SeekBrain. a, Overall workflow. SeekBrain processes research query along with neuroscience data through two core engines. The Research Planning Engine generates a treestructured hierarchical research plan, while the Analysis Execution Engine performs the planned analysis tasks and produces figures and statistical evidence, which are synthesized into a report. Methodological rigor is enforced by the Neuroscience Analysis Repertoire, which contains a wide range of domain-specific analysis recipes crystallized from published studies and their associated codebase. b, Research Planning Engine. After receiving a broad research query, the engine conducts an iterative debate between the Orchestrator and the Reviewer Agents to decompose the query into smaller tasks. This process yields a tree-structured research plan where terminal nodes represent specific analysis tasks for downstream execution. c, Analysis Execution Engine. A Formulation Agent converts each analysis query into a structured computational schema, which the Code Agent then translates into Python scripts. A Validation Agent evaluates the outputs across four core dimensions. This closed-loop refinement process continues until the results pass validation. d, Taxonomy of the Neuroscience Analysis Repertoire, which includes both unimodal and multimodal analysis recipes. Representative recipe types are presented. FC, functional connectivity; SC, structural connectivity. e, SeekBrain’s researcher-in-the-loop paradigm, which enables users to flexibly guide the analysis process and supports the continuous refinement of analysis recipes based on expert interaction traces. 3

facilitate comprehensive scientific research, this framework supports a broad spectrum of multiple data modalities, including behavior, neural activity, structural connectivity, and transcriptomic/genomic profiles, which enables diverse needs in cross-modal analytical exploration (Fig. 1a, left). To ensure methodological rigor, SeekBrain is anchored by a Neuroscience Analysis Repertoire, which contains analysis recipes derived from published studies and refined through expert input. This repertoire crystallizes domain knowledge into reusable analysis recipes that guide automated data analysis (Fig. 1a, middle). Using these recipes, SeekBrain operates through two agent-based engines: a Research Planning Engine, which breaks each general research query into specific analysis tasks organized in a hierarchical plan (Fig. 1b), and an Analysis Execution Engine, which orchestrates analytical workflows, producing empirical results and scientific findings through iterative cycles of code generation, execution, and validation (Fig. 1c). Researchers interact with SeekBrain through a natural-language interface and can guide the analysis process flexibly, establishing a collaborative researcher-in-the-loop paradigm (Fig. 1e). Neuroscience analysis repertoire. Data analysis in neuroscience requires specialized knowledge of data processing, modeling, statistics, and visualization. However, this knowledge is often dispersed across research papers, code repositories, and expert analysis traces rather than organized into standardized, reusable workflows. This makes it difficult for general-purpose agents to develop reliable analysis procedures from task instructions alone. To combine information from these sources, we developed an LLM-driven crystallization agent. The agent automatically creates analytical recipes by extracting the scientific context from research papers and the implementation details from their associated code repositories. Each recipe is a reusable procedure that links a scientific aim to the required data, processing and modeling steps, visualization methods, and corresponding code. Together, these recipes form the Neuroscience Analysis Repertoire, allowing SeekBrain to automatically retrieve, adapt, and combine them for new research questions. The initial repertoire contains 68 foundational recipes covering analysis stages from data preprocessing to visualization (Supplementary Table 1). These recipes are organized by data modality into unimodal recipes for neural dynamics, behavior, structure, and molecular data, and multimodal recipes that combine two or more modalities (Fig. 1d). The repertoire is also designed to improve through use. After an adapted workflow has been reviewed and validated, its execution trace can used to update the corresponding recipe or create a new variant. This co-evolving process allows the repertoire to retain useful expert decisions and accumulate analysis experience over time. Research planning engine. To handle high-level requests, our hierarchical planning engine performs scientific problem decomposition through a debate between an Orchestrator and a Reviewer agent (Fig. 1b). With scientific background provided by a Survey Agent, the Orchestrator breaks the initial query down into targeted sub-questions (working hypotheses), which are further resolved into atomic analysis tasks. This creates a tree-structured plan: the root node represents the original query, while terminal nodes represent specific tasks mapped to distinct biological hypotheses. This structure makes the reasoning process transparent by linking evidence from each branch directly to the main question. To ensure scientific rigor, the Reviewer assesses the plan’s relevance to the research questions, logical coherence, data adequacy, and redundancy across tasks. The Orchestrator then refines the plan based on this feedback. This review cycle continues until both agents agree on a concrete research plan that connects the biological question to downstream analyses. Analysis execution engine. To autonomously execute each atomic analysis task in the research plan, the Analysis Execution Engine orchestrates a four-phase workflow: formulation, coding, validation, and interpretation (Fig. 1c). First, a Formulation Agent evaluates the assigned task alongside the 4

available data to generate a methodological schema. This schema details the necessary processing steps, algorithms, and visualization strategies. By retrieving and adapting recipes from the Neuroscience Analysis Repertoire, the agent aligns its approach with established domain knowledge, minimizing unsupported methodological choices. Next, a Code Agent translates this schema into Python scripts. Through an iterative cycle of data exploration, parameter tuning, and error resolution, it generates figures and statistics to produce a reliable analytical workflow. To prevent methodological hallucinations, a Validation Agent acts as an internal reviewer. It assesses the outputs for inquiry alignment, methodological fidelity, quantitative integrity, and visual quality. This feedback loop continues until the results are scientifically sound. Finally, an Interpretation Agent analyzes the validated figures and statistics to formulate a self-contained interpretation summary. Throughout this process, the system transparently reports any limitations, explicitly documenting uncertainty whenever the evidence or diagnostics prove insufficient. Researcher-in-the-loop paradigm. SeekBrain implements a researcher-in-the-loop paradigm (Fig. 1e), enabling scientists to inject feedback at every stage of the workflow. Users can refine the initial research plan, assess methodological feasibility, and fine-tune analytical pipelines. Crucially, SeekBrain captures these user traces after each session and distills them into updated analysis recipes. By seamlessly reintegrating them into the Neuroscience Analysis Repertoire, the system continuously evolves with the broader research community.

2.2. Assessing SeekBrain’s analysis capability with BrainArena Introducing BrainArena. Autonomous discovery in neuroscience requires LLM agents to independently execute end-to-end data analysis. To evaluate this capability, we introduce BrainArena, a benchmark designed to test how well these agents handle complex analyses of multimodal neuroscience data (Fig. 2a). BrainArena includes 32 tasks constructed from selected neuroscience studies (Supplementary Table 2). Each task includes an expert-annotated analysis query and a scoring rubric (Fig. 2b). The query sets the scope of the task by specifying the input data, overall objective, and required output. Each rubric is based on the original study and includes three to five scoring items, with the total score set at 100 points. The rubric evaluates three aspects: method correctness, figure fidelity (whether the visual patterns match the reference), and interpretive validity (whether the results are interpreted in a scientifically valid way). Across the 32 tasks, the rubrics contain over 700 deduction subitems that map observable errors to point deductions. BrainArena includes data from five common model organisms (Fig. 2c): C. elegans (16%), Drosophila (6%), zebrafish (31%), mouse (34%), and macaque monkey(12%). It covers diverse data modalities, including neural activity, behavioral, anatomical and structural, and molecular and genetic data. Most tasks integrate multiple modalities: 28% use a single modality, whereas 59%, 9%, and 3% combine two, three, and four modalities, respectively; behavioral and neural activity data constitute the most common combination (53%; Fig. 2d). These tasks span seven types of analysis: neural response characterization, neural decoding, population representation analysis, computational modeling, statistical analysis, brain network and anatomical mapping, and molecular and developmental analysis. Based on task complexity, the tasks were classified as easy or hard. Together, BrainArena combines broad coverage of common model organisms, data modalities, analysis types, and difficulty levels with expert-annotated rubrics, making it a rigorous domain-specific benchmark for neuroscience data analysis. Benchmark results. We evaluated SeekBrain against two leading agent systems on BrainArena: Claude Code (Claude Opus 4.7 [25]) and Codex (GPT-5.5 [26]). SeekBrain uses Claude Opus 4.7 5

a

c

b

d 3%

9% 6%

Macaque

Mouse

9% 6%

Triple

Single

6% 6%

C. elegans Dual

Zebrafish

53% M = Molecule, B = Behavior, N = Neural activity, S = Structure Species Mouse (34%) Modality count Combination detail Zebrafish (31%) Single modality (28%) M (9%) B + N (53%) C.elegans (16%) Dual modality (59%) B (6%) M + S (6%) Macaque monkey (12%) N (6%) Triple modality (9%) B + N + S (9%) Drosophila (6%) Quadruple modality (3%) S (6%) B + M + N + S (3%) 53%

e

f ***

Population Repr. Analysis Computational Modeling

50

Statistical Analysis Brain Network & Anatomical Mapping

25

Molecular & Developmental Analysis

SeekBrain

SeekBrain 100

Codex

Claude Code

Codex

*

*** ***

75

Mean score

Claude Code

50

25

0

Easy

Hard

0

h

SeekBrain 100

Mean normalized score (%)

g

Codex

Neural Decoding

75

0

Claude Code

Neural Response Characterization

75

**

20

Claude Code

**

40

Codex

80

i

100

****

100

*** *

50

75

50

25

0

60

Mean score

Mean score

Mean score

100

SeekBrain

****

Figure fidelity Method correctness Interpretive validity

25

w/ recipe

w/o recipe

Figure 2 | BrainArena evaluation. a, Overview of BrainArena construction and evaluation process. Expert-curated tasks are generated from neuroscience publications and evaluated with ground-truth annotations and scoring rubrics. b, A task example, including the task query and evaluation scoring rubrics annotated by domain experts. c, Tasks in BrainArena span diverse model organisms in neuroscience. d, BrainArena includes single- and multi-modality analysis tasks across molecule, behavior, neural activity, and structure modalities. 6

Figure 2 | e, Performance (mean scores) of SeekBrain and other LLM agents on BrainArena. Points indicate scores on individual tasks. Differences among the three agents were assessed using a Friedman test. Following a significant result ( 𝑃 < 0.05), post hoc Dunn’s multiple-comparisons test was used to compare SeekBrain with Claude Code and Codex. Error bars indicate s.e.m. g, Mean scores grouped by task difficulty: easy-level (12 tasks) and hard-level (20 tasks). Within each difficulty level, statistical analyses were performed as described for panel e. h, Component-level performance across rubric dimensions. Mean normalized scores (% of maximum rubric points) for each rubric component are presented. The three components assess figure fidelity (visual pattern consistency with the reference), method correctness, and interpretive validity (scientific soundness of result interpretation). Within each rubric component, statistical analyses were performed as described for panel e. i, Performance comparison between with and without recipes from neuroscience analysis repertoire. Wilcoxon signed-rank test. Statistical significance: * 𝑃 < 0.05; ** 𝑃 < 0.01; *** 𝑃 < 0.001, **** 𝑃 < 0.0001.

as base LLMs for most agents (details in Section 4). The Research Planning Engine is ablated in this experiment. Certain AI scientist frameworks are excluded from this comparison because of architectural incompatibilities: Robin [9] and AI Scientist [10] do not support native-format scientific data, whereas AI Co-scientist [11] focuses on theoretical ideation and hypothesis generation. Under these settings, SeekBrain achieved the highest mean score (75.8), compared with 64.1 for Claude Code and 58.1 for Codex; performance differed significantly among the three agents (Friedman test, 𝑃 < 0.0001), and post hoc Dunn’s multiple-comparisons test showed that SeekBrain outperformed both Claude Code (adjusted 𝑃 = 0.0002) and Codex (adjusted 𝑃 < 0.0001; Fig. 2e). SeekBrain also achieved the highest mean score in all seven analysis types (Fig. 2f and Supplementary Figure 1). Its largest gains over Codex were in neural decoding (28.7 points) and brain network and anatomical mapping (28.5 points). When the tasks were grouped by difficulty, SeekBrain’s advantage was greater on hard tasks than on easy tasks (Fig. 2g). Performance declined more on hard tasks for the two baseline agents, whereas SeekBrain remained more consistent across the two difficulty levels. Score decomposition showed that SeekBrain outperformed both baselines across all three evaluation dimensions (Fig. 2h). Its largest advantage was in interpretive validity, which requires accurate statistical inference and biologically grounded interpretation beyond code execution. To further validate the role of the analysis recipes, we removed them from the workflow in an ablation study. This ablation led to issues in key analysis steps (Supplementary Table 4) and reduced the mean score from 75.8 to 63.2 (Wilcoxon signed-rank test, 𝑃 < 0.0001; Fig. 2i). For example, in a network visualization task, the agent failed to correctly align the spatial axes of brain region projections with those of the network node graph (Supplementary Figure 2). Overall, these results suggest that SeekBrain’s domain-grounded engine automatically performs diverse, end-to-end analyses, consistently producing faithful figures and scientifically meaningful interpretations.

2.3. Case study 1: SeekBrain analyzes a tri-modal zebrafish dataset through a sequence of researcher-led questions To show how SeekBrain supports complex scientific investigations that require multiple linked stages of analysis, we applied it to a multimodal dataset capturing brain-wide neuronal calcium dynamics and free-swimming behavior in larval zebrafish. We used SeekBrain to perform integrated analyses linking behavior, neural activity, and anatomy. At each stage, the researcher posed a question based on earlier results, and the agent ran the corresponding analysis pipeline. First, to determine whether swimming bouts could be mapped into a behavioral space based on tail kinematics, the researcher asked the agent to run an adapted behavioral clustering pipeline 7

Figure 3 | Structured and distributed brain-wide neural representations of behavioral types in free-moving larval zebrafish. a, Bout-level behavioral space. Each point denotes one bout (𝑛 = 1,546 bouts from 18 larvae). Colors denote the three bout clusters. Faint points indicate long-duration bouts (>500 ms). Yellow markers a–i locate 9 representative skeletons in Panel b. b, Tail features of 3 bout clusters. For each cluster, the upper trace plot shows all member bouts’ tail-tip angle trajectories; the kinematic heatmap displays the mean time-varying angles of 9 tail segments from rostral (top) to caudal (bottom), using a fixed color scale spanning [−90◦ , +90◦ ]; the right skeleton plot exhibits the representative tail posture in different bout clusters. ros., rostral; cau. caudal. c, Region-aware spatial-functional supervoxel (SV) topology from an example larva. Each point is one of 3,000 supervoxel centroids. Colors indicate major anatomical blocks. Forebrain, orange; midbrain, green; hindbrain, blue; spinal cord, purple; ganglia, red. d, Example calcium activity heatmap of brain-wide supervoxels sorted by major anatomical block and region label. Bout-locked activity peaks are denoted at the top of the heatmap. S., spinal cord; G., ganglion. e, InfoNCE training loss of the bout type supervised CEBRA model fitted in one example larva. Dashed line, uninformative-baseline reference. f, CEBRA latent manifold for the same larva in Panel e. Gray, full-session 8D latent trajectory projected onto a 2D plotting basis derived by SVD of bout-mean latents; colored points, accepted bout-mean latents colored by bout type. g, Aggregated 3NN confusion matrix across larvae. 8

Figure 3 | h, Cross-fish cluster centroids aligned by orthogonal Procrustes. Within each fish, the mean cosine distance to same-cluster centroids of other fish was compared to the mean cosine distance to different-cluster centroids. Boxes show fish-level summary statistics; overlaid points are fish; lines connect paired same-cluster and different-cluster values for each fish. Paired Wilcoxon signed-rank test across 𝑛 = 13 larvae, 𝑝 = 2.44 × 10−4 . Diff., different. i, Pairwise cosine-distance matrix among aligned fish-by-cluster centroids, grouped by cluster and fish. Color encodes pairwise cosine distance on the aligned 8D centroid vectors. j, PCA projection of the same aligned centroids, colored by bout cluster type. k, Region-aware significant supervoxels for each bout type. Top and side brain projections; gray points show all supervoxel centroids pooled across larvae; colored points show supervoxels belonging to spatial clusters that were significant in the 𝐾 -vs-rest contrast for each type. Sig. SVs, significant supervoxels. l, m Per-region significant-supervoxel counts and ratios across larvae. The stacked display is visual only. Ratios for clusters are computed independently, not mutually exclusive. See Supplementary Table 6 for the full names of the brain region abbreviations.

based on an established framework [27]. The agent generated and ran all scripts required for this multistep analysis. The analysis included decomposing time-varying tail-segment angles into principal components (PCs) of posture, calculating polarity-invariant dynamic time warping (DTW) distances, embedding the resulting distance matrix with Isomap, and identifying bout types with Ward’s hierarchical clustering. Consistent with previous studies, larval swimming bouts occupied a continuous behavioral space (Fig. 3a). Across the tested values of 𝐾 , the silhouette score was highest at 𝐾 = 3, supporting a three-cluster partition (Fig. 3a,b). Cluster 1 contained large-angle turns, Cluster 2 corresponded to J-turns, and Cluster 3 consisted mainly of low-amplitude forward swimming (Fig. 3b). To link high-dimensional, brain-wide neural dynamics to these behavioral clusters, we next sought to identify a latent neural manifold and isolate candidate brain regions associated with each cluster. The raw dataset contained 64,032 ± 8,575 ROIs per fish (mean ± s.d.). To reduce the data size while retaining spatial and functional information, the agent grouped ROIs across the brain into region-aware supervoxels (SVs). Within each anatomical region, k-means clustering was applied to a feature matrix that combined PCs of neural activity with the spatial coordinates of the ROIs. This procedure produced a fixed total of 3,000 SVs per larva (Fig. 3c). The activity traces from these region-aware SVs were then aligned to bout epochs for subsequent cross-modal analyses (Fig. 3d). To test whether a low-dimensional representation representation of brain-wide activity could capture the three behavioral clusters, the researcher asked the agent to train a CEBRA model supervised by bout type [28] using behaviorally relevant frames. The trained model was then applied to the full time series from each fish (Fig. 3e,f). For every fish, training converged to an InfoNCE loss below the uninformative baseline; the final loss was 5.703 ± 0.070 (mean ± s.d.). When averaged within bouts, the resulting eight-dimensional (8D) neural representations formed clearly separated clusters (Fig. 3f). In cross-validation, a 3-nearest-neighbor (3NN) classifier decoded the behavioral cluster labels from the CEBRA representations with an accuracy of 0.949 ± 0.034 (mean ± s.d.; Fig. 3g). These results show that the learned neural representation captured behavior-related differences within individual fish. To compare the latent cluster geometry across fish, we asked the agent to align the cluster centroids using orthogonal Procrustes analysis. Only clusters with at least three bouts were included. After alignment, the agent calculated, for each fish, the mean cosine distances between its cluster centroids and matching or nonmatching centroids from the other fish. The mean distance between matching centroids was significantly smaller than that between nonmatching centroids (paired Wilcoxon signedrank test, 𝑝 = 2.44 × 10−4 , 𝑛 = 13 larvae; Fig. 3h). This shared cluster organization across fish was 9

Figure 4 | A shared axis of regional decoding strength and task-specific residual patterns in a mouse decision-making task. a, Pearson correlation matrix of stimulus, choice, feedback, speed, and velocity decoder maps across 201 regions. b, Spearman rank-correlation matrix for the same five decoders. c, Observed PC1 variance fraction (red line; 0.632) compared with a 5,000-permutation column-shuffle null; empirical 𝑃 ≈ 0.0002. d, PC1 loadings for the metric-value PCA and rank-PCA decompositions; all decoder loadings are positive. e, 𝑅2 from ordinary least-squares regression of each task decoder on speed and velocity. f, Main PC1 scores for the 201 included regions, rendered bilaterally and mirror-symmetrically on the Swanson mouse-brain flatmap. Higher scores occurred in reticular, motor-related brainstem, and thalamic regions. g, Rank-PCA PC1 score projected onto the same atlas, confirming robustness of the shared axis. h, Choice-decoder residual map after removing speed- and velocity-linked variance. i, Stimulus-decoder residual map after removing speed- and velocity-linked variance. j, Feedback-decoder residual map after removing speed- and velocity-linked variance. See Supplementary Table 6 for the full names of the brain region abbreviations. also visible in the pairwise cosine-distance matrix and the PC projection of the aligned centroids (Fig. 3i,j). Finally, to map the spatial distribution of neural activity associated with each behavioral cluster, we asked the agent to run a region-aware permutation test to identify SVs whose bout-averaged calcium activity differed significantly between that cluster and the other two clusters combined (Fig. 3k). These SVs were distributed throughout the brain, from the forebrain to the hindbrain. They were found in 47, 39, and 43 of the 81 anatomical regions for Clusters 1, 2, and 3, respectively (Fig. 3l,m). Some regions contained a relatively high proportion of tuned SVs, including the thalamus, pretectum (PT), several hypothalamic subregions, the nucleus of the medial longitudinal fasciculus (nMLF), the superior and inferior raphe nuclei, and several hindbrain rhombomeres (Fig. 3m). Taken together, this case study discovered the structured and distributed brain-wide neural representations of behavioral types in free-moving larval zebrafish. Throughout this data exploration, SeekBrain reliably carried out the required analyses, helping the researcher investigate the neural basis of behavior.

10

2.4. Case study 2: SeekBrain reanalyzes a mouse decision-making dataset through agent-led hypothesis testing While the first case study focused on SeekBrain’s ability to carry out a sequence of analyses defined by the researcher, the second tested whether it could propose testable hypotheses and plan analyses to evaluate them. We therefore used SeekBrain to reanalyze the published International Brain Laboratory (IBL) mouse decision-making dataset [6] in an agent-led process supervised by the researcher. SeekBrain proposed candidate hypotheses and analysis plans, while the researcher selected promising directions and gave specific guidance on implementation. Findling et al. used this dataset to examine how task priors are represented in the brain during visual decision-making [29]. The mice inferred the changing probability of a stimulus appearing on either side and used this information to guide their choices. This internal estimate, termed the subjective prior, was encoded in 20–30% of brain regions. The widespread representation of this prior was consistent with distributed recurrent Bayesian inference, rather than with prior integration limited to decision-making regions. SeekBrain proposed another region-level hypothesis using the same dataset: decoder performance for stimulus, choice, feedback, speed, and velocity might share a dominant pattern across brain regions rather than vary independently. To test this hypothesis, the agent used data from 279 unique Beryl-labeled regions in the original region-level table from Findling et al. [29]. It retained the 201 regions with finite values for all five decoder variables, forming a 201 × 5 region-by-decoder matrix. Pairwise correlations first suggested that the five decoder maps shared a common pattern across brain regions. Both Pearson and Spearman analyses showed generally positive correlations: regions with high decoding strength for one variable also tended to have high decoding strength for the others (Fig. 4a,b). SeekBrain then standardized each decoder column using z-scores and applied principal component analysis (PCA) to the resulting region-by-decoder matrix. The PC1 explained 63.2% of the total variance. This value was higher than expected under a null distribution generated by independently shuffling the values within each decoder column 5,000 times (empirical 𝑃 ≈ 2.0 × 10−4 ; Fig. 4c). The loadings of all five decoder variables on PC1 had the same sign, supporting a shared axis of decoding strength rather than an opposing pattern among the variables (Fig. 4d). To assess whether the shared-axis pattern depended on the scale of the decoder values or on extreme values, SeekBrain repeated the PCA after replacing the values in each decoder map with their ranks within that map. After matching the orientation of PC1 across the two analyses, the regional PC1 scores from rank-based PCA were strongly correlated with those from the original value-based PCA (𝑟 = 0.96). The two analyses also produced similar spatial patterns of regional PC1 scores in the brain atlas (Fig. 4f,g), showing that the shared axis was robust to differences in scale and to extreme values. High shared-axis scores were found in reticular regions, motor-related brainstem areas, and thalamic regions, whereas several olfactory and piriform-amygdalar regions had lower scores. SeekBrain next asked whether this shared axis reflected the widespread influence of movement on neural activity. It used ordinary least-squares regression to predict each task decoder map from the speed and velocity decoder strengths. These movement variables explained 61% of the regional variance in choice decoding ( 𝑅2 = 0.61), 39% in feedback decoding (𝑅2 = 0.39), and 15% in stimulus decoding ( 𝑅2 = 0.15; Fig. 4e). After accounting for movement, the residual maps still showed distinct regional patterns for choice, stimulus, and feedback (Fig. 4h-j). Thus, although movement was associated with the shared decoding pattern, it did not fully explain the regional patterns specific to each task variable. In summary, this case study demonstrates that SeekBrain can propose and test new hypotheses using published datasets, helping researchers explore existing data in new ways.

11

3. Discussion In this work, we developed SeekBrain, a specialized multi-agent framework that accelerates discovery in neuroscience through domain-grounded research planning and cross-modal data analysis. Crucially, our approach overcomes a fundamental bottleneck in deploying AI for neuroscience: the severe mismatch between generic agent architectures and the fragmented nature of neuroscientific research. To bridge this gap, we constructed a dynamic Neuroscience Analysis Repertoire that autonomously crystallizes dispersed expertise into analysis workflows. Driven by hierarchical planning and analysis execution engines, this architecture successfully transforms general-purpose LLMs into scientific collaborators capable of proposing plausible exploratory plans and orchestrating analysis workflow. We validated this capability through comprehensive benchmarking and two real-world cross-modal discovery cases, highlighting its potential to automate complex research pipelines. SeekBrain diverges from existing LLM-based scientific agents by prioritizing the internalization of domain-specific methodology. Whereas frameworks like AI Co-Scientist [11] or AI Scientist [10] rely on generic reasoning without grounding execution in field-specific constraints, and bioinformatics agents like CellVoyager [12] inherit established software abstractions (e.g., scanpy [19]), our system must construct its expertise largely from scratch. By distilling analysis recipes directly from literature and open-source codebases, our framework creates an indexable, evolvable knowledge base. Future iterations will focus on broadening this repertoire to include spatial transcriptomics, large-scale connectomics, and closed-loop behavior. The BrainArena, as a domain-specific benchmark to evaluate automatic analysis ability on real-world data, complements general academic-reasoning benchmarks such as Humanity’s Last Exam [30] and code-centric benchmarks such as SWE-bench [31] and MLE-bench [32]. Its 32 expert-annotated tasks are scored using rubrics comprising more than 700 deduction subitems across figure fidelity, method correctness, and interpretive validity, making BrainArena a rigorous benchmark for assessing agents against the demands of end-to-end scientific workflows. On this benchmark, SeekBrain outperformed Claude Code and Codex across all seven analytical domains, with its advantage becoming more pronounced on harder tasks. Among the three scoring dimensions, its largest advantage was in interpretive validity, which assesses whether agents can draw scientifically valid and meaningful conclusions from their analytical results. Ablating the repertoire from the framework caused a significant performance drop of over 10 points, demonstrating that agents struggle to derive valid biological interpretations without explicit domain knowledge to enforce scientific rigor. Beyond the quantitative benchmark, two real-world case studies demonstrate SeekBrain’s capability to drive scientific exploration in different modes. Under researcher guidance, SeekBrain analyzed the multimodal dataset from free-moving larval zebrafish. Building on Mearns et al.’s study of swimming behavior [27], it examined how different types of swimming bouts were represented across the brain. Through a sequence of instructed analyses, SeekBrain revealed a low-dimensional structure linking neural activity and behavior. In a more autonomous setting, SeekBrain reanalyzed decoder outputs [29] from the IBL mouse dataset [6] and tested whether decoding strength shared a common pattern across brain regions. Previous studies had shown that movement broadly affects neural dynamics [33, 34, 35]. On this background, SeekBrain represented five published region-wise decoder maps (stimulus, choice, feedback, speed, and velocity) along a shared low-rank axis of decoding strength. It then tested whether this shared axis was related to movement and separated the movement-related component from the remaining task-specific regional patterns. Together, these cases demonstrate that SeekBrain supports real-world research across varying levels of automation, enabling researchers to maintain control over key scientific decisions while delegating routine analysis and exploration. More broadly, large-scale brain initiatives and data-sharing efforts of the community have created an expanding ecosystem of open-access, multimodal datasets [1, 3, 4, 5, 6, 36], yet much

12

of which remains underexplored. Equipped with dynamic domain-specific repertoire, SeekBrain serves as a evolving system to generate new hypotheses and systematically reanalyze these public datasets, ultimately helping uncover biological principles hidden within the growing body of neuroscience data resources. Looking ahead, SeekBrain points toward a continuous development loop between AI-driven neuroscience and biologically inspired AI. Within neuroscience, the framework provides a foundation for a scalable, collaborative ecosystem where research groups can seamlessly contribute and curate specialized analytical workflows. Conversely, on the AI development side, extracting the computational principles underlying biological intelligence offers vital inductive biases for next-generation artificial models. We envision that this feedback loop will create a future where intelligent systems accelerate our understanding of the brain, and in turn, use those insights to guide their own evolution.

4. Methods 4.1. Research planning engine We address this limitation by introducing a hierarchical planning workflow that explicitly models this progressive decomposition. Every stage of the process incorporates an automated peer-review mechanism that evaluates scientific validity, data adequacy, and execution feasibility, providing iterative feedback to refine the plan. We use Gemini-3.1-Pro [37] for sub-question generation, question-level synthesis and reporting, whereas Claude-Opus-4.7 [25] is used for review, analysis task generation and sub-question-level conclusion synthesis. Research background survey. Before research planning, a Survey Agent constructs the required scientific background from the literature corpus. Given a research query, the agent identifies relevant publications, synthesizes key findings and unresolved gaps into a concise background summary via retrieval-augmented generation. Literature-grounded planning produces scientifically motivated research plans that build on prior work and target meaningful open questions. The system also supports user-provided background context. Sub-question decomposition and task generation. The Orchestrator Agent converts a broad scientific objective into a structured research plan. It takes three structured inputs: a user-defined scientific objective, the target dataset metadata, and optional background information. If the metadata is absent, a Data Analysis Agent is invoked to autonomously profile the dataset, extracting data structures, variable types, sample sizes, and etc. The Orchestrator Agent then divides the main question into related sub-questions and creates concrete analysis tasks for each one. These subquestions and tasks must be testable, biologically meaningful, non-redundant, and feasible with the available data. The resulting plan has a tree structure in which the original query forms the root node, and the analysis tasks form terminal nodes. Users may specify the number of sub-questions and tasks or allow the system to determine it based on task complexity. Automated iterative peer-review. To ensure the workflow’s robustness, we implemented a multiagent iterative refinement loop following both the sub-question and analysis task decomposition stages. An independent Reviewer Agent evaluates the Orchestrator’s outputs based on three criteria: scientific validity, data adequacy, and execution feasibility. The Reviewer assigns a quantitative score to each criterion. If a proposed question or task falls below the acceptance threshold, the Reviewer synthesizes targeted feedback and revision suggestions. These comments are then fed back into the Orchestrator’s context to guide the next round of generation. This generation-review-revision loop 13

continues until the output meets the acceptance criteria or reaches a predefined limit (e.g., 𝑁 = 3 iterations).

4.2. Analysis execution engine SeekBrain coordinates a set of specialized language agents that convert analysis task query into executable data analysis workflows. The coding and core reasoning agents in this engine are powered by Claude Opus 4.7 [25] with adaptive extended thinking enabled; long-text auxiliary tasks are routed to Gemini 3.1 Pro [37]. Analytical formulation. Given an neuroscience analysis inquiry, we design a Formulation Agent to produce a step-by-step detailed methodological schema. The agent evaluates the inquiry against the empirical data and retrieves analysis recipes from constructed neuroscience analysis repertoire. Under the guidance of analysis recipes, it decomposes the request into a structured sequential methodological schema, explicitly detailing the data inputs, expected output and execution specifications. Acting as a strict computational constraint, the generated schema anchors the subsequent workflow, ensuring downstream agents accurately translate it into executable scripts. Closed-loop generation and validation. To automate the neuroscience data analysis, we design a generation and validation loop to generate robust analysis workflow scripts. In the generation phase, the Code Agent translates the generated analysis plan into executable Python scripts. The Code Agent employs a ReAct (Reasoning and Acting) paradigm [38] to autonomously write, execute, and debug the analytical code under restricted reasoning and execution iterations. The process iterates between data probing, parameter tuning, statistical evaluation, and runtime error resolution before ultimately yielding robust scripts grounded in comprehensive data exploration. The scripts are executed to output scientific figures and key quantitative statistics. In the validation phase, SeekBrain includes a Validation Agent that enforces research rigor by evaluating whether the outputs constitute a scientifically valid neuroscience analysis, which distinguishes it from general-purpose agentic systems. Following each successful execution round, the agent assesses the generated results across four dimensions: (1) alignment with the original scientific inquiry, (2) methodological fidelity to the analysis plan, (3) quantitative data integrity, and (4) visual quality of the generated figures. Any failure to meet these criteria triggers an iterative refinement loop, with the Validation Agent supplying its diagnostic feedback to the Coding Agent. This generation-and-validation loop iterates across several rounds until the analytical outputs achieve full compliance, ensuring the final results are both methodologically robust and biologically meaningful. Result interpretation. The Interpretation Agent evaluates the generated figures and quantitative statistics to extract core scientific insights. It then integrates these insights and full execution traces to formulate a comprehensive scientific summary of the entire analytical process. It explicitly reports method limitations whenever analysis fails, statistical evidence proves insufficient, or numerical diagnostics trigger warnings.

4.3. Report generation Completing the analysis tasks yields a large set of intermediate outputs, including statistics, figures, tables, and text. To handle this, a Report Agent aggregates the evidence hierarchically to form traceable conclusions. It begins by summarizing each analysis task results into key findings, supporting statistics, visuals, and limitations. It then merges the evidence from related plans to answer each sub-question, 14

keeping clear references to the original plans. Finally, the agent combines these sub-question answers to address the main scientific objective. Using this aggregated evidence, the agent generates a formatted academic PDF report with standard sections.

4.4. Construction of the Neuroscience Analysis Repertoire To bridge the knowledge gap between SeekBrain and neuroscience expert, the Neuroscience Analysis Repertoire acts as a dynamic knowledge base by continually expanding through two ways: crystallization analytical procedures from open-sourced literature, and the subsequent distillation of validated execution traces guided by user feedback. Recipe crystallization. To construct the repertoire autonomously, a Crystallization Agent extracts and synthesizes methodological knowledge from paired inputs comprising open-access publications and their corresponding source-code repositories. For a specified analytical objective, the agent extracts task-relevant methodological descriptions from the manuscript and simultaneously scans the repository to identify relevant executable components. The agent then performs a cross-modal alignment, explicitly mapping theoretical analytical steps described in the text to specific functional implementations and code blocks in the repository. Once this correspondence is established, the agent compiles the theoretical context and the executable code into a standardized recipe. This formatted recipe integrates the scientific rationale with specific computational operations, which is subsequently reintegrated into the repertoire for future neuroscience analysis. Distillation of user trajectories. The repertoire could continuously evolves with its users through a Distillation Agent. Upon the successful completion of an analytical workflow, the system invokes a Distillation Agent to process the computational history. The agent receives the complete execution trace, which comprises the finalized code blocks, runtime logs, and any qualitative feedback provided by the human user. To derive a reusable analysis recipe, the Distillation Agent parses these inputs to abstract the core analytical logic. Specifically, the agent is prompted to retain the sequence of methodological operations and expected outputs, while stripping away incidental, run-specific artifacts such as local file paths or dataset-dependent variable names. The abstracted procedure is then formatted into a standardized recipe and appended to the Neuroscience Analysis Repertoire. This update mechanism allows the system to progressively expand its methodological database using validated researcher-in-the-loop trajectories.

4.5. BrainArena Benchmark construction. BrainArena built on source studies with publicly available neuroscience data and reproducible figure-level analyses. The benchmark contained 32 tasks from 10 source studies [39, 29, 6, 40, 41, 42, 43, 44, 45, 46], with each task corresponding to a key figure in the original publication. Each task was annotated with an expert-written analysis query, data paths, a reference figure and a ground-truth scientific conclusion. Domain experts checked the task specification against the source study and constructed a task-specific 100-point scoring rubric. Supplementary Table 3 provides the complete task metadata. Experimental setup. For inference on BrainArena, each system received the same task package, consisting only of the analysis task query and raw data paths. To ensure rigorous and unbiased evaluation, all tested systems were strictly blinded to the source paper identity, reference figures, and ground-truth conclusions. We evaluated SeekBrain, Claude Code and Codex on BrainArena. 15

To mitigate the risk of benchmark data leakage, BrainArena included tasks derived from a bioRxiv preprint that was published after the official release date of the Claude model. The task execution prompt used for benchmark baselines is provided in Supplementary Figure 3. Evaluation. For evaluation, we used Claude Opus 4.6 [47] as a rubric-based LLM judge to assess whether the model outputs reproduced the reference figures and scientific conclusions. The judge received the task query, model output, and returned item-level scores based on task-specific rubrics annotated by experts. The LLM-as-judge evaluation prompt is provided in Supplementary Figure 4.

16

References [1] BRAIN Initiative Cell Census Network (BICCN). “A multimodal cell census and atlas of the mammalian primary motor cortex”. In: Nature 598.7879 (2021), pp. 86–102. [2] Emily S. Finn, Russell A. Poldrack, and James M. Shine. “Functional neuroimaging as a catalyst for integrated neuroscience”. In: Nature 623.7986 (2023), pp. 263–273. [3] Zizhen Yao, Cindy T. J. van Velthoven, Michael Kunst, Meng Zhang, Delissa McMillen, et al. “A high-resolution transcriptomic and spatial atlas of cell types in the whole mouse brain”. In: Nature 624.7991 (2023), pp. 317–332. [4] The MICrONS Consortium. “Functional connectomics spanning multiple areas of mouse visual cortex”. In: Nature 640.8058 (2025), pp. 435–447. [5] Sven Dorkenwald, Arie Matsliah, Amy R. Sterling, Philipp Schlegel, Szi-chieh Yu, et al. “Neuronal wiring diagram of an adult brain”. In: Nature 634.8032 (2024), pp. 124–138. [6] International Brain Laboratory, Dora Angelaki, Brandon Benson, Julius Benson, Daniel Birman, et al. “A brain-wide map of neural activity during complex behaviour”. In: Nature 645.8079 (2025), pp. 177–191. [7] Anton Arkhipov, Nuno da Costa, Saskia de Vries, Trygve Bakken, Corbett Bennett, et al. “Integrating multimodal data to understand cortical circuit architecture and function”. In: Nature Neuroscience 28.4 (2025), pp. 717–730. [8]

Mackenzie Weygandt Mathis and Alexander Mathis. “Joint modelling of brain and behaviour dynamics with artificial intelligence”. In: Nature Reviews Neuroscience 27.2 (2026), pp. 87–100.

[9] Ali Essam Ghareeb, Benjamin Chang, Ludovico Mitchener, Angela Yiu, Caralyn J. Szostkiewicz, et al. “A multi-agent system for automating scientific discovery”. In: Nature 655.8122 (2026), pp. 497–505. [10] Chris Lu, Cong Lu, Robert Tjarko Lange, Yutaro Yamada, Shengran Hu, et al. “Towards end-toend automation of AI research”. In: Nature 651.8107 (2026), pp. 914–919. [11] Juraj Gottweis, Wei-Hung Weng, Alexander Daryin, Tao Tu, Petar Sirkovic, et al. “Accelerating scientific discovery with Co-Scientist”. In: Nature 655.8122 (2026), pp. 487–496. [12] Samuel Alber, Bowen Chen, Eric Sun, Alina Isakova, Aaron J. Wilk, et al. “Cellvoyager: Ai compbio agent generates new insights by autonomously analyzing biological data”. In: Nature Methods 23.4 (2026), pp. 749–759. [13] Eser Aygün, Anastasiya Belyaeva, Gheorghe Comanici, Marc Coram, Hao Cui, et al. “An AI system to help scientists write expert-level empirical software”. In: Nature 654.8120 (2026), pp. 909–916. [14] Shiyang Feng, Runmin Ma, Xiangchao Yan, Yue Fan, Yusong Hu, et al. “InternAgent-1.5: A Unified Agentic Framework for Long-Horizon Autonomous Scientific Discovery”. In: arXiv preprint arXiv:2602.08990 (2026). [15] InternAgent Team, Bo Zhang, Shiyang Feng, Xiangchao Yan, Jiakang Yuan, et al. “InternAgent: When Agent Becomes the Scientist–Building Closed-Loop System from Hypothesis to Verification”. In: arXiv preprint arXiv:2505.16938 (2025). [16] Lei Bai, Zongsheng Cao, Yang Chen, Zhiyao Cui, Shangheng Du, et al. Scaling the Horizon, Not the Parameters: Reaching Trillion-Parameter Performance with a 35B Agent. 2026. arXiv: 2606.30616 [cs.CL]. url: https://arxiv.org/abs/2606.30616.

17

[17] Thomas Marwitz, Alexander Colsmann, Ben Breitung, Christoph Brabec, Christoph Kirchlechner, et al. “Predicting new research directions in materials science using large language models and concept graphs”. In: Nature Machine Intelligence 8.4 (2026), pp. 535–544. [18] Xiaoliang Luo, Akilles Rechardt, Guangzhi Sun, Kevin K. Nejad, Felipe Yáñez, et al. “Large language models surpass human experts in predicting neuroscience results”. In: Nature Human Behaviour 9.2 (2025), pp. 305–315. issn: 2397-3374. [19] F Alexander Wolf, Philipp Angerer, and Fabian J Theis. “SCANPY: large-scale single-cell gene expression data analysis”. In: Genome Biology 19.1 (2018), p. 15. [20] Kexin Huang, Serena Zhang, Hanchen Wang, Yuanhao Qu, Yingzhou Lu, et al. “Autonomous biomedical research with an artificial intelligence agent”. In: Science (2026), eadz4351. doi: 10.1126/science.adz4351. [21] Keyan Ding, Jian Yu, Jiaxin Huang, Yuxin Yang, Qi Zhang, et al. “SciToolAgent: a knowledgegraph-driven scientific agent for multitool integration”. In: Nature Computational Science 5 (2025), pp. 962–972. [22] Zhizheng Wang, Qiao Jin, Chih-Hsuan Wei, Shubo Tian, Po-Ting Lai, et al. “GeneAgent: selfverification language agent for gene-set analysis using domain databases”. In: Nature Methods 22 (2025), pp. 1677–1685. [23] Xiangru Tang, Qiao Jin, Kunlun Zhu, Tongxin Yuan, Yichi Zhang, et al. “Risks of AI scientists: prioritizing safeguarding over autonomy”. In: Nature Communications 16 (2025), p. 8317. [24] Lisa Messeri and M. J. Crockett. “Artificial intelligence and illusions of understanding in scientific research”. In: Nature 627 (2024), pp. 49–58. [25] Anthropic. Claude Opus 4.7 System Card. Tech. rep. Accessed: 2026-06-29. Anthropic, 2026. url: https://www.anthropic.com/claude-opus-4-7-system-card. [26] OpenAI. GPT-5.5 System Card. Tech. rep. Accessed: 2026-06-29. OpenAI, 2026. url: https: //openai.com/index/gpt-5-5-system-card/. [27] Duncan S. Mearns, Joseph C. Donovan, António M. Fernandes, Julia L. Semmelhack, and Herwig Baier. “Deconstructing Hunting Behavior Reveals a Tightly Coupled Stimulus-Response Loop”. In: Current Biology 30.1 (2020), 54–69.e9. [28] Steffen Schneider, Jin Hwa Lee, and Mackenzie Weygandt Mathis. “Learnable latent embeddings for joint behavioural and neural analysis”. In: Nature 617.7960 (2023), pp. 360–368. [29] Charles Findling, Félix Hubert, International Brain Laboratory, Luigi Acerbi, Brandon Benson, et al. “Brain-wide representations of prior information in mouse decision-making”. In: Nature 645.8079 (2025), pp. 192–200. [30] Center for AI Safety, Scale AI, and HLE Contributors Consortium. “A benchmark of expert-level academic questions to assess AI capabilities”. In: Nature 649.8099 (2026), pp. 1139–1146. [31] Carlos E. Jimenez, John Yang, Alexander Wettig, Shunyu Yao, Kexin Pei, et al. “SWE-bench: Can Language Models Resolve Real-World GitHub Issues?” In: arXiv preprint arXiv:2310.06770 (2023). [32] Jun Shern Chan, Neil Chowdhury, Oliver Jaffe, James Aung, Dane Sherburn, et al. “MLE-bench: Evaluating Machine Learning Agents on Machine Learning Engineering”. In: International Conference on Learning Representations. 2025. [33] Carsen Stringer, Marius Pachitariu, Nicholas Steinmetz, Charu Bai Reddy, Matteo Carandini, et al. “Spontaneous behaviors drive multidimensional, brainwide activity”. In: Science 364.6437 (2019), eaav7893. 18

[34] Simon Musall, Matthew T. Kaufman, Ashley L. Juavinett, Steven Gluf, and Anne K. Churchland. “Single-trial neural dynamics are dominated by richly varied movements”. In: Nature Neuroscience 22.10 (2019), pp. 1677–1686. [35] Nicholas A. Steinmetz, Peter Zatka-Haas, Matteo Carandini, and Kenneth D. Harris. “Distributed coding of choice, action and engagement across the mouse brain”. In: Nature 576.7786 (2019), pp. 266–273. [36] Oliver Rübel, Andrew Tritt, Ryan Ly, et al. “The Neurodata Without Borders ecosystem for neurophysiological data science”. In: eLife 11 (2022), e78362. [37] Google. Gemini 3.1 Pro Model Card. Tech. rep. Accessed: 2026-06-29. Anthropic, 2026. url: https://storage.googleapis.com/deepmind-media/Model-Cards/Gemini-31-Pro-Model-Card.pdf. [38] Shunyu Yao, Jeffrey Zhao, Dian Yu, Nan Du, Izhak Shafran, et al. “ReAct: Synergizing Reasoning and Acting in Language Models”. In: International Conference on Learning Representations. 2023. [39] Adam A. Atanas, Jungsoo Kim, Ziyu Wang, et al. “Brain-wide representations of behavior spanning multiple timescales and states in C. elegans”. In: Cell 186.19 (2023), 4134–4151.e31. [40] Ryosuke Tanaka and Ruben Portugues. “Plastic landmark anchoring in zebrafish compass neurons”. In: Nature 650.8102 (2026), pp. 673–680. [41] Melanie A. Basnak, Anna Kutschireiter, Tatsuo S. Okubo, Albert Chen, Pavel Gorelik, et al. “Multimodal cue integration and learning in a neural representation of head direction”. In: Nature Neuroscience 28.8 (2025), pp. 1729–1740. [42] Yuan Gao, Cindy T. J. van Velthoven, Changkyu Lee, et al. “Continuous cell-type diversification in mouse visual cortex development”. In: Nature 647.8088 (2025), pp. 127–142. [43] Mikhail Genkin, Krishna V. Shenoy, Chandramouli Chandrasekaran, and Tatiana A. Engel. “The dynamics and geometry of choice in the premotor cortex”. In: Nature 645.8079 (2025), pp. 168–176. [44] Bin Yu, Yuxuan Yue, Chi Ren, Rui Yun, Byungkook Lim, et al. “Cholinergic feedback for modalityand context-specific modulation of sensory representations”. In: Science 388.6753 (2025), pp. 1324–1329. [45] Antoine Légaré, Mado Lemieux, Vincent Boily, Sandrine Poulin, Arthur Légaré, et al. “Structural and genetic determinants of zebrafish functional brain networks”. In: Science Advances 11.28 (2025), eadv7576. [46] Emmanuel Marquez-Legorreta, Greg M. Fleishman, Luuk W. Hesselink, et al. “Whole-Brain Co-Mapping of Gene Expression and Neuronal Activity at Cellular Resolution in Behaving Zebrafish”. In: bioRxiv (2026). [47] Anthropic. Claude Opus 4.6 System Card. Tech. rep. Accessed: 2026-06-29. Anthropic, 2026. url: https://www.anthropic.com/claude-opus-4-6-system-card.

19

A. Supplementary Information A.1. Supplementary Notes A.1.1. Implementation details of SeekBrain This Supplementary Note provides the implementation details of analysis execution engine. For each task, SeekBrain receives a natural-language neuroscience question together with optional dataset metadata and data-path information. The system then (1) retrieves a task-relevant analysis recipe from the Neuroscience Analysis Repertoire, (2) formulates a data-grounded analysis schema, (3) generates and executes Python code, (4) validates the output, and (5) interprets the validated figures, statistics, and execution traces to produce a self-contained scientific conclusion. In the configuration used for all reported experiments, analytical formulation, code execution, validation, and result interpretation are routed to Claude Opus 4.7. The architecture supports modular replacement of base LLMs for each agent. Analysis execution workflow. Before analytical formulation, the workflow retrieves a targeted recipe from the Neuroscience Analysis Repertoire. If no verified recipe is matched, the workflow initiates a task-specific recipe crystallization process: it identifies relevant publications with associated code repositories, selects task-relevant paper and repository context, and crystallizes the paper-code correspondence into a reusable analysis recipe. The retrieved or newly generated recipe is then supplied to the Formulation Agent as a methodological reference. The Formulation Agent produces a data-grounded methodological schema by integrating the user question, dataset metadata, recipe, and literature background. The Code Agent writes, executes, and debugs Python analysis code within a bounded working directory, producing required code, outputting the generated scripts, execution logs, figures, statistics, and execution logs. The workflow explicitly prohibits the use of dummy, placeholder, or proxy data. Instead, missing data errors are either resolved via user-supplied paths or formally recorded as limitations. The Validation Agent evaluates whether the output answers the original question and conforms to the schema, drawing on the figure, statistics, optional reference figure, code-agent confidence tag, and literature-grounded background. The Validation Agent can reject a result by flagging specific mismatches. This feedback is routed back to the Code Agent for a maximum of three validation cycles before final interpretation. Finally, the Interpretation Agent reports scientific conclusions extracted from the result figure, quantitative statistics, execution log, methodological schema, plausibility warnings, and literature background. Analysis recipe grounding. Within the analysis workflow, the Neuroscience Analysis Repertoire functions as procedural memory. Any retrieved recipe serves as a methodological reference rather than a strict directive. The Formulation Agent evaluates its alignment with the current question, data modality, and validation requirements before deciding whether to adopt, adapt, or discard it. This approach ensures that the user’s question and dataset metadata remain the primary constraints, while literature-validated procedures provide adaptable guidance.

20

A.2. Supplementary Figures

Spatial response mapping

Ring attractor phase mapping

Structural network visualization

Brain-wide modulation mapping

Temporal offset decoding

Neural spiking dynamics analysis

Neural embedding visualization

Head direction measure

Functional connectivity estimation

Supplementary Figure 1 | Examples of SeekBrain’s analysis results on BrainArena.

21

Task topic: Structural network visualization

Evaluation Summary: Successfully reproduced bilateral SC network with proper anatomical alignment and appropriate structural interpretation.

Evaluation Summary: Reproduced network visualization but with side-view misalignment and missing bottom-quartile thresholding; failed to interpret edge weights as structural wiring estimates.

Task topic: Neural prior2 decoding Supplementary Figure | Cases with and without established analysis recipes.

Prompt for Benchmark Task Execution Your task is to complete a scientific data analysis task based on the provided specifications. Task Parameters: • Question file: {QUESTION_FILE} • Data directory: {DATA_DIR} Output Requirements: Evaluation Summary: Successfully Bayesian prior Save all required outputs to thereproduced {OUTPUT_DIR} directory: Evaluation Summary: Correctly identified vlOFC prior decoding

tracking in vlOFC with proper decoding methodology and cautious during ITI but failed to distinguish uncorrected from corrected R² and interpretation distinguishing graded belief states from discrete block omitted the movement-minimization rationale for ITI analysis. 1. conclusion.md: Write 1 to 3 scientific conclusions. Each conclusion must be a single, self-contained sentence encoding.

or short paragraph focusing on key scientific findings, supported by quantitative evidence where available. 2. analyze.md: Document the complete analysis process, including data exploration, methodology, result interpretation, and your reasoning trace. 3. analyze.png: Generate a figure that closely matches the visual and analytical specifications in the question file. 4. analyze.py: Provide the complete executable analysis code script.

Execution Strategy: • First, explore the data structure. Next, write and execute the analysis code. Finally, draft analyze.md and conclusion.md. • Prioritize executing an end-to-end pipeline over exhaustive data exploration. • Ensure the figure, analysis document, and conclusions are generated before the timeout, even if the results are partial or imperfect. • Begin execution immediately without asking for user confirmation.

Supplementary Figure 3 | Prompt for benchmark task execution. Each coding agent receives this prompt alongside a task-specific question file and data directory, instructing it to perform data exploration, write analysis code, generate a figure, and synthesize scientific conclusions within a strict timeout

22

Prompt for LLM-as-Judge Evaluation You are a strict but fair evaluator for a scientific benchmark. Your task is to evaluate the submitted neuroscience benchmark results against the provided rubric. Grade only the submitted artifacts. Do not credit claims or figure content that are absent from the submission. Do not penalize phrasing differences if the scientific meaning is correct. Inputs Provided: 1. Ground-truth figure from the original paper (Image) 2. Model’s reproduced figure (Image) 3. Rubric (JSON): Ground-truth conclusion, details scoring criteria, point allocations (max_point), and deduction rules (deduction_rule) 4. Model’s reproduced code (Python): The executable analysis and visualization code used to generate the reproduced figure 5. Model’s submission (Text): The conclusion or methods summary generated by the model Grading Rules: • The maximum total score is {max_score}. Assign each rubric item a score from 0 up to its specified max_point. • Award partial credit if the submission partially satisfies a criterion. • Apply each deduction_rule strictly when applicable. • Treat the rubric’s ground-truth conclusion as the scoring reference, not as text automatically credited to the submission. • Scores may be integers or decimals. Task & Output Format: Score the model’s conclusion from 0 to 100 using the rubric. Provide your step-by-step reasoning and the final score. Return ONLY valid JSON with no Markdown fences or additional text, adhering strictly to this schema: { "reasoning": "score":

"<your detailed step-by-step reasoning>",

<final numerical score>

}

Supplementary Figure 4 | Prompt for LLM-as-Judge evaluation. The evaluator receives multimodal inputs: reproduced codes, reproduced figures, structured analysis metadata, a detailed rubric, and the model’s written conclusions. It then assigns a score from 0 to 100 following strict grading rules with partial credit and specified deduction policies. Example of a Neuroscience Analysis Recipe Recipe File Tree:

. |– RECIPE.md |– metadata.json |– references | |– implementation_notes.md | ‘– scientific_context.md ‘– scripts |– run_analysis.py ‘– utils.py Key Contents of RECIPE.md: • Objective: Analyze how head-direction cell tuning transforms when a zebrafish experiences a symmetric virtual scene. Align tuning curves to each neuron’s pre-learning preferred orientation and compare pre-learning, learning,

23

and post-learning epochs to identify systematic bimodal or double-peaked responses caused by visual ambiguity. • Input: Synchronized behavioral data containing scene orientation 𝜃 and epoch boundaries; neural data containing ROI-matched Δ𝐹 / 𝐹 traces; stable ROI_IDs across all three epochs; and MATLAB circular-statistics and plotting utilities, including Kent_func, stats_plot, and im_dissim. • Processing: Identify head-direction cells in the pre-learning epoch using a Rayleigh significance criterion ( 𝑝 < 0.05), estimate each cell’s preferred orientation 𝑏, and center its tuning curve with 𝜃 − 𝑏 wrapped to [−𝜋, 𝜋]. Sort ROIs once by pre-learning 𝑏 and preserve this order across all epochs. Construct three tuning-curve heatmaps with a shared color scale, then fit a two-component mixture model to learning-epoch curves to quantify the secondary peak near ±𝜋. • Output: (1) A synchronized three-panel heatmap for the Pre, Learning, and Post epochs; (2) centered_tuning_matrices.mat; (3) preferred_orientations.mat; (4) ROI_order.mat; and (5) mixture-model statistics describing bimodality. A valid result preserves identical ROI ordering and color limits across panels and reveals whether symmetric-scene learning produces a population-level second tuning peak.

Supplementary Figure 5 | Example of a neuroscience analysis recipe. The complete skill directory tree is shown first, followed by the essential objective, input specification, and output specification.

A.3. Supplementary Tables Supplementary Table 1 | Neuroscience analysis recipes. Name

Description

ana-mouse-behavior-markovepoch_interval

Decomposes trial-by-trial choices into interpretable laBehavioral Modality tent behavioral strategies using a logistic-emission GLMHMM fitted by expectation-maximization and decoded with Viterbi inference. It produces state-specific GLM weights, posterior state-probability traces, decoded sequences, and AIC/BIC model-selection scores. Fits Bayes-optimal, stimulus-kernel, and action-kernel mod- Behavioral Modality els to trial-by-trial choice data, estimates session-level log evidence, and applies random-effects Bayesian model selection. It outputs model frequencies, exceedance probabilities, evidence matrices, and summary comparisons of competing behavioral strategies. Behavioral Modality Analyzes trial-aligned behavioral and calcium-imaging data using OASIS deconvolution, learning-curve fitting, and groupwise statistical tests. It produces performance time courses and significance estimates for CNO-induced olfactory learning impairment and post-washout recovery. Pools binary response masks to calculate percentages of Dynamics Modality excited, suppressed, and non-responsive granule cell–odor pairs under task and passive conditions. It outputs summary statistics and comparative side-by-side pie charts to reveal sparse, state-dependent odor-response distributions. Detects statistically significant brain-wide co-fluctuation Dynamics Modality events in zebrafish calcium data using edge-time-series RMS, region-wise circular-shift null distributions, and Benjamini–Hochberg FDR. It outputs event-level regional activity vectors for downstream brain-state clustering, plus RMS and null statistics. Tests cross-day zebrafish larva re-identification by correlat- Dynamics Modality ing vectorized upper triangles of functional-connectivity matrices and symmetrizing pairwise similarities. It outputs identification accuracy, a larva-by-larva similarity matrix, and a heatmap marking each larva’s strongest match.

ana-mouse-behavior-glmperformance

ana-mouse-behavior-anovatask_performance

multi-mouse-caodor_tuning_ratio_chi_squarepie_plot

ana-zebrafish-ca-circular_nullcofluctuation_event_trace

ana-zebrafish-cacross_day_correlationidentity_fingerprint_matrix

Modality

24

Name

Description

Modality

ana-zebrafish-caTests zebrafish functional-connectivity fingerprint signifiDynamics Modality fingerprint_mannwhitney-box_plot cance by comparing cross-day intra- and inter-individual similarities with a two-sided Mann–Whitney U test. It outputs the similarity matrix and group vectors plus a jittered min-to-max boxplot annotated with the p-value. Dynamics Modality ana-zebrafish-ca-pearsonBuilds per-larva and group-average zebrafish functional fc_matrix connectivity matrices by Gaussian-smoothing regional Δ𝐹 / 𝐹0 traces and computing pairwise Pearson correlations after quality-control exclusions. It outputs region-by-region matrices quantifying brain-wide functional coupling and synchronization. Dynamics Modality ana-zebrafish-ca_fc_matrixQuantifies inter-larva functional-connectivity stereotypy cross_individual_correlationby correlating vectorized upper triangles of individual network_stereotypy_matrix zebrafish FC matrices. It outputs an ordered larva-by-larva Pearson similarity matrix, representative-larva ranking, summary statistics, and lower-triangular heatmap. pre-zebrafish-caBuilds a cross-fish bilateral zebrafish brain-state matrix by Dynamics Modality detrend_minfilter_zscoresymmetrically excluding unusable regions, applying minactivity_state_matrix filter Δ𝐹 / 𝐹0 detrending, and z-scoring regional calcium traces before time-wise concatenation. It produces standardized cohort activity arrays for heatmaps, co-fluctuation detection, and brain-state clustering. Dynamics Modality pre-zebrafish-caConverts motion-corrected zebrafish fluorescence to baseline_correlation_minfilterper-neuron Δ𝐹 / 𝐹0 using a Gaussian-smoothed runningdff_region_timeseries minimum baseline, then averages activity across atlas regions. It saves neuronal and regional time series and produces an anterior-to-posterior stacked regional-trace plot with the stimulus window marked. ana-zebrafish-ca-tsne_hdbscanClusters high-amplitude zebrafish whole-brain coactivation Dynamics Modality correlation_matrix_heatmap events using t-SNE and HDBSCAN, reassigning noise by correlation to native-space cluster centroids. It outputs discrete brain-state labels, embeddings, and a cluster-sorted event correlation heatmap revealing block-diagonal state structure. vis-zebrafish-fc_atlasProjects a group-averaged Pearson functional connectivDynamics Modality edge_projection-network_map ity matrix onto Mapzebrain top and side anatomy using bilateral regional centroids as nodes and correlationweighted edges, hiding the weakest quartile. It produces a publication-quality spatial network figure for assessing the anatomical distribution of strong connections. Dynamics Modality ana-zebrafish-ca-eigen_spectrum- Visualizes the diffusion-map eigenvalue spectrum of zefc_scree_plot brafish whole-brain functional connectivity as a scree plot. It reveals spectral decay and the elbow supporting retention of leading functional gradients, including the 50-mode cutoff used for structural–functional gradient matching. ana-zebrafish-ca-tsne_hdbscanEmbeds high-amplitude zebrafish whole-brain coactivaDynamics Modality cluster_scatter_plot tion events with t-SNE, clusters them using HDBSCAN, and reassigns noise by correlation to centroids in the original regional space. It outputs cluster labels, embedding coordinates, and a brain-state-coloured 2-D scatter plot. ana-mouse-omicsQuantifies developmental changes in scRNA-seq cell-type Molecular Modality hierarchical_clusteringdiversity using cross-age integration, resolution-consistent type_temporal_trajectory clustering, downsampling, and cluster stability analysis. It produces temporal plots and tables of cluster counts to identify emerging populations and diversification plateaus.

25

Name

Description

Modality

ana-mouse-omicsdifferential_genetemporal_expression_profile

Quantifies subclass-specific developmental gene expresMolecular Modality sion changes in mouse visual cortex using replicate-aware pseudobulk aggregation and DESeq2 or MAST. It produces per-gene log2 fold changes and adjusted significance values, DEG summaries, volcano plots, and up/downregulated gene counts across subclasses. ana-mouse-omics-gamModels lineage-specific marker expression across develop- Molecular Modality expression_profile mental time using generalized additive models and spline derivatives, classifying trajectories as Up, Down, Transient, or Constant. It outputs a classification table and faceted expression curves with confidence intervals for excitatory, inhibitory, and glial lineages. ana-mouse-omicsAnalyzes developmental scRNA-seq using integrated hierar- Molecular Modality hierarchical_clusteringchical clustering, temporal linkage, and Jaccard similarity scatter_type_temporal_trajectory to track cell-type emergence, branching, and stability. It outputs age-resolved subcluster counts, stability metrics, and alluvial diagrams of lineage trajectories. Molecular Modality vis-zebrafish-gene-zscoreVisualizes larval zebrafish spatial transcriptomics as a expression_matrix_heatmap grayscale heatmap of genes across retained brain regions, z-scoring genes within each region. It produces a figure and descriptive statistics to inspect regional marker patterns before gene–connectivity comparisons. Structural Modality ana-zebrafish-swc-nblastClassifies zebrafish cerebellar efferent neurons by commorphological_metric puting all-by-all NBLAST morphological similarities and applying Ward hierarchical clustering. It produces an eightcluster dendrogram, neuron classification table, and silhouette scores distinguishing Purkinje and eurydendroid morphologies. ana-zebrafishQuantifies distance-dependent connectivity in the larval Structural Modality connectome-pearsonzebrafish brain using Euclidean centroid distances, volumecorrcoef_distance_scatter_plot normalized neurite lengths, and Pearson correlation. It also applies hierarchical consensus clustering to produce statistical evidence, a distance–strength plot, and a dendrogram of stable brain communities. ana-zebrafish-swc-nblastComputes all-by-all NBLAST similarities for registered SWC Structural Modality morphology neuronal skeletons, applies hierarchical clustering to define morphological groups, and renders them in a standard 3D brain template. Outputs a dendrogram, neuron-to-cluster map, and cluster visualization gallery. Structural Modality ana-zebrafish-sc_matrixIdentifies zebrafish structural-connectivity modules by hierarchical_consensus_clustering- applying Ward hierarchical clustering to co-assignment consensus_matrix_dendrogram probabilities from multi-resolution Louvain partitions. It outputs module labels reordered by anterior-posterior position, the co-assignment matrix, and a dendrogram-matrix visualization. vis-zebrafish-sc_atlasVisualizes a bilateral zebrafish structural connectivity net- Structural Modality edge_projection-network_map work by overlaying atlas-centroid nodes and thresholded, weight-coded edges on dorsal and sagittal MapzebrainAtlas projections. It outputs a dual-panel figure and JSON statistics summarizing regions, edges, weights, and coordinate ranges. ana-zebrafishApplies unscaled PCA to structural predictors derived Structural Modality structure_connectivityfrom directed and undirected zebrafish connectomes and graph_predictor_pca-scree_plot inter-regional Euclidean distances. It reports and plots percomponent and cumulative explained variance for the first 10 PCs to assess predictor redundancy and dimensionality.

26

Name

Description

Modality

pre-zebrafishcell_coordinate_anatomyregion_censusatlas_contour_map_cell_density

Censuses neurons per larva across Mapzebrain regions and Structural Modality excludes retina plus regions with fewer than 25 neurons in at least half of larvae. It reports whole-brain counts, writes the downstream QC mask, and generates anatomical projection backgrounds. ana-zebrafish-sc_matrix_anatomy- Projects consensus-Louvain structural-connectivity module Structural Modality hierarchical_consensus_clustering- labels onto bilateral zebrafish atlas centroids in color-coded module_topography top and side views. It produces an anatomical topography figure and module counts to show anterior–posterior organization and spatial contiguity. ana-zebrafish-sc_matrixPlots the larval zebrafish structural-connectivity diffusion- Structural Modality eigen_spectrum-sc_scree_plot map eigenvalue spectrum as a scree curve to reveal its decay and elbow. It outputs a figure and summary statistics for the eigenvalue count, first value, mean, median, minimum, and maximum. vis-zebrafish-sc-module_sortedVisualizes the zebrafish structural connectome by sorting Structural Modality sc_block_heatmap its adjacency matrix by four consensus-Louvain modules, squaring weights for contrast, and outlining diagonal blocks. It produces a heatmap and descriptive statistics showing within- versus between-module connectivity structure. multi-zebrafish-ca_behaviorQuantifies zebrafish neural bump–scene phase coupling Dynamics-Behavior vonmises-phase_locking during jump and noise epochs by fitting von Mises concen- Coupling tration to angular offsets and comparing observed 𝜅 values with circular-shift null distributions. It outputs per-fish coupling statistics, 95th-percentile significance thresholds, and summary plots. multi-zebrafish-ca_behaviorQuantifies zebrafish scene–neural bump alignment disrup- Dynamics-Behavior circular_statistics-phase_locking tion using circular angular differences and thresholding Coupling to calculate each fish’s out-of-phase fraction across prelearning, learning, and post-learning epochs. It outputs per-subject fractions and a paired longitudinal trajectory plot. multi-drosophila-ca_behaviorQuantifies Drosophila E-PG head-direction encoding by Dynamics-Behavior circular_statisticspreprocessing fluorescence data, estimating activity-bump Coupling encoding_accuracy phase, and calculating vector strength relative to heading. It statistically compares visual cue conditions and produces an SVG grouped scatter plot of per-fly accuracy, means, error bars, and significance. multi-zebrafish-ca_behaviorQuantifies neural bump stability relative to visual scenes by Dynamics-Behavior vonmises_shufflingfitting von Mises distributions to circular phase offsets and Coupling distribution_bar_plot comparing observed concentration parameters with shuffled null distributions. It outputs epoch-level stability statistics and visualizations identifying significant phase-locking during perturbations. multi-zebrafish-ca_behaviorAnalyzes circular correlations between decoded neural Dynamics-Behavior circular_correlation-tuning_curve bump phase and scene orientation across pre-learning, Coupling learning, and post-training epochs. It outputs epoch-specific slopes and scatter plots to reveal shifts between 1:1 and 2:1 neural mappings. multi-c_elegans-ca_behaviorProjects pooled multi-animal neural activity into a 2D Dynamics-Behavior umap-embedding_scatter_plot UMAP and overlays CePNEM-derived feeding-tuning inCoupling dices. It produces a neural-manifold plot and per-neuron coordinate table revealing graded transitions from feedinginhibited to feeding-activated encoding.

27

Name

Description

Modality

multi-mouse-spike_behavior-glm- Decodes trial-by-trial prior probability from binned ORBvl Dynamics-Behavior decoding_metrics inter-trial population activity using cross-validated linear Coupling models benchmarked against a Bayesian observer. It outputs decoded and reference prior traces plus session-level Pearson correlation and noise-ceiling-corrected R2. multi-c_elegans-ca_behaviorAnalyzes neuronal velocity encoding by correlating binned Dynamics-Behavior correlation-timeseries_trace firing rates with causally EWMA-filtered velocity across Coupling decay constants. It outputs per-neuron correlation profiles, optimal integration timescales, significance against shuffled spike-time nulls, and representative plots. multi-c_elegans-ca_behaviorFits Ridge regression encoding kernels and exponential Dynamics-Behavior half_decay-distribution decays to estimate neuronal half-decay times, retaining Coupling neurons significant against a shuffle-based null. It outputs an ECDF of population timescales with ±1 SD inter-animal shading and per-neuron summary data. multi-c_elegans-ca_behavior-ridge- Uses time-lagged ridge regression to decode past, current, Dynamics-Behavior timeseries_trace and future behavioral velocity from neural activity, comCoupling paring smoothed prediction errors against shuffled-label controls. It outputs performance across temporal offsets and plots mean decoding curves with variability across animals. multi-zebrafish-ca_behaviorAnalyzes head-direction cell tuning across pre-learning, Dynamics-Behavior circular_statistics-tuning_curve learning, and post-learning epochs using circular fitting, Coupling preferred-orientation alignment, and Gaussian-mixture bimodality assessment. It outputs synchronized, consistently sorted heatmaps and tuning-shift statistics revealing symmetry-induced double-peaked responses. multi-macaque-spike_behaviorQuantifies prior-dependent shifts in superior colliculus Dynamics-Behavior glm-performance firing by epoching synchronized post-stimulus spikes and Coupling fitting logistic regressions with contrast, prior, and neural activity. It produces regression coefficients, prior-stratified neurometric curves and PSE shifts, and stimulus-aligned PSTHs. multi-mouse-spike_behaviorFits exponential decays to neural latent-state autocorreDynamics-Behavior decay_correlation-scatter_plot lations and behavioral dwell times, then correlates their Coupling inverse time constants across quality-controlled, synchronized sessions. It produces a session-level scatter plot and a table of the correlation coefficient, p-value, and sample size. multi-c_elegans-ca_behavior-pca- Performs PCA on z-scored neural activity and Ward hierDynamics-Behavior pc_feature_heatmap archical clustering of neurons using the top five PC loadCoupling ings. It aligns clusters with CePNEM behavioral tuning weights and produces a dendrogram and heatmaps revealing behavior-associated functional ensembles. multi-macaque-spike_behaviorGenerates stimulus-aligned PSTHs for 14 PMd neurons Dynamics-Behavior classified_firing_rate-psth_trace using a 75 ms sliding window and 10 ms step, grouped Coupling by choice and difficulty. It outputs a multi-panel figure with trial-averaged firing rates and SEM shading to reveal condition-dependent temporal heterogeneity. multi-drosophila-ca_behaviorAnalyzes movement-gated Drosophila EPG imaging and Dynamics-Behavior pearson-scatter_plot behavior data using circular phase decoding, von Mises Coupling bump-width fitting, and Pearson correlation. It outputs perfly accuracy and width metrics, correlation statistics, and a regression scatter plot of heading-encoding accuracy versus bump width.

28

Name

Description

multi-macaque-spike_behaviorclassified_firing_ratepsth_distribution

Computes stimulus-aligned, trial-averaged single-neuron Dynamics-Behavior PMd PSTHs by choice and difficulty using a 75-ms sliding Coupling boxcar, 10-ms steps, reaction-time masking, and SEM. It produces condition-wise plots and temporal-profile classifications to characterize ramping, transient, non-monotonic, and mixed decision-related activity. Applies CP tensor decomposition to PMd spike-rate data Dynamics-Behavior organized by neuron, time, and trial, then relates latent Coupling trial factors to choice and stimulus coherence. It outputs shared ramping or transient temporal motifs, behavioral sensitivities, and reconstructed population activity. Applies canonical polyadic decomposition to trial-binned Dynamics-Behavior PMd firing-rate tensors, extracting temporal factors and Coupling trial loadings linked to stimulus coherence and choice. It produces latent PSTH visualizations and condition-sorted loading plots validated against trial-averaged neural activity. Computes event-aligned population PSTHs by binning and Dynamics-Behavior Gaussian-smoothing spike trains, then averaging and option- Coupling ally baseline-normalizing activity by stimulus difficulty and behavioral choice. It outputs condition-segregated firingrate arrays and plots with SEM to reveal decision-related temporal dynamics. Estimates session-specific trial-history decay constants from Dynamics-Behavior behavioral logistic kernels and neural GLMs, then measures Coupling their cross-session alignment using modality-stratified Pearson correlation. Outputs paired decay estimates, correlation statistics, and neural-versus-behavioral scatter plots. Aggregates functionally classified neuron centroids across Dynamics-Behavior zebrafish in a common 3D atlas and uses permutationCoupling based spatial p-values to identify locations consistent across individuals. It outputs filtered centroid and p-value arrays plus density-colored whole-brain maps for motor and visual response populations. Creates a single-larva sensorimotor figure aligning binaDynamics-Behavior rized swimming, light/dark stimuli, and calcium traces Coupling for the 50 most extreme neurons in four precomputed correlation-ranked ensembles. It overlays individual traces and ensemble means to qualitatively show motor- and stimulus-linked activity without population-level inference. Integrates atlas registration and NBLAST clustering with Dynamics-Behaviorstimulus-aligned calcium analysis and neural–swimming Anatomy-Gene correlations to characterize glyt2–gfra1a zebrafish neurons. Coupling It produces 3D brain maps, peri-stimulus calcium traces, and comparative violin plots. Compares baseline-corrected neurometric shifts from elec- Dynamics-Behaviortrophysiology and widefield imaging across matched corStimulus Coupling tical regions using Spearman rank correlation. It produces an annotated scatter plot with a regression line, 95% confidence interval, correlation coefficient, and p-value. Benchmarks latent dynamical models against trial-averaged Dynamics-BehaviorPSTH baselines using held-out marginal log-likelihood Stimulus Coupling differences across neurons. It outputs the improvement distribution, mean population gain, and a Wilcoxon p-value testing whether gains exceed zero.

multi-macaque-spike_behaviortensor_decompositiontemporal_motif

multi-macaque-spike_behaviorclassified_firing_rate-raster_psth

multi-mouse-spike_behaviorclassified_firing_ratepsth_heatmap

multi-mouse-spike_behaviorexponential_kernel_glmdecay_distribution

multi-zebrafish-ca_behaviorpermutation-spatial_density

multi-zebrafish-ca_behaviorlinear_regressionsensorimotor_trace

multi-zebrafishca_behavior_anatomy_genecorrelation_visual_responseneuronal_subtype_features multi-mousespike_behavior_visual_stimspearman-scatter_plot

multi-macaquespike_behavior_visual_stimlatent_dynamical_modellog_likelihood_ratio

Modality

29

Name

Description

Modality

multi-zebrafish-gene_anatomyexpression_level_clusteringheatmap

Generates a publication-quality heatmap of larval zebrafish Molecule-Structure gene expression after excluding selected brain regions and Coupling z-scoring expression across genes within each region. Outputs the heatmap and summary statistics to reveal regionspecific spatial expression patterns relevant to functional connectivity. multi-zebrafish-gene_anatomyMaps eomesa- and pvalb7-expressing zebrafish neurons by Molecule-Structure expression_count_mappingdeformably registering EASI-FISH volumes to MapZebrain, Coupling spatial_density assigning RNA spots to segmented cells, and thresholding expression. It produces atlas-aligned 3D point clouds and per-cell tables for individual and co-expressing populations. multi-mouse-spike_visual_stimAnalyzes Allen Neuropixels data using responsiveness Stimulus-Dynamics svm-performance tests and cross-validated linear SVM decoding with unitCoupling count subsampling across brain regions. It outputs regional responsive-neuron fractions, decoding accuracy statistics, null comparisons, and comparative visualizations. multi-mouse-spike_anatomyAnalyzes Neuropixels firing-rate differences between left Structure Function significant_neuron_fractionand right stimuli using Mann–Whitney U tests with BenCoupling atlas_contour_map jamini–Hochberg FDR correction, aggregated by Beryl brain region. It produces a right-hemisphere Swanson flatmap showing each region’s fraction of significantly modulated neurons. multi-zebrafish-sc_fc_matrixClassifies zebrafish region pairs as unconnected, unidirec- Structure Function directionality_anovational, or bidirectional from the directed structural connec- Coupling matrix_histogram tome, pools per-larva functional connectivity by class, and compares distributions using one-way ANOVA. It outputs category matrices, FC arrays, density plots, and F and P statistics. multi-zebrafishComputes structural and functional connectivity gradients Structure Function sc_fc_matrix-diffusion_mapin the larval zebrafish brain using diffusion map embedding, Coupling connectivity_gradient then correlates the first two gradients across modalities. It saves gradient and eigenvalue arrays and produces bilateral atlas overlays showing their anatomical organization. multi-zebrafish-sc_fc_matrixComputes approximately graph-theoretic predictors from ze- Structure Function graph_predictorsbrafish structural connectivity and ranks them by edgewise Coupling fc_variance_explained squared Pearson correlation with functional connectivity. It reports variance explained relative to undirected structural connectivity, identifying matching-index measures as the strongest predictors. multi-zebrafish-sc_fc_matrixMatches the zebrafish structural and functional connecStructure Function hungarian_matchingtivity gradients using all-pairs Pearson correlations and Coupling gradient_correspondence sign-invariant Hungarian assignment. It outputs a oneto-one mode mapping, reordered correlation matrix, and absolute-similarity heatmap for assessing SC–FC gradient correspondence. multi-zebrafish-sc_fc_matrixQuantifies zebrafish structure–function hub coupling by Structure Function degree_correlation_polyfitcorrelating weighted structural degree with positive-only Coupling scatter_plot functional degree and fitting an ordinary least-squares line. It outputs degree arrays, correlation and regression statistics, and a scatter plot with the best-fit line. multi-zebrafishExtracts and reshapes the leading PCA score from zebrafish Structure Function structural_functional_matrixstructural-connectivity predictors, sign-aligns it to funcCoupling graph_predictor_pca_correlation- tional connectivity, and measures edge-wise Pearson corpc1_fc_matrix relation across unique region pairs. It outputs a structural PC1 matrix, coupling statistics, and a composite structure–function heatmap.

30

Name

Description

Modality

multi-zebrafish-sc_fc_matrixnode_degree_correlationdegree_atlas_map

Computes weighted structural degree and absolute functional degree across filtered zebrafish brain regions, zscores them, and maps them onto opposite hemispheres of the dorsal Mapzebrain atlas. It outputs degree arrays, an anatomical hub map, and the structural–functional degree correlation. Creates an edge-wise scatter plot of structural PC1 versus functional connectivity and fits an ordinary least-squares linear regression, reporting Pearson r, R2, slope, and intercept. It visualizes structure–function coupling and the proportion of FC variance explained by structural PC1.

Structure Function Coupling

multi-zebrafish-sc_fc_matrixcorrelation_polyfit-scatter_plot

Structure Function Coupling

Supplementary Table 2 | Source publications used in BrainArena. Only publications contributing tasks to the current 32-task benchmark are included. Num. tasks denotes the number of figure-reproduction tasks drawn from each publication. Title

Venue

Year

A brain-wide map of neural activity during complex behaviour

Nature

2025 mouse

2

Neuropixels electrophysiology, Brain-wide modulation mapping; visual decision behavior and population trajectory analysis brain-atlas annotations

Brain-wide representations of prior information in mouse decision-making

Nature

2025 mouse

4

Neuropixels electrophysiology, widefield calcium imaging and block-structured visual decision behavior

Neural prior decoding; neurometric curve construction; Bayesian model comparison; neural– behavioral decay correlation

Plastic landmark anchoring in zebrafish compass neurons

Nature

2026 zebrafish

3

Calcium imaging of headdirection neurons in virtualreality visual environments

Compass alignment quantification; circular tuning-curve analysis; ring-attractor phase mapping

Brain-wide representations of behavior spanning multiple timescales and states in C. elegans

Cell

2023 C. elegans 5

Brain-wide calcium imaging and freely moving behavioral recordings

Multi-timescale filtering; timescale distribution characterization; temporal-offset decoding; manifold and PCA population analyses

Whole-brain co-mapping of bioRxiv gene expression and neuronal activity at cellular resolution in behaving zebrafish

2026 zebrafish

Whole-brain calcium imaging, Molecular subtype classification; multiplexed gene-expression marker-gene spatial distribution; mapping, behavior and visual-response profiling anatomical coordinates

Multimodal cue integration and Nature learning in a neural representa- Neurotion of head direction science

2025 Drosophila 2

Population calcium imaging in multimodal virtual-reality environments

Continuous cell-type diversification in mouse visual cortex development

2025 mouse

3

Single-cell transcriptomic and Developmental expression trajecepigenomic atlas of developing tories; cell-type diversity quanvisual cortex tification; differential expression analysis

The dynamics and geometry of Nature choice in the premotor cortex

2025 macaque

4

Premotor-cortex electrophysiology during perceptual decision-making

Single-neuron PSTH characterization; spike-raster visualization; latent-dynamics model fitting; simultaneous population recording

Science Cholinergic feedback for modality- and context-specific modulation of sensory representations

2025 mouse

2

Cholinergic feedback recordings, olfactory-task behavior and chemogenetic manipulation

Odor-response classification; chemogenetic learningmanipulation analysis

Structural and genetic determi- Science nants of zebrafish functional Advances brain networks

2025 zebrafish

4

Whole-brain calcium imaging, structural-connectivity matrices and spatial/genetic annotations

Functional connectivity estimation; structural network visualization; community detection; spatial response mapping

Nature

Organism Num. Data modality tasks

3

Analyses

Cue-dependent encoding comparison; bump-width and encodingaccuracy correlation

31

Supplementary Table 3 | BrainArena task metadata.

Task Topic

Category

Multi-timescale filtering Timescale distribution characterization Temporal offset decoding

Neural Response Characteri- Easy zation Neural Response Characteri- Easy zation Neural Decoding Hard

Neural manifold embedding Population dynamics decomposition

Population Representation Analysis Population Representation Analysis

Easy

Neural prior decoding

Neural Decoding

Hard

Neurometric curve conNeural Decoding struction Bayesian model compariComputational Modeling son Cross-modal decay correla- Statistical Analysis tion

Difficulty Algorithms / methods

Easy

Hard Hard Hard

EWMA filtering, Pearson correlation, Z-score normalization ECDF computation, Median pooling, Cross-animal statistics (SD) Ridge regression (CV), Contiguous-block CV, Circular shift shuffling UMAP, Hierarchical clustering (single-linkage), Behavioral encoding classification PCA, Hierarchical clustering (single-linkage), Weighted PC loadings Ridge regression (CV), Bayes-optimal prior (HMM), Circular shift null Logistic decoder (CV), Psychometric curve fitting, Priorconditioned neurometric Bayesian model selection (BMS/BIC), Logistic regression, Exceedance probability Action kernel fitting, Exponential decay fitting, Pearson correlation

Odor response classification Statistical Analysis Chemogenetic learning Statistical Analysis manipulation

Easy Easy

Response classification (threshold-based), Fraction computation Two-way ANOVA, SEM computation, Learning curve analysis

Compass alignment quantification Circular tuning curve analysis Ring attractor phase mapping

Hard

Population vector decoding, Circular offset calculation, Wilcoxon signed-rank test Circular histogram binning, Rayleigh test, Tuning curve construction Population vector decoding, Circular correlation, Circular harmonic fitting

Cue-dependent encoding comparison Bump-accuracy correlation

Statistical Analysis

Neural Response Characteri- Hard zation Neural Response Characteri- Hard zation Statistical Analysis

Hard

Statistical Analysis

Hard

Circular variance, Von Mises bump fitting, Linear mixed-effects model, Tukey post-hoc Pearson correlation, Von Mises bump fitting, Linear regression

Brain Network & Anatomical Easy Mapping Brain Network & Anatomical Hard Mapping Brain Network & Anatomical Hard Mapping Brain Network & Anatomical Hard Mapping

Pearson correlation (pairwise), Group averaging (nan-mean), Inter-subject reproducibility Matrix symmetrization, Coordinate transformation, Weighted network construction Louvain consensus clustering (iterated), Co-assignment matrix, Ward hierarchical clustering Voxel grid construction, Monte Carlo spatial threshold, Gaussian smoothing

Molecular subtype classification Marker gene spatial distribution Visual response profiling

Molecular & Developmental Hard Analysis Brain Network & Anatomical Easy Mapping Population Representation Hard Analysis

Combinatorial gene expression classification, Log1p transformation, Cross-animal filtering Gene expression thresholding, 3D spatial projection, Marker co-expression classification Percentile thresholding, Hierarchical clustering (response similarity), Multi-stimulus heatmap sorting

Single-neuron PSTH characterization Latent dynamics model fitting Spike raster visualization

Neural Response Characteri- Easy zation Computational Modeling Hard

PSTH computation, Trial averaging, SEM calculation

Functional connectivity estimation Structural network visualization Community detection Spatial response mapping

Neural Response Characteri- Easy zation Neural Response Characteri- Easy zation

Latent dynamical model, Poisson log-likelihood, Crossvalidated model comparison, Wilcoxon signed-rank test Spike raster construction, PSTH computation, Conditiongrouped trial averaging PSTH computation, Multi-neuron simultaneous recording, Trial averaging

Brain-wide modulation Brain Network & Anatomical Hard mapping Mapping Population trajectory analy- Population Representation Hard sis Analysis

Mann-Whitney U test, FDR correction, Permutation test (blockstratified), Swanson flatmap rendering PCA, Euclidean distance (population trajectory), Permutation test, Latency/amplitude extraction

Simultaneous population recording

Continued on next page

32

Supplementary Table 3 | BrainArena task metadata. Continued. Task Topic

Category

Difficulty Algorithms / methods

Developmental trajectory Molecular & Developmental Hard tracking Analysis Cell-type diversity quantifi- Molecular & Developmental Easy cation Analysis Differential expression Molecular & Developmental Hard analysis Analysis

Mean expression per age, LOWESS smoothing, Developmental trajectory classification Unique cluster counting, Developmental time series, Event annotation Pseudobulk aggregation, Differential expression testing, FDR correction, Log2 fold-change

Supplementary Table 4 | System-level ablation of the Co-evolving Neuroscience Analysis Repertoire on BrainArena. We report task-level scores of SeekBrain with and without the Repertoire. Diff. denotes the score difference between the full system and the ablated system. Tasks are grouped by difficulty level. Task

w/ Rep.

w/o Rep.

Diff.

76 74 92 76 71 77 92 89 81 84 91 94 74

67 32 75 48 69 62 89 56 67 78 83 71 63

9 42 17 28 2 15 3 33 14 6 8 23 11

87 89 89 72 92 84 85 59 74 66 81 63 58 65 46 75 58 46 66

76 82 83 53 79 83 81 55 51 51 76 39 46 58 42 52 57 42 58

11 7 6 19 13 1 4 4 23 15 5 24 12 7 4 23 1 4 8

Easy Multi-timescale filtering Odor response classification Chemogenetic learning manipulation Ring attractor phase mapping Cue-dependent encoding comparison Bump-accuracy correlation Functional connectivity estimation Structural network visualization Marker gene spatial distribution Single-neuron PSTH characterization Spike raster visualization Simultaneous population recording Cell-type diversity quantification Hard Timescale distribution characterization Temporal offset decoding Neural manifold embedding Population dynamics decomposition Neural prior decoding Neurometric curve construction Bayesian model comparison Cross-modal decay correlation Compass alignment quantification Circular tuning curve analysis Community detection Spatial response mapping Molecular subtype classification Visual response profiling Latent dynamics model fitting Brain-wide modulation mapping Population trajectory analysis Developmental trajectory tracking Differential expression analysis

33

Supplementary Table 5 | Representative rubric-level errors after removing the Repertoire. Task

Error introduced or worsened after removing the Repertoire

Multi-timescale filtering

The overall analysis remained executable and retained the EWMA interpretation, but neuron selection became less faithful: all three example correlations were negative and were interpreted as locomotion suppression, whereas the source examples were positively correlated. The short, intermediate, and long half-decay estimates also moved farther from the reference values, slightly weakening both the timescale and biological-conclusion scores.

Odor response classification

The ablated analysis replaced the required paired pie charts with bars and produced substantially incorrect response fractions, especially an inflated suppressed fraction. It then emphasized a significant state-dependent redistribution rather than the intended similarity in the overall responsive fraction between passive and task conditions. These errors were much more severe than the pre-ablation overinterpretation of a small task-related increase.

Chemogenetic learning manipula- The ablated figure collapsed the post-CNO/no-CNO recovery condition into a continuation of the tion inhibited-group trace instead of displaying the required third group separately, and it omitted the principal significance annotations. It also replaced the two-way group-by-day ANOVA with a Mann–Whitney test, obscuring the experimental design even though the impairment and recovery conclusions remained broadly correct. Ring attractor phase mapping

The learning panel lost the slope-2 reference and was interpreted as a two-branch, 180-degree “flip-flop” state. This is the prediction that the experiment was meant to go beyond: the key result is a stretched mapping in which 180 degrees of scene orientation spans the full 360-degree bump phase. The ablated analysis also misplaced the twofold signature into the post-training epoch and did not clearly recover the post-learning return to the diagonal mapping.

Cue-dependent encoding compari- The score changed only slightly because the pre-ablation analysis had already missed the son significant dim-versus-bright comparison. After ablation, an additional methodological error appeared: HD accuracy was no longer defined as one minus circular variance (equivalently, resultant-vector length), and the pairwise comparison brackets were omitted. The erroneous conclusion that dim and bright cues do not differ therefore persisted. Bump-accuracy correlation

The ablated analysis replaced HD accuracy (resultant-vector length of the circular offset) with the absolute circular–circular correlation, and its bump-width calculation yielded implausibly narrow values. Consequently the negative across-fly relationship became only marginal and non-significant ( 𝑝 = 0.068), whereas the pre-ablation result recovered a significant negative correlation consistent with the intended narrower-bump/higher-accuracy conclusion.

Functional connectivity estimation The matrix construction remained largely correct and the total score decreased only slightly. The principal new loss was interpretive: the ablated conclusions did not clearly explain the matrix as mesoscopic functional co-activity among larval zebrafish regions, and the figure omitted the major brain-subdivision labels. The computation of region-averaged Pearson correlations itself was preserved. Structural network visualization

The ablated pipeline plotted every positive edge rather than suppressing the bottom quartile, producing a much denser and less interpretable network, and the side-view coordinate handling was inconsistent. It also failed to state clearly that edge weights estimate relative structural wiring strength. These were substantial new degradations beyond the pre-ablation presentationlevel limitations.

Marker gene spatial distribution

The ablated maps weakened the expected rostral-telencephalic enrichment of eomesa, used unverified hard-coded gene-column assignments, and implemented the lateral view with a questionable axis choice. The conclusions continued to omit lateral dorsal-pallium pvalb7, Purkinje-cell identity, and hippocampal homology; thus the pre-existing anatomical omissions remained while the spatial mapping and method became less reliable.

Single-neuron PSTH characteriza- All four ablated examples showed broadly similar rise-and-fall profiles, further reducing the tion intended contrast between canonical ramps and non-monotonic responses. The analysis also used 20-ms bins plus Gaussian smoothing instead of 75-ms windows advanced by 10 ms and gave a slightly less bounded heterogeneity conclusion. Spike raster visualization

The ablated analysis changed the PSTH construction from the specified 75-ms sliding window to 20-ms bins with Gaussian smoothing. More importantly, its interpretation largely ignored the trial-by-trial spike-timing variability visible in the rasters, discussing reaction-time ordering without extracting the raster’s main single-trial information.

Simultaneous population record- Although the ablated figure retained 14 simultaneously recorded units, incompatible temporal ing binning/smoothing and a broader time window reduced fidelity. The largest loss was interpretive: the diverse temporal profiles were no longer clearly tied to intrinsic heterogeneity within the simultaneously recorded local PMd population, producing a much lower population-heterogeneity score.

34

Task

Error introduced or worsened after removing the Repertoire

Cell-type diversity quantification The counting procedure and overall developmental trend were preserved, but the ablated interpretation became substantially less precise. It failed to describe the early decline and the multiple bouts of subcluster expansion, omitted the P20–P28 critical-period context, and suggested activity-dependent diversification without adequately maintaining the associationversus-causation boundary. Timescale distribution characteri- The ablated analysis pooled neurons and bootstrapped an artificial variability band instead of zation computing ECDFs per animal, while also reporting 42 rather than 14 animals. Its recovered distribution was truncated at about 16.5 s, losing the long-timescale tail, and it introduced speculative higher-order/state-tracking interpretations rather than clearly stating a short-timescale majority with a long tail. Temporal offset decoding

The basic decoder analysis survived, but the ablated curve became noisier, with very wide intervals and neural performance exceeding the purported best-possible reference near zero offset. The conclusion also softened the key asymmetry to “past more strongly than future” instead of clearly stating that future velocity is not reliably decoded above baseline.

Neural manifold embedding

The ablated UMAP omitted the background of unclassified neurons and appeared substantially sparser and less continuous than the pre-ablation embedding. It also became less explicit that distance in the CePNEM-derived manifold represents behavioral-encoding similarity, although the activated/inhibited gradient and continuous-landscape conclusion were retained.

Population dynamics decomposi- The ablated analysis no longer established that the displayed PCs were obtained from neural tion activity, instead taking quantities from hierarchical model parameters; it also used populationlevel tuning scalars rather than clearly defined per-neuron CePNEM categories. The combined heatmap lost fidelity, and the already incorrect discrete-module interpretation persisted rather than recovering complex mixtures with limited cluster structure. Neural prior decoding

The ablated decoder did not clearly distinguish ordinary from pseudosession-corrected 𝑅2 , reported a mismatched value, and weakened the explanation of why the inter-trial window minimizes movement contamination. The scatter orientation and visual match also degraded slightly, although decoding remained correctly localized to vlOFC/ORBvl and retained the anticipatory-prior interpretation.

Neurometric curve construction

There was little net change. The ablated analysis simplified the Bayes-optimal prior to block probability and pooled sessions, while retaining the pre-existing failure to recover the significant dim-to-bright/zero-contrast prior shift and the emerging-choice caveat. Thus the score was nearly unchanged rather than showing a distinct new conclusion-level failure caused by ablation.

Bayesian model comparison

It added an acceptable companion panel and improved the framing of action-history models, while the already present limitations—incomplete model-selection details and minor imprecision about subjective priors—remained.

Cross-modal decay correlation

No new rubric-level error was introduced; the ablated score increased from 49 to 60. It corrected the widefield correlation’s sign, but both modalities were represented by sparse or surrogate data with severely inflated correlations, incorrect statistics, and no sufficiently specific interpretation in terms of matched neural and behavioral past-action decay.

Compass alignment quantification The ablated curve lost the strong learning-specific increase and post-learning decrease in outphase time. It then attributed the change to nonspecific time-dependent degradation rather than the two-sun manipulation and experience-dependent plasticity, directly reversing the intended biological interpretation. Simplified HD-cell/bump estimation compounded this conclusion-level error. Circular tuning curve analysis

The ablated heatmaps assigned bimodality to the post-training rather than learning epoch and described learning as degraded tuning. It completely omitted the ordered rotation of tuning curves as a function of pre-learning preferred orientation, so it missed the population-level, experience-dependent remapping that distinguishes the result from simple response duplication.

Community detection

The ablated consensus used only 200 Louvain runs rather than the source-style thousands, yielding a more uniform and overly blocky matrix. It continued to present four modules as a definitive partition and became even less explicit about multilevel hierarchical organization, although the basic consensus-matrix and dendrogram construction remained intact.

Spatial response mapping

The ablated maps used point displays and inconsistent normalization that obscured the source density patterns, including the left-habenula light-responsive concentration. The conclusions omitted the distributed hindbrain motor-positive pattern, the restricted motor-negative organization, accurate dark-responsive anatomy, and light-response lateralization; the underlying exponential-kernel regressor definition was also not explained.

35

Task

Error introduced or worsened after removing the Repertoire

Molecular subtype classification

The ablated method relied directly on a precomputed binary-gene matrix instead of constructing it from continuous counts with gene- and fish-specific thresholds, and omitted spatial filtering, cross-fish reproducibility, and cell-support criteria. It also failed to report the distribution of subtypes by number of positive genes and overemphasized broad anatomical distribution rather than non-random regional enrichment.

Visual response profiling

The ablated figure retained a large sorted population and improved some visual structure, but supplied essentially no scientific conclusion. It therefore failed to interpret the stimuluslocked ensembles, identify forward visual motion as the dominant swimming-coupled response, or connect molecular subtypes to distributed functional responses; first-trial exclusion and reproducible neuron ordering were also under-specified.

Latent dynamics model fitting

The already non-source latent model was further replaced by an incompatible latent AR(1) GLM, eliminating the remaining model-fitting/model-selection credit. Cross-validation omitted reaction-time likelihood removal and used a per-spike rather than source-scale likelihood ratio. The direction of model advantage remained correct, but medians and sample sizes were still substantially wrong.

Brain-wide modulation mapping

The ablated analysis shifted the response window from 0–100 ms to 50–150 ms, pooled neurons across sessions, and replaced session-level Fisher combination with a regional binomial/BH procedure at the wrong threshold. It consequently reported 34/223 rather than 19/201 significant regions and omitted both the concentration in named visual cortical/thalamic areas and significant non-classical regions.

Population trajectory analysis

The overall score was nearly unchanged, but the ablated trajectories had non-zero baselines and non-standard normalization. Its latency permutation did not fully preserve correlated task strata, it placed LP too late, and it omitted VISpm/VISam and the later significant midbrain/hindbrain sequence; thus the intended region-by-region temporal hierarchy remained incomplete.

Developmental trajectory tracking Both versions omitted prenatal-antecedent mapping and GAM smoothing and under-described within-class similarity, Tbr1 layer divergence, L5 ET downregulation, and complex transcriptionfactor trajectories; the ablated version modestly improved figure structure without resolving these persistent errors. Differential expression analysis

The ablated analysis retained the pre-existing wrong developmental windows and non-source Wilcoxon workflow, but its biological interpretation worsened: it characterized deep-layer glutamatergic neurons as downregulation-dominant, contrary to the expected neuronal upregulation bias. It also failed to generalize downregulation across non-neuronal subclasses, missed the strong microglial effect, and weakened the causation boundary.

Supplementary Table 6 | Brain region name. The full names of the brain region abbreviations in Figs. 3 and 4. Brain region

Abbreviation

Zebrafish brain regions Olfactory Bulb Pallium Subpallium Preoptic Region Retina Optic Tract and Accessory Optic System Epiphysis (Pineal gland) Habenula Pretectum Migrated posterior tubercular region Dorsal Thalamus Ventral Thalamus Posterior Tuberculum Rostral Hypothalamus Diffuse Nucleus of the Intermediate Hypothalamus Intermediate Hypothalamus Caudal Hypothalamus Pituitary

OB P S Po R Otr Ep H Pr mPT dTh vTh PT rH diH iH cH Pit Continued on next page

36

Supplementary Table 6 | Brain region name. (continued) Brain region Torus Longitudinalis Neuropil Layers (optic tectum) Periventricular Cell Layer (optic tectum) Torus Semicircularis Tegmentum Nucleus of Medial Longitudinal Fascicle Oculomotor Nerve Nucleus Trochlear Nerve Motor Nucleus Interpeduncular Nucleus Superior Raphe Valvula Cerebelli Corpus Cerebelli Eminentia Granularis Medial octavolateralis nucleus Rhombomere 1 Rhombomere 2 Rhombomere 3 Inferior Raphe Rhombomere 4 Medial Vestibular Nucleus Facial nerve motor nucleus and octavolateralis efferent neurons Facial nerve motor nucleus and octavolateralis efferent neurons Rhombomere 5 Rhombomere 6 Rhombomere 7 Glossopharyngeal Vagal nerve motor nucleus Rhombomere 8 Vertebra 01 ventral lateral grey matter Vertebra 01 white matter Olfactory Epithelium Trigeminal Ganglion Anterior Lateral Line Ganglion Facial Ganglion Octaval ganglion

Abbreviation TL Np Ot TS T nMLF nIII nIV IPN sR VC CC EG MON R1 R2 R3 iR R4 mVe nVII’ nVII R5 R6 R7 nX R8 lV01 wV01 OE gT gaLL gF gO

Mouse brain regions in the Allen Common Coordinate Framework template gigantocellular reticular nucleus magnocellular reticular nucleus lateral preoptic area posterior amygdalar nucleus central lateral nucleus of the thalamus main olfactory bulb facial motor nucleus substantia nigra, reticular part primary somatosensory area, nose visceral area olivary pretectal nucleus pontine reticular nucleus, caudal part intertrigeminal nucleus primary auditory area inferior colliculus nucleus of the brachium of the inferior colliculus anterolateral visual area

GRN MARN LPO PA CL MOB VII SNr SSp-n VISC OP PRNc I5 AUDp IC NB VISal

37