Hierarchical Distribution of Reward Representation in the Cortical and Hippocampal Regions

Animals

All experiments were approved by the Animal Research Ethics Committee of Tamagawa University (animal experiment protocol, H22/27-32) and were carried out in accordance with the Fundamental Guidelines for Proper Conduct of Animal Experiment and Related Activities in Academic Research Institutions (Ministry of Education, Culture, Sports, Science, and Technology of Japan) and the Guidelines for Animal Experimentation in Neuroscience (Japan Neuroscience Society). All surgical procedures were performed under appropriate isoflurane anesthesia (see below). All effort was made to minimize suffering. The procedures for the animal experiments were established in previous studies and are described in detail below (Soma et al., 2017; Rios et al., 2019). This study is based on data from channelrhodopsin-2 (ChR2)–expressing (Thy1-ChR2) transgenic rats (W-TChR2V4; N = 33 rats; male; >3 months) abundantly expressing ChR2-Venus fusion protein under the control of the Thy1.2 promoter in cortical and other neurons (Tomita et al., 2009; Saiki et al., 2018). These animals were kept in their home cage under an inverted light schedule (lights off at 9 A.M., lights on at 9 P.M.).

Surgery

Rats were handled briefly by the experimenter (10 min, twice) before the day of surgery. For head plate implantation, rats were anesthetized with isoflurane (4.5% for induction and 2.0–2.5% for maintenance; Pfizer) using an inhalation anesthesia apparatus (Univentor 400 anesthesia unit, Univentor) and placed on a stereotaxic frame (SR-10R-HT, Narishige). In addition, lidocaine jelly (AstraZeneca) was administered around surgical incisions for local anesthesia. During anesthesia, body temperature was maintained at 37°C using an animal warmer (BWT-100, Bio Research Center). The head plate (CFR-2, Narishige) was attached to the skull with small anchor screws and two combination of dental resin cements (Super-Bond C&B, Sun Medical; Unifast II, GC Corporation). Reference and ground electrodes (Teflon-coated silver wires, A-M Systems; 125 µm in diameter) were implanted above the cerebellum. Analgesics and antibiotics were applied after the operation (meloxicam, 1 mg/kg, s.c., Boehringer Ingelheim; gentamicin ointment, 0.1% ad usum externum, Merck Sharp & Dohme).

Water deprivation was started after full recovery from surgery (6 d postoperatively). The rats had ad libitum access to water during weekends, but during the rest of the week, they obtained water only by performing the task correctly. When necessary, an agar block (containing 15 ml water) was given to the rats in their home cage to maintain them at >85% of their original body weight (Soma et al., 2021; Suematsu et al., 2025).

Behavioral task

A self-initiated left–right choice task (the LR pedal task) was used in the original system (custom made by O'HARA; Fig. 1A,B) to examine the timing of neural representation of two distinct behavioral events related to action and outcome in six brain regions (M1, M2, PPC, dCA1, vCA1, and LEC). In this task, the rats had to manipulate left and right pedals with the corresponding forelimb in a head-fixed condition. They spontaneously started each trial by pushing both pedals down with both forelimbs and holding them down for a short period (“holding period,” at least 1 s; Fig. 1B). After completing the holding period, the rats had to release either the left or the right pedal, depending on the context without any instruction cue, to obtain 0.1% saccharin water (10 µl) as a reward. The reward was dispensed from the tip of a spout by a micropump with a 0.3–0.7 s delay (0.1 s steps at random). If the rats incorrectly released the other pedal (nonrewarded trial), then they were given an error sound and were not rewarded (3 kHz, 0.3 ms). If they did not complete the holding period (immature trial), then they did not receive feedback. This task consisted of two blocks, right pedal-rewarded and left pedal-rewarded blocks. Each block lasted until the rat performed >30 correct (rewarded) trials and achieved 80% correct performance in the 10 most recent trials or until 100 rewards had been obtained (Fig. 1A). The rats typically learned this task within 2 weeks (2–3 h per day).

Once the rats completed task learning, they underwent a second surgery under isoflurane anesthesia for later recording experiments. Tiny holes (1.0–1.5 mm in diameter) were made in the skull and dura mater above the vCA1 (5.0 mm posterior, 5.0 mm lateral from the bregma), the dCA1 (4.5 mm posterior, 2.0 mm lateral from the bregma), the LEC (6.0 mm posterior, 6.8 mm lateral from the bregma), the PPC (3.8 mm posterior, 2.4 mm lateral from the bregma), the M2 (3.5 mm anterior, 2.4 mm lateral from the bregma), and the M1 (1.0 mm anterior, 2.5 mm lateral from the bregma). All holes were immediately covered with silicone sealant (DentSilicone-V, Shofu) until the recording experiments.

In vivo electrophysiological recording

Extracellular multineuronal (multiple, isolated, single-unit) recordings from individual neurons were performed while the rats were performing behavioral tasks. A supportive layer of agarose gel (2% agarose-HGT, NACALAI TESQUE) was placed on the brain, and then 32-channel silicon probes (Iso_3x_tet-A32 or Iso_4x_tet-A32; NeuroNexus Technologies) were precisely inserted into one or two of six targeted brain regions. Insertions were performed using fine micromanipulators (SM-15 or SMM-200B, Narishige) at least 1 h before the start of each recording experiment.

Wide-band signals were amplified and filtered (FA64I, Multi Channel Systems; final gain, 2,000; bandpass filter; 0.5 Hz–10 kHz) through a 32-channel head stage (MPA32I, Multi Channel Systems; gain, 10). These signals were digitized at 20 kHz and recorded with three 32-channel hard-disc recorders (LX-120, TEAC) that simultaneously digitized the pedal positions tracked by angle encoders.

Spike isolation

Wide-band signal data were processed offline to isolate spike events of individual neurons in each tetrode of the silicon probes. Briefly, wide-band signals were bandpass filtered (800 Hz–3 kHz) for spike detection, and spike candidates were detected and clustered by a semiautomatic spike-sorting software, EToS (Takekawa et al., 2010, 2012). Using open-source software (Klusters clustering software and NeuroScope viewing software; Hazan et al., 2006), spike clusters were further combined, divided, and/or discarded manually to refine single-neuron clusters based on the presence of refractory periods (>2 ms) in their own autocorrelograms (ACG) and from the absence of refractory periods in cross-correlograms with other clusters. Single-neuron clusters were included for further analysis.

Spike analysis

For each neuron (spike cluster), basal spiking properties and neuronal dynamics related to behavioral task performance were analyzed using MATLAB (MathWorks; Fig. 2). The ongoing spike rate, spike duration, and spike width for individual spike clusters were defined following the previous studies (Isomura et al., 2009). Briefly, the spike duration was the time from spike onset to the first positive peak, and the spike width referred to the elapsed time above the half amplitude of the positive spike waveform.

A temporal feature in its ACG was evaluated by defining classical “ACG bias” as the median value in a time-window from 0 to +100 ms in ACG (Saiki et al., 2014). Spike burstiness and postspike suppression were assessed by two additional indexes (Saiki et al., 2018), peak ACG bias (median at 0–50 ms) and baseline ACG bias (median at 50–250 ms).

Other temporal features of spike, interspike interval (ISI), were evaluated using the conventional coefficient of variation (Cv) and local variation (Lv) as described as follows (Shinomoto et al., 2009):Cv=SDofISIsMeanISIs, where ISI is the interspike interval and SD represents the standard deviation.Lv=3n−1∑i=1n−1(Ii−Ii+1Ii+Ii+1), where Ii and Ii+1 are the ith and i + 1 ISIs and n is the number of ISIs. Both Cv and Lv adopt a value of 0 for a sequence of perfectly regular intervals and are expected to take a value of 1 for a Poisson random series of events with ISIs that are independently exponentially distributed.

Next, a perievent time histogram (PETH) was constructed to examine neuronal dynamics associated with self-initiated actions (pedal-releasing) and outcomes (reward delivery/error sound). For action-related activity, spike trains were analyzed in relation to unilateral forelimb movements during task performance (contralateral and ipsilateral trials), aligning them to the onset (0 s) of pedal release (after ≥1 s of holding; window, onset ±500 ms). The pedal's range of motion was 0–100%, with the holding area defined as 0–30% (Fig. 1B). Task progression was marked by pedal release, detected when the pedal moved beyond the holding area. To more precisely determine release onset for neuronal analysis, pedal release was defined as the moment the pedal position exceeded 5% before reaching full release (compared with the standard 30% detection threshold; see above). For reward-related activity, spike trains were aligned to the onset of the reward delivery (0 s; window, onset ±500 ms). For nonrewarded trials, this was the onset of the auditory error cue, as this cue signaled the trial's outcome to the animal. To distinguish trial types, they were referred to as action-contralateral (AC), action-ipsilateral (AI), outcome-contralateral (OC), and outcome-ipsilateral (OI).

To assess changes in spike distribution related to action or outcome, a Kolmogorov–Smirnov (KS) test was performed, as previously described (Saiki et al., 2014; Soma et al., 2017). The cumulative distribution of all spike positions within each trial was compared with a uniformly distributed set of spike positions of the same number. The KS statistic was calculated separately for rewarded and nonrewarded trials.

The firing rate change (FRc) index, which quantifies changes in action- and outcome-related activity, was calculated using the following equation:FRcindex=(FRpost−FRpre)/(FRpost+FRpre), where FR_pre is the mean firing rate (mean FR) during the period before the action or outcome (−500 to 0 ms) and FR_post is the mean FR during the period after the action or outcome (0–500 ms). This index was calculated separately for rewarded and nonrewarded trials. A positive FRc index (>0) indicates that neural activity is positively modulated by action or outcome.

Histological observations

After the recording experiments, animals were deeply anesthetized with urethane (2–3 g/kg, i.p.; NACALAI TESQUE) and transcardially perfused with cold saline followed by 4% formaldehyde in 0.1 M phosphate buffer. Whole brains were postfixed and sliced coronally into 50 µm serial sections using a microslicer (VT1000S, Leica). Electrode tracks labeled with 1,1′-dioctadecyl-3,3,3′,3′-tetramethylindocarbocyanine perchlorate (DiI, Thermo Fisher Scientific) were observed in six brain regions under a fluorescence microscope (BX51N, Olympus).

Dataset

The final dataset for this study comprised electrophysiological data from 71 recording sessions conducted in a total of 33 Thy1-ChR2 male rats. This dataset was used to examine reward-related representations across six brain regions. This comprehensive dataset was compiled by combining newly acquired data with previously published datasets (Soma et al., 2017, 2019, 2023), which is justified by the use of identical experimental setups, behavioral tasks, recording methods, and spike-sorting criteria across all experiments. The detailed breakdown of sessions and animals contributing to each region is as follows: M1, 23 sessions/16 rats, mean ± SD, 1,064 neurons, 35 ± 25 neurons/session; M2, 28 sessions/15 rats, 1,463 neurons, 44 ± 26 neurons/session; PPC, 13 sessions/8 rats, 694 neurons, 43 ± 28 neurons/session; dCA1, 12 sessions/11 rats, 744 neurons, 37 ± 29 neurons/session; vCA1, 7 sessions/5 rats, 1,096 neurons, 157 ± 45 neurons/session; and LEC, 15 sessions/15 rats, 1,619 neurons, 70 ± 107 neurons/session. It is important to note that the total number of unique animals (33) is less than the sum of rats per region because some animals contributed data from multiple regions or participated in multiple recording sessions (single-session/day). For instance, recordings were taken from one to three brain regions per day from a single animal. To construct the feature matrix, we first filtered at the trial level: we used all trials except those in which rats moved both forelimbs during the post-action period. In other words, only action trials that were clearly defined as ipsilateral or contralateral movements were included. The detailed breakdown of sample numbers and further experimental conditions is summarized in Table 1.

Machine learning classification

To determine whether each of the six domains retains information relevant to the reward task, we employed a “neuron-by-neuron” classification approach based on the dataset described above. The analysis was constructed from the spiking dynamics of individual neurons, treating each neuron as an independent statistical sample. Specifically, for the binary classification task, two distinct data points (i.e., rows in our feature matrix) were created from each neuron: one representing its averaged activity profile under the “rewarded” condition and one for the “nonrewarded” condition. Theoretically, this would result in a total sample size exactly twice the number of recorded units. However, if a neuron lacked sufficient trial data for robust feature calculation under a given condition, that data point was excluded. Therefore, the final, actual sample size used for classification, which reflects these exclusions, is reported in Table 1.

This “neuron-by-neuron” approach was intentionally chosen for two main reasons. First, from a scientific perspective, it allows us to investigate the general, nontransient firing properties (e.g., firing rate vs spike timing patterns) that characterize a typical neuron's response to reward within each brain region's population. This stands in contrast to standard population decoding strategies, which typically rely on simultaneously recorded ensembles to capture information encoded in neuronal correlations and synergistic network interactions. By treating neurons as isolated samples, our analysis specifically targeted the intrinsic information capacity of the average neuron in each region, rather than momentary collective states.

Second, from a technical and practical perspective, this approach addresses critical constraints related to trial availability and data sparsity. Specifically, due to the high behavioral performance of the rats, the number of nonrewarded (error) trials was often limited within individual sessions. This scarcity made standard simultaneous population decoding—which requires sufficient trial numbers for both conditions within each session for robust cross-validation—impractical for many sessions. By pooling neurons across sessions, we could maximize the utility of the dataset. Furthermore, this averaging approach ensures that complex temporal features (e.g., Cv, spike timing skewness) can be robustly calculated, as detailed in the feature vector construction section below. A “trial-by-trial” analysis, in contrast, would suffer from sparse firing in short time windows, leading to numerous missing values (NAs) for these critical features. The limitations of this approach and the complementary insights offered by a “trial-by-trial” analysis are discussed further in the Discussion section.

Following this data construction, binary classification models were developed and evaluated. All machine learning analyses were performed on a host computer running Windows 11 Pro. To ensure a consistent and reproducible computational environment, we executed all analyses within a Docker container. The container environment was built from a python:3.10-slim (Debian-based) image and included key libraries such as PyCaret (Ali, 2020) for model training and SHapley Additive exPlanations (SHAP; Lundberg and Lee, 2017) for model interpretation. The complete environment is defined in the Dockerfile available in the code repository. All steps described below, from data preprocessing to model comparison, model tuning, and model interpretation, were performed using PyCaret (Ali, 2020).

Feature vector construction

The feature vector for each data point (i.e., for each neuron's rewarded or nonrewarded condition, as defined in the machine learning classification section above) was composed of a wide array of 88 features (Fig. 2, Extended Data Fig. 2-1), which can be broadly categorized into three main groups for analysis.

The first group, “fundamental spiking properties,” includes features characterizing the intrinsic firing patterns and waveform of each neuron, independent of specific task events (Fig. 2, Extended Data Fig. 2-1). These task-independent features were included based on the hypothesis that a neuron's intrinsic properties could be predictive of its role in reward-related processing. For example, it is well known that interneurons often exhibit narrower spike widths and shorter ISIs compared with pyramidal cells (Isomura et al., 2009). A model could potentially leverage such combinations of intrinsic features and task-dependent activity to learn that a specific cell type is more likely to be involved in encoding reward-related information. These properties, calculated once per neuron, comprised measures of the action potential waveform (spike duration and spike width), ISI variability (Cv and Lv) used to assess firing regularity, and features from the ACG function (namely, classical ACG bias, peak ACG bias, and baseline ACG bias) used to assess temporal firing patterns like bursting or refractoriness. This provided the model with a comprehensive baseline characterization of each neuron, which was treated as constant across rewarded and nonrewarded conditions in this analysis.

The second group, “task-related spike–timing properties,” focuses on the precise timing and distribution of spikes relative to action and outcome triggers, derived from timestamp data within the ±500 ms analysis window. This category includes moments of the spike timing distribution such as the mean spike timing, SD of spike timing, skewness (asymmetry), and kurtosis (tailedness). It also incorporates the time points corresponding to the first (Q1, 25%), second (Q2, 50%, median), and third (Q3, 75%) quartiles of the spike time distribution. In addition, the KS statistic was calculated to quantify the difference between the observed cumulative distribution of spike times and a uniform distribution, assessing nonuniformity in spike timing. These timing properties were computed separately for contralateral (AC, OC) and ipsilateral (AI, OI) trials relative to the action or outcome trigger (Fig. 2, Extended Data Fig. 2-1).

The third group consists of “task-related firing-rate–based properties,” which were derived from PETHs aligned to action onset and outcome onset. These features quantify the magnitude and change in the firing rate around these task events. They included the mean FR calculated within various time windows relative to the trigger (e.g., 0–50 ms, 50–100 ms, 100–250 ms), the peak firing rate (Peak FR) observed within specific time windows, and the FRc index quantifying the relative change in firing rate between the 500 ms periods before and after the trigger (Eq. 3). Similar to the spike timing properties, these rate-based features were calculated separately for contralateral (AC, OC) and ipsilateral (AI, OI) trials relative to the action or outcome trigger, respectively (see Fig. 2 and Extended Data Fig. 2-1 for all specific windows and trial types).

The final, comprehensive feature vector for each data point was a concatenation of all the properties described above. This included task-dependent features derived from all combinations of epoch (action vs outcome) and movement direction (ipsilateral vs contralateral), i.e., AC, AI, OC, and OI conditions. This approach, rather than building separate decoders for each specific condition, allowed a single classifier to leverage all available information simultaneously—both the neuron's intrinsic properties and its diverse task-modulated responses—to find the most predictive patterns for the overall trial outcome.

To focus exclusively on neural information related to reward processing, we did not analyze activity during periods likely to be strongly influenced by forelimb movements (actions)—specifically, the post-action and pre-outcome epoch. This period includes neuronal activity immediately following pedal pressing and immediately preceding outcome delivery, which may primarily reflect the physical act of movement rather than reward-related or predictive processing (Fig. 1B). Therefore, these time windows were excluded from the final feature set used for modeling except for the calculation of the FRc (see above).

Data preprocessing

Following the data preparation phase, a series of analytical procedures were conducted using PyCaret (Ali, 2020). These procedures included preprocessing, model comparison, model training (focusing on the top three models), hyperparameter tuning (for the top three models), and model evaluation.

The dataset underwent several preprocessing steps as follows: missing value imputation was conducted by imputing missing values with the mean. Feature scaling was performed using standard scaling, specifically z-score normalization. Class imbalance was addressed through the application of synthetic minority oversampling technique (SMOTE), which generates new synthetic samples by selecting examples that are proximate in the feature space and interpolating between them. Multicollinearity was assessed by examining the absolute Pearson's correlation coefficient between any two features; if this coefficient exceeded 0.80, the features were considered highly collinear. In such cases, the feature with the lower correlation to the response variable was removed. Feature selection was implemented to reduce dimensionality, thereby minimizing overfitting and enhancing computational efficiency. An automatic feature selection process was employed, which involved the following: automatic evaluation of all features based on their importance to the target variable. A Light Gradient Boosting Machine (LightGBM) model was trained on the preprocessed data to compute feature importance scores, determined by the frequency and effectiveness of features in decision splits. LightGBM is recognized for its rapid training and high accuracy, particularly with structured data, making it a robust choice for evaluating feature importance. It is capable of capturing complex, nonlinear relationships, ensuring the retention of the most informative features. The importance scores generated by LightGBM were subsequently used to rank the features. Only the top 20% of features, as ranked by their estimated importance, were retained for model training. These selected features were then utilized in subsequent steps, including model comparison, training, tuning, and testing.

Model training and comparison

After the preprocessing and the preparation of the training data (80% of the full dataset), the performance of numerous classification algorithms was compared. These encompassed gradient boosting machines [CatBoost, Extreme Gradient Boosting (XGBoost), LightGBM], other ensemble methods [adaptive boosting (AdaBoost), decision tree, extra trees, random forest, gradient boosting], linear and discriminant models [linear discriminant analysis (LDA), logistic regression, quadratic discriminant analysis, ridge, support vector machine (SVM)], K-nearest neighbors, naive Bayes, and a dummy baseline. Each model was evaluated using a stratified 10-fold cross-validation strategy, ensuring performance was robust and generalizable. For each candidate model, key performance metrics [e.g., accuracy, F1 score, the area under the receiver operating characteristic curve (AUC), etc.] were computed during cross-validation. The results were averaged across the validation folds for each model and the models were ranked based on their performance (i.e., accuracy). A subset of models was identified for further tuning and analysis for each brain area. The candidate models were chosen based on the mean accuracy (the proportion of correctly predicted instances).

Tune model

Subsequent tuning was conducted on the subset of models previously identified. This phase involved refining hyperparameters to further optimize the models. The tuned models were subsequently reevaluated. The procedure was executed as follows: a range of hyperparameters was assessed using stratified fivefold cross-validation to determine the optimal configurations that enhance model performance based on accuracy. This process was iterated for a predetermined number of iterations (i.e., 20), sampling various hyperparameter combinations. Optuna, an open-source library for hyperparameter optimization, was employed to explore the hyperparameter space. During this iterative search, the function compared performance metrics to select the optimal set of hyperparameter values. Upon completion of the tuning process, a new model object with the best-found hyperparameters was returned. This tuned model is anticipated to deliver superior performance compared with the original model with default parameters. In the event that the tuning process does not yield improvements, the original model was retained, ensuring the selection of the best-performing model for the task. To ensure robustness, this tuning procedure was repeated three times for each model and region combination.

The tuned model evaluation

To thoroughly assess the generalization performance of the tuned models, we employed the held-out test set, comprising 20% of the total dataset, which was completely isolated from the training and tuning phases. For each brain region, the model with the optimal hyperparameters identified during tuning was retrained on the entire training set and subsequently evaluated on the test set. Importantly, to ensure the robustness of our findings and account for potential variability due to data splitting and model initialization, this entire evaluation process—from the train–test split to hyperparameter tuning and final testing—was repeated three times with different random seeds. The final performance metrics, specifically classification accuracy and the AUC, were calculated as the mean and SD across these three independent repetitions. These aggregated results reflect the model's capability to classify unseen neural data.

Statistical verification of model performance: shuffle control analysis and ANOVA

To rigorously validate that the classification performance of our models was significantly above chance and to statistically compare performance across brain regions and model architectures, we performed two distinct statistical analyses: a shuffle control analysis and a standard two-way ANOVA.

First, to validate that the classification performance of our models was significantly above chance and genuinely reflected reward-related information in the neural activity, we performed a shuffle control analysis for each brain region and each of the tuned models. For each analysis, the hold-out test dataset (20% of the data) was used. While keeping the feature vectors and the model's predictions constant, the corresponding trial outcome labels were randomly shuffled 1,000 times. For each iteration, the decoding accuracy was calculated by comparing the model's original predictions against the shuffled labels. This procedure generated a null distribution of accuracy scores, representing the performance expected by chance. To account for multiple comparisons across the 18 conditions (6 brain regions × 3 model architectures), we calculated p values based on the aggregated accuracy across the three independent repeats. A Bonferroni’s correction (correction factor, 18) was applied to the uncorrected p values.

Second, to examine the hierarchical differences in representation of reward-task–related information, we performed a standard two-way ANOVA with the brain region (six levels) and model architecture (three levels) as factors. For this analysis, individual accuracy scores from the three independent repeats were treated as observations (N = 54). This approach was adopted to explicitly incorporate the variability arising from data splitting and model initialization into the statistical evaluation, thereby ensuring that the assessed regional differences are robust to these stochastic factors. The main effect of the brain region was evaluated to determine whether classification performance consistently differed across regions regardless of the model used.

Interpretation of the final models via SHAP analysis

To interpret the predictions of the final “black box” models, SHAP was employed, a game-theoretic approach that explains the output of any machine learning model (Lundberg and Lee, 2017). Grounded in Shapley values from cooperative game theory, SHAP provides a unified framework for model interpretation with several key advantages. It not only offers global feature importance by averaging the absolute SHAP values across all samples but also provides local, instance-level explanations, revealing how each feature contributes to a specific, individual prediction. This dual capability allows for a comprehensive understanding of model behavior.

For the analysis, SHAP values were computed for each feature and data point in the test set. We visualized the global feature importance using summary plots, which display the distribution of SHAP values for each feature and rank them by their overall impact on the model. To understand the drivers of specific classifications, we used waterfall plots to illustrate how individual feature contributions sum up to push the model's output from the base value to the final prediction for representative single-neuron examples. This approach allowed us to identify the key neural properties that our models learned to associate with rewarded versus nonrewarded trial outcomes.

Since SHAP values are specific to a trained model instance, we selected a representative model for visualization purposes. First, the best-performing model architecture for each brain region was identified based on the mean test accuracy across three independent repetitions. Then, from the three repetitions of this architecture, the specific model instance achieving the highest test accuracy was selected to generate the SHAP summary plots (beeswarm and bar plots) and representative waterfall plots for example neurons. We confirmed that the top-contributing features were largely consistent across the three repetitions.

Statistical analysis

Standard replication of measurements was performed for this study. The reported findings were reproduced across animals. All quantifications were conducted at the single-neuron level. Sample sizes (the numbers of animals, sessions, and neurons) were estimated according to previous studies (Rios et al., 2019; Sakairi et al., 2025). The KS test was used to extract the test statistic, rather than the p value, as a feature for the model. These statistical tests were conducted with MATLAB's Statistics and Machine Learning Toolbox (MathWorks). Blinding and randomization were not performed. Models were tuned via cross-validation on training data, and all reported metrics and interpretability are computed on a held-out test set using the tuned (nonfinalized) model and training-fitted preprocessing.

Data availability

Data are available upon reasonable request. Data will be made available upon publication in a public repository.

Code availability

The code/software described in the paper is freely available online at https://github.com/mokamotosan/soma_okamoto_lr_01_public.git. The code is available as Extended Data 1.