Reproducibility of in vivo electrophysiological measurements in mice

Abstract

Understanding brain function relies on the collective work of many labs generating reproducible results. However, reproducibility has not been systematically assessed within the context of electrophysiological recordings during cognitive behaviors. To address this, we formed a multi-lab collaboration using a shared, open-source behavioral task and experimental apparatus. Experimenters in 10 laboratories repeatedly targeted Neuropixels probes to the same location (spanning secondary visual areas, hippocampus, and thalamus) in mice making decisions; this generated a total of 121 experimental replicates, a unique dataset for evaluating reproducibility of electrophysiology experiments. Despite standardizing both behavioral and electrophysiological procedures, some experimental outcomes were highly variable. A closer analysis uncovered that variability in electrode targeting hindered reproducibility, as did the limited statistical power of some routinely used electrophysiological analyses, such as single-neuron tests of modulation by individual task parameters. Reproducibility was enhanced by histological and electrophysiological quality-control criteria. Our observations suggest that data from systems neuroscience is vulnerable to a lack of reproducibility, but that across-lab standardization, including metrics we propose, can serve to mitigate this.

Keywords: decision-making; electrophsyiology; mouse; neuropixels; neuroscience; open science; reproducibility; vision.

Figure 1.. Standardized experimental pipeline and apparatus; location of the repeated site.

(a) The pipeline for electrophysiology experiments. (b) Drawing of the experimental apparatus. (c) Location and brain regions of the repeated site. VISa: Visual Area A/AM; CA1: Hippocampal Field CA1; DG: Dentate Gyrus; LP: Lateral Posterior nucleus of the thalamus; PO: Posterior Nucleus of the Thalamus. Black lines: boundaries of sub regions within the five repeated site regions. (d) Gray lines: Actual repeated site trajectories shown within a 3D brain schematic. Red line: planned trajectory. (e) Raster plot of all measured neurons (including some that ultimately failed quality control) from one example session. (f) A comparison of neuron yield (neurons/channel of the Neuropixels probe) for this dataset (IBL), Steinmetz et al., 2019 (STE), and Siegle et al., 2021 (ALN) in three neural structures. Bars: SEM; the center of each bar corresponds to the mean neuron yield for the corresponding study.

Figure 1—figure supplement 1.. Detailed experimental pipeline for the Neuropixels experiment.

The experiment follows the main steps indicated in the colored squares at left in chronological order from top to bottom. Within each main step, actions are undertaken from left to right; diamond markers indicate points of quality control.

Figure 1—figure supplement 2.. Electrophysiology data quality examples.

(a) Example raster (left) and raw electrophysiology data snippet (right) for a recording that passes quality control. The black vertical line on the raster indicates the completion of the behavioral task and transition to spontaneous and passive stimulus periods. Red dots on the raw data snippets indicate detected spikes from multi-unit activity; green dots indicate spikes from ‘good’ units. Colored columns to the right of each plot indicate the brain regions at each depth as in Figure 1c and e. (b) Example raster and raw data snippets for four recordings that fail quality control; either because of the presence of epileptiform activity (top-left), pronounced drift (top-right), artifacts (bottom-left), or large number of noisy channels (bottom-right).

Figure 1—figure supplement 3.. Detailed comparison of yield between our dataset and published reference datasets.

Detailed comparison between the current dataset (IBL, left), Steinmetz et al., 2019 (STE, middle), and Siegle et al., 2021 (ALN, right) for measurements within cortex (Top; blue graphs), hippocampus (Middle; green graphs), and thalamus (Bottom; pink graphs). The plots show, in order: the number of insertions included in each dataset; the number of sites in the brain region per insertion; the number of ‘neurons’ produced by the spike sorting algorithm without any quality metric filters, per insertion; the number of neurons (no QC filters) per site; the fraction of neurons passing the refractory period metric; the fraction of neurons passing the noise cutoff metric; the fraction of neurons passing the minimum amplitude metric; and, the number of neurons passing all three QC metrics, per site (as shown in Figure 1f). In all cases, error bars show mean ± SE.

Figure 1—figure supplement 4.. Visual inspection of datasets by three observers blinded to data identity yielded similar metrics for all three studies.

Each plot reflects manually determined scores from a single observer (letters at top indicate their initials), who was blinded to the origin of the dataset at the time of scoring. Labels on the horizontal axis indicate the dataset: the current dataset (IBL, left), Steinmetz et al., 2019 (STE, middle), and Siegle et al., 2021 (ALN, right). Each point is the score assigned to a single snippet. Error bars represent mean ± SE. A two-way ANOVA was performed with rater identity and source dataset as categorical variables. A two-way ANOVA was performed with rater identity and source dataset as categorical variables. We found a significant effect of rater ID (p < 10^-6), and a significant interaction between rater ID and dataset (p < 0.05), reflecting statistically reliable differences between datasets that were inconsistent from rater to rater. There was no significant effect of dataset (p = 0.08) overall.

Figure 2.. Histological reconstruction reveals resolution limit of probe targeting.

(a) Histology pipeline for electrode probe track reconstruction and assessment. Three separate trajectories are defined per probe: planned, micro-manipulator (based on the experimenter’s stereotaxic coordinates), and histology (interpolated from tracks traced in the histology data). (b) Example tilted slices through the histology reconstructions showing the repeated site probe track. Plots show the green auto-fluorescence data used for CCF registration and red cm-DiI signal used to mark the probe track. White dots show the projections of channel positions onto each tilted slice. Scale bar: 1 mm. (c) Histology probe trajectories, interpolated from traced probe tracks, plotted as 2D projections in coronal and sagittal planes, tilted along the repeated site trajectory over the allen CCF, colors: laboratory. Scale bar: 1 mm. (d, e, f) Scatterplots showing variability of probe placement from planned to: micro-manipulator brain surface insertion coordinate (d, targeting variability, N=88), histology brain surface insertion coordinate (e, geometrical variability, N=98), and histology probe angle (f, angle variability, N=99). Each line and point indicates the displacement from the planned geometry for each insertion in anterior-posterior (AP) and mediolateral (ML) planes, color-coded by institution. (g, h, i) Assessment of probe displacement by institution from planned to: micro-manipulator brain surface insertion coordinate (g, N=88), histology brain surface insertion coordinate (h, N=98), and histology probe angle (i, N=99). Kernel density estimate plots (top) are followed by boxplots (bottom) for probe displacement, ordered by descending median value. A minimum of four data points per institution led to their inclusion in the plot and subsequent analysis. Dashed vertical lines display the mean displacement across institutions, indicated in the respective scatterplot in (d, e, f).

Figure 2—figure supplement 1.. Tilted slices of histology along the probe insertion for all insertions used in assessing probe placement.

Histology slices were tilted to be off-coronal along the linearly interpolated best-fit to the probe insertion vector (green: auto-fluorescence data for image registration; red: CM-DiI fluorescence signal marking probe tracks). Traced probe tracks are highlighted in white. Note that tissue damage at the brain surface visible in many sections generally occurs at the time of perfusion and brain removal, as scar tissue related to the craniotomy is detached from the brain. Scale bar: 1mm.

Figure 3.. Electrophysiological features are mostly reproducible across laboratories.

(a) Number of experimental sessions recorded; number of sessions used in analysis due to exclusion criteria. Upper branches: exclusions based on RIGOR criteria (Table 1); lower branches: exclusions based on experiment-specific criteria. For the rest of this figure, an additional targeting criterion was used in which an insertion had to hit two of the target regions to be included. (b) Power spectral density between 20 and 80 Hz of each channel of each probe insertion (vertical columns) shows reproducible alignment of electrophysiological features to histology. Insertions are aligned to the boundary between the dentate gyrus and thalamus. Tildes indicate that the probe continued below –2.0 mm. Location of channels reflect the data after the probes have been adjusted using electrophysiological features CSHL: Cold Spring Harbor Laboratory [(C): Churchland lab, (Z): Zador lab], NYU: New York University, SWC: Sainsbury Wellcome Centre, UCL: University College London, UCLA: University of California, Los Angeles, UW: University of Washington. (c) Firing rates of individual neurons according to the depth at which they were recorded. Colored blocks indicate the target brain regions of the repeated site, grey blocks indicate a brain region that was not one of the target regions. If no block is plotted, that part of the brain was not recorded by the probe because it was inserted too deep or too shallow. Each dot is a neuron, colors indicate firing rate. Probes within each institute are sorted according to their distance from the planned repeated site location. (d) p-values for five electrophysiological metrics, computed separately for all target regions, assessing the reproducibility of the distributions over these features across labs. p-values are plotted on a log-scale to visually emphasize values close to significance. (e) A Random Forest classifier could successfully decode the brain region from five electrophysiological features (neuron yield, firing rate, LFP power, AP band RMS and spike amplitude), but could not decode lab identity. The red line indicates the decoding accuracy and the grey violin plots indicate a null distribution obtained by shuffling the labels 500 times. The decoding of lab identity was performed per brain region. (* p<0.05, *** p<0.001).

Figure 3—figure supplement 1.. Recordings that did not pass QC can be visually assessed as outliers.

(a, b) Probe plots as in Figure 3b and c. The mouse name (top) is colored according to whether the recording passed QC (green) or failed (red). The probes are ordered according to the distance from the planned repeated site from left to right. (c) Probe plots showing the power spectrum of the LFP prior to using electrophysiological features to align the location of the recording channels.

Figure 3—figure supplement 2.. High LFP power in dentate gyrus was used to align probe locations in the brain.

Power spectral density between 20 and 80 Hz recorded along each probe shown in Figure 3 overlaid on a coronal atlas slice (from the Allen Common Coordinate Framework template volume). Each coronal slice has been rotated so that the probe lies along the vertical axis. Colors correspond to probe insertions belonging to a single lab (color legend at bottom). Scale bar (red) is 1 mm.

Figure 3—figure supplement 3.. Bilateral recordings assess within- vs across-animal variance.

Bilateral recordings of the repeated site in both hemispheres show within-animal variance is often smaller than across-animal variance. (a, b) Power spectral density and neural activity of all bilateral recordings. L and R indicate the left and right probe of each bilateral recording. Each L/R pair is recorded simultaneously. The color indicates the lab (lab-color assignment as in Figure 3). (c) Comparison of within-animal variance to across-animal variance for LFP power in CA1. The across-animal variance is depicted as the distribution of all pair-wise absolute differences in LFP power between any two recordings in the entire dataset (blue shaded violin plot). The black horizontal ticks indicate where the bilateral recordings (within-animal differences) fall in this distribution. (d, e) Violin plots for firing rate and spike amplitude in CA1, similar analysis as in (c). (f) The ratio of the mean of the differences within and across animals. Whether within or across animal variance is larger is dependent on the metric and brain region; red colors indicate that within < across and green colors within > across.

Figure 3—figure supplement 4.. Values used in the decoding analysis, per metric and per brain region.

All electrophysiological features (rows) per brain region (columns) that were used in the permutation test and decoding analysis of Figure 3. Each dot is a recording, colors indicate the laboratory.

Figure 4.. Neural activity is modulated during decision-making in five neural structures and is variable between laboratories.

(a) Raster plot (top) and firing rate time course (bottom) of an example neuron in LP, aligned to movement onset, split for correct left and right choices.The firing rate is calculated using a causal sliding window; each time point includes a 60ms window prior to the indicated point. (b) Peri-event time histograms (PETHs) of all LP neurons from a single mouse, aligned to movement onset (only correct choices in response to right-side stimuli are shown). These PETHs are baseline-subtracted by a pre-stimulus baseline. Shaded areas show standard error of mean (and propagated error for the overall mean). The thicker line shows the average over this entire population, colored by the lab from which the recording originates. (c,d) Average PETHs from all neurons of each lab for LP (c) and the remaining four repeated site brain regions (d). Shaded regions indicate one standard error of the mean. (e) Schematic defining all six task-modulation tests (top) and proportion of task-modulated neurons for each mouse in each brain region for an example test (movement initiation) (bottom). Each row within a region correspond to a single lab (colors same as in (d), points are individual sessions). Horizontal lines around vertical marker: SEM around mean across lab sessions. (f) Schematic to describe the power analysis. Two hypothetical distributions: first, when the test is sensitive, a small shift in the distribution is enough to make the test significant (non-significant shifts shown with broken line in grey, significant shift outlined in red). By contrast, when the test is less sensitive, the vertical line is large and a corresponding large range of possible shifts is present. The possible shifts we find usually cover only a small range. (g) Power analysis example for modulation by the stimulus in CA1. Violin plots: distributions of firing rate modulations for each lab; horizontal line: mean across lab; vertical line at right: how much the distribution can shift up- or downwards before the test becomes significant, while holding other labs constant. (h) Permutation test results for task-modulated activity and the Fano Factor. Top: tests based on proportion of modulated neurons (or in the case of Fano Factor, proportion of neurons with a value <1); Bottom: tests based on the distribution of firing rate modulations (or Fano Factor values). Comparisons performed for correct trials with non-zero contrast stimuli. (Figure analyses include collected data that pass our quality metrics and have at least four good units in the specified brain region and three recordings from the specified lab, to ensure that the data from a lab can be considered representative.).

Figure 4—figure supplement 1.. Proportion of task-modulated neurons, defined by six statistical comparisons, across mice, labs, and brain regions.

(a-f) Schematics of six statistical comparisons in specific task-related time windows to quantify neuronal task modulation. Each panel includes a schematic depiction of the time windows analyzed in the comparison, where the example trials show potential caveats of each method. For instance, in (a) the stimulus period may or may not include movement, depending on the reaction time of the animal. In (c), to avoid any overlap between pre-stimulus and reaction periods, we calculated the reaction period only in trials with >50 ms between the stimulus and movement onset. We also limited the reaction period to a maximum of 200 ms before movement onset. Below each test schematic, the corresponding proportion of task-modulated neurons across brain regions is shown, colored by lab ID (points are individual sessions; horizontal lines around vertical marker are mean ± SEM across sessions).

Figure 4—figure supplement 2.. Power analysis of permutation tests.

(a) For every test and every lab, we performed a power analysis to test how far the values of that lab would have to be shifted upwards or downwards to cause a significant permutation test (tests that were already significant in the absence of such shifts are crossed out). The p-values shown indicate the p-value of the unperturbed test. Horizontal lines indicate the means of the lab distributions and vertical bars indicate the range over which up- and downwards perturbations do not result in a significant test (p-value < 0.01). These typically show a rather small range of permissible values, which means that our permutation testing procedure is sensitive to detecting small differences in the measured quantities in individual labs. The most notable exception is the Fano Factor, meaning our ability to detect systematic differences in Fano Factor across labs is limited. (b) The correlation between the within-lab SEM of the measured quantity and the minimal shift size required to yield a significant test, showing that, as expected, measured modulations that are precisely distributed (small standard error) are very sensitive to small shifts. (c) Histogram of the magnitude of minimal shifts in units of the SD of the corresponding distribution. Most shifts are below 1 SD of the corresponding lab distribution.

Figure 5.. Principal component embedding of trial-averaged activity uncovers differences that are clear region-to-region and more modest lab-to-lab.

(a) PETHs from two example cells (black, fast reaction times only) and 2-PC-based reconstruction (red). Goodness of fit ${\displaystyle R^2}$ indicated on top with an example of a poor (top) and good (bottom) fit. (b) Histograms of reconstruction goodness of fit across all cells based on reconstruction by 1–3 PCs. Since PETHs are well approximated with just the first 2 PCs, subsequent analyses used the first 2 PCs only. (c) Two-dimensional embedding of PETHs of all cells colored by region (each dot corresponds to a single cell). (d) Mean firing rates of all cells per region, note visible pink/green divide in line with the scatter plot. Error bands are SD across cells normalized by the square root of the number of cells in the region (e) Cumulative distribution of the first embedding dimension (PC1, as shown in (c)) per region with inset of KS statistic, measuring the maximum absolute difference between the cumulative distribution of a region’s first PC values and that of the remaining cells; asterisks indicate significance at p = 0.01. (f) same data as in (c) but colored by lab. Visual inspection does not show lab clusters. (g) Mean activity for each lab, using cells from all regions (color conventions the same as in (f)). Error bands are SD across cells normalized by square root of number of cells in each lab. (h) same as (e) but grouping cells per lab. (i) p-values of all KS tests without sub-sampling; white squares indicate that there were too few cells for the corresponding region/lab pair. The statistic is the KS distance of the distribution of a target subset of cells’ first PCs to that of the remaining cells. Columns: the region to which the test was restricted and each row is the target lab of the test. Bottom row ‘all’: p-values reflecting a region’s KS distance from all other cells. Rightmost column ‘all’: p-values of testing a lab’s KS distance from all other cells. Small p-values indicate that the target subset of cells can be significantly distinguished from the remaining cells. Note that all region-targeting tests are significant while lab-targeting tests are less so. (j) Cell-number-controlled test results; for the same tests as in i the the minimum cell number across compared classes was sampled from the others before the test was computed and p-values combined using Fisher’s method across samplings. Fewer tests are significant after this control. (Figure analyses include collected data that pass our quality metrics and have at least four good units in the specified brain region and three recordings from the specified lab, to ensure that the data from a lab can be considered representative.).

Figure 5—figure supplement 1.. Lab-grouped average PETH, CDF of the first PC and 2-PC embedding, separate per brain region.

Mean firing rates across all labs per brain region (VISa/am, CA1, DG, LP, and PO shown in panels a, d, g, j, m; error bands: SD). In addition, the second column of panels panels b, (e, h, k, n) shows, for each region, the cumulative distribution function (CDF) of the first embedding dimension (PC1) per lab. The insets show the Kolmogorov-Smirnov distance (KS) per lab from the distribution of all remaining labs pooled, annotated with an asterisk if p < 0.01 for the KS test. The third column of panels (c, f, i, l, o) displays the embedded activity of neurons from VISa/am, CA1, DG, LP, and PO. For 12 region/lab pairs, dynamics were significantly different from the mean of all remaining labs using the KS test.

Figure 6.. High-firing and task-modulated LP neurons have slightly different spatial positions and spike waveform features than other LP neurons, possibly contributing only marginally to variability between sessions.

(a) Spatial positions of recorded neurons in LP. Colors: session-averaged firing rates on a log scale. To enable visualization of overlapping data points, small jitter was added to the unit locations.(b) Spatial positions of LP neurons plotted as distance from the planned target center of mass, indicated with the black x. Colors: session-averaged firing rates on a log scale. Larger circles: outlier neurons (above the firing rate threshold of 13.5 sp/s shown on the colorbar). In LP, 114 out of 1337 neurons were outliers. Only histograms of the spatial positions and spike waveform features that were significantly different between the outlier neurons (yellow) and the general population of neurons (blue) are shown (two-sample Kolmogorov-Smirnov test with Bonferroni correction for multiple comparisons; * and ** indicate corrected p-values of <0.05 and<0.01). Shaded areas: the area between 20th and 80th percentiles of the neurons’ locations. (c) (Left) Histogram of firing rate changes during movement initiation (Figure 4e, Figure 4—figure supplement 1d) for task-modulated (orange) and non-modulated (gray) neurons. 697 out of 1337 LP neurons were modulated during movement initiation. (Right) Spatial positions of task-modulated and non-modulated LP neurons, with histograms of significant features (here, z position) shown. (d) Same as (c) but using the left vs. right movement test (Figure 4e, Figure 4—figure supplement 1f) to identify task-modulated units; histogram is bimodal because both left- and right-preferring neurons are included.292 out of 1337 LP neurons were modulated differently for leftward vs. rightward movements. (Figure analyses include all collected data that pass our quality metrics, regardless of the number of recordings per lab or number of repeated site brain areas that the probes pass through.).

Figure 6—figure supplement 1.. High-firing and task-modulated VISa/am neurons.

(a–d) Similar to Figure 6 but for VISa/am. High-firing and some task-modulated VISa/am neurons had slightly different spatial positions than other VISa/am neurons, as shown, but had similar spike waveform features. Out of 877 VISa/am neurons, 550 were modulated during movement initiation (in (c)) and 170 were modulated differently for leftward vs. rightward movements (in (d)).

Figure 6—figure supplement 2.. High-firing and task-modulated CA1 neurons.

(a–d) Similar to Figure 6 but for CA1. High-firing and some task-modulated CA1 neurons had slightly different spatial positions than other CA1 neurons, as shown, but had similar spike waveform features. Out of 548 CA1 neurons, 366 were modulated during movement initiation (in (c)) and 109 were modulated differently for leftward vs. rightward movements (in (d)).

Figure 6—figure supplement 3.. High-firing and task-modulated DG and PO neurons.

Spatial positions of high-firing and some task-modulated neurons in PO, but not DG, were different from other neurons. (a, b), Spatial positions of DG neurons plotted as distance from the planned target center of mass, indicated with the black x. Spatial positions were not significantly different between high-firing neurons (yellow) and the general population of neurons (blue hues) in (a), nor between task-modulated (orange) and non-modulated (gray) neurons in (b). High-firing neurons and task-modulated neurons (for only some tests) tended to have lower spike durations, possibly related to cell subtype differences. 262 out of 448 DG neurons were modulated during movement initiation (in (b)). (c, d), Same as (a, b) but for PO neurons. Spatial positions and spike waveform features were significantly different between outliers and the general population of neurons (in (c)). Task-modulated and non-modulated PO neurons had small differences in spatial position for only some tests (not the test shown in (d)) and no difference in spike waveform features. 833 out of 1529 PO neurons were modulated during movement initiation (in (d)).

Figure 6—figure supplement 4.. Time-course and spatial position of neuronal Fano Factors.

(a) Left column: Firing rate (top) and Fano Factor (bottom) averaged over all VISa/am neurons when aligned to movement onset after presentation of left or right full-contrast stimuli (correct trials only; Fano Factor calculation limited to neurons with a session-averaged firing rate >1 sp/sec). Error bars: standard error means between neurons. Right column: Neuronal Fano Factors (averaged over 40–200 ms post movement onset after right-side full-contrast stimuli, i.e. shaded green section from the left column) and their spatial positions. Larger circles indicate neurons with Fano Factor <1. (b–e) Same as (a) for CA1, DG, LP, and PO. Spatial position between high vs. low Fano Factor neurons was only significantly different in VISa/am (deeper neurons had lower Fano Factors). In VISa/am, spike duration between high and low Fano Factor neurons was also significantly different, possibly due to cell subtype differences (neurons with shorter spike durations tended to have higher Fano Factors; histograms not shown). Lastly, in CA1 and PO, neurons with larger spike amplitudes had slightly higher Fano Factors.

Figure 6—figure supplement 5.. Neuronal subtypes and firing rates.

High-firing neurons do not belong to a specific cell subtype. To identify putative Fast-Spiking (FS) and Regular-Spiking (RS) neuronal populations,we examined spike peak-to-trough durations (Jia et al., 2019). This distribution was bimodal in Visa/am, CA1, and DG, but not LP and PO (as expected from Jia et al., 2019). This bimodality (indicated with the black and blue boxes) suggests distinct populations of FS and RS neurons only in cortical and hippocampal regions, which should have narrow (black) and wide (blue) spike widths, respectively. To confirm the distinct populations of FS and RS neurons, we next plotted the cumulative probability of firing rate for these two putative neuronal categories. Indeed, in cortex and hippocampus, neurons with narrow spikes tend to have higher firing rates (in black) while neurons with wider spikes have lower firing rates (in blue). In contrast, in LP and PO, we did not identify specific populations of neuronal subtypes using the spike waveform (to our knowledge, this has not been done in previous work either). Importantly, even in cortex/hippocampus where putative RS and FS neurons are distinguishable, there is still a large firing rate overlap between these two groups, especially for firing rates above 11–19 sp/s (the firing rate thresholds from Figure 6 and Figure 6—figure supplements 1–4). Hence, high-firing neurons do not seem to belong to only a specific neuronal subtype.

Figure 7.. Brain region, but not lab identity is decodable from single-neuron response profiles.

Results of decoding analysis for different scenarios: (a) brain region decoding using all neurons, (b) lab decoding using all neurons, and (c) lab decoding using neurons from specific brain regions. The histogram shows the distribution of macro-averaged F1 for perturbed null datasets. The null distribution is compared to the macro-averaged F1 of the original dataset (indicated by the dashed red line). The scatter plot shows the cross-validated two-dimensional linear discriminant analysis (LDA) projection. Neurons are randomly and evenly split into two groups: one for training the LDA model and the other for testing. The plots are generated exclusively using the test set.

Figure 8.. Single-covariate, leave-one-out, and leave-group-out analyses show the contribution of each (group of) covariate(s) to the MTNN model.

Lab and session IDs have low contributions to the model. (a) We adapt a MTNN approach for neuron-specific firing rate prediction. The model takes in a set of covariates and outputs time-varying firing rates for each neuron for each trial. See Table 2 for a full list of covariates. (b) MTNN model estimates of firing rates (50ms bin size) of a neuron in VISa/am from an example subject during held-out test trials. The trials with stimulus on the left are shown and are aligned to the first movement onset time (vertical dashed lines). We plot the observed and predicted PETHs and raster plots. The blue ticks in the raster plots indicate stimulus onset, and the green ticks indicate feedback times. The trials above (below) the black horizontal dashed line are incorrect (correct) trials, and the trials are ordered by reaction time. The trained model does well in predicting the (normalized) firing rates. The MTNN prediction quality measured in ${\displaystyle R^2}$ is 0.45 on held-out test trials and 0.94 on PETHs of held-out test trials. (c) We plot the MTNN firing rate predictions along with the raster plots of behavioral covariates, ordering the trials in the same manner as in (b). We see that the MTNN firing rate predictions are modulated synchronously with several behavioral covariates, such as wheel velocity and paw speed. (d) Single-covariate analysis, colored by the brain region. Each dot corresponds to a single neuron in each plot. (e) Leave-one-out and leave-group-out analyses, colored by the brain region. The analyses are run on 1133 responsive neurons across 32 sessions. The leave-one-out analysis shows that lab/session IDs have low effect sizes on average, indicating that within and between-lab random effects are small and comparable. The ‘noise’ covariate is a dynamic covariate (white noise randomly sampled from a Gaussian distribution) and is included as a negative control: the model correctly assigns zero effect size to this covariate. Covariates that are constant across trials (e.g., lab and session IDs, neuron’s 3D spatial location) are left out from the single-covariate analysis.

Figure 8—figure supplement 1.. MTNN prediction quality; MTNN slightly outperforms GLMs on predicting the firing rates of held-out test trials; PETHs and MTNN predictions for held-out test trials.

(a) For each neuron in each session, we plot the MTNN prediction quality on held-out test trials against the firing rate of the neuron averaged over the test trials. Each lab/session is colored/shaped differently. ${\displaystyle R^2}$ values on concatenations of the held-out test trials are shown on the left, and those on PETHs of the held-out test trials on the right. (b) MTNN slightly outperforms GLMs on predicting the firing rates of held-out trials when trained on movement/task-related/prior covariates. (c) The left half shows for each neuron the trial averaged activity for left choice trials and next to it right choice trials. The vertical green lines show the first movement onset. The horizontal red lines separate recording sessions while the blue lines separate labs. The right half of each of these images shows the MTNN prediction of the left half. The trial-averaged MTNN predictions for held-out test trials capture visible modulations in the PETHs.

Figure 8—figure supplement 2.. MTNN prediction quality on the data simulated from GLMs is comparable to the GLMs’ prediction quality.

The effect sizes computed by the MTNN leave-one-out analysis are similar to the effect sizes computed by the GLMs’ leave-one-out analysis. To verify that the MTNN leave-one-out analysis is sensitive enough to capture effect sizes, we simulate data from GLMs and compare the effect sizes estimated by the MTNN and GLM leave-one-out analyses. We first fit GLMs to the same set of sessions that are used for the MTNN effect size analysis and then use the inferred GLM kernels to simulate data. (a) We show the scatterplot of the GLM and MTNN predictive performance on held-out test data, where each dot represents the predictive performance for one neuron. The MTNN prediction quality is comparable to that of GLMs. (b) We run GLM and MTNN leave-one-out analyses and compare the estimated effect sizes for eight covariates. The effect sizes estimated by the MTNN and GLM leave-one-out analyses are comparable.

Figure 8—figure supplement 3.. Lab IDs have a negligible effect on the MTNN prediction.

MTNN leave-one-out analysis captures artificially inflated lab ID effects. (a) (Left) Observed per-lab PETH of held-out test trials for PO. (Center) Per-lab PETH of held-out test trials with MTNN predictions for PO. (Right) Per-lab PETH of held-out test trials with MTNN predictions for PO, after fixing the lab IDs of the input to CCU. (b) Here, we verify that the MTNN leave-one-out analysis is sensitive enough to capture variability induced by lab IDs. First, given an MTNN model trained with the full set of covariates listed in Table 2, we created four different MTNN models by perturbing the weights for the lab IDs. One is a model with the originally learned weights for the lab IDs (‘Original lab weights’). The other three models have (randomly sampled) orthogonal weights for the lab IDs, where the l2 norms of the weights for the first and last labs (i.e. the labs with lab IDs 1 and 7) are set to ${\displaystyle a}$ and ${\displaystyle b}$, respectively (‘Perturbed lab weights’, Varying l2 norm: ${\displaystyle a-b}$"). The l2 norms of the weights for the labs in between (i.e. the labs with lab IDs from 2 to 6) are set to values equally spaced between ${\displaystyle a}$ and ${\displaystyle b}$. We then generated four different sets of datasets by running the models forward with the original features that are used to train, validate and test the model. Finally, we reran the MTNN leave-one-out analysis for the lab IDs with each simulated dataset. Each dot in each column represents the change in prediction quality (measured by ${\displaystyle R^2}$) for each neuron in the dataset, when the lab IDs are left out from the set of covariates.

Figure 8—figure supplement 4.. MTNN single-covariate effect sizes are highly correlated across sessions and labs.

We plot pairwise scatterplots of MTNN single-covariate effect sizes. Each dot represents the effect sizes of one neuron and is colored by lab. Many of the effect sizes are highly correlated across sessions and labs.

Figure 9.. The decodability of task variables from the population is consistent across labs, but varies by brain region.

(a) Task-specific decoding performance per region with per session results. The decoding scores over chance level for choice, stimulus side, and reward represent decoding accuracy, while the score for wheel speed is measured in ${\displaystyle R^2}$ units. Each dot indicates the decoding performance for a single insertion, color-coded by lab identity. A one-way ANOVA is conducted to determine if there are significant differences in mean decoding scores across labs. The F-statistic and p-value are reported. After multi-comparison correction, the reported p-values reveal no statistically significant differences in decoding scores across labs. (b) The scatter plot shows the decoding score over chance level plotted against the number of neurons used. Each dot represents a single insertion, color-coded by the region (top) and lab (bottom) identity. The scores for choice, stimulus side, and reward represent decoding accuracy, while the score for wheel speed is measured in ${\displaystyle R^2}$ units. We conduct a quantitative test of the decodability of lab and region identity using decoding scores and the number of units as covariates. A support vector classifier is used, and a permutation test with 5000 permutations is applied. We report the macro F1 scores and p-values from the permutation test. After multiple-comparison correction, lab identity is less decodable than region identity, and its observed decodability may occur by chance.

Author response image 1.. PCA analysis applied to a stimulus-aligned window ([-0.5, 0.1] sec relative to stim onset).

Figure conventions as in main text Fig 5. Results are comparable to the post-movement window analysis, however regional differences are weaker here, possibly because fewer cells were active in the pre-movement window. We added panel j here and in the main figure, showing cell-number-controlled results. I.e. for each test, the minimum neuron number of the compared classes was sampled from all classes (say labs in a region), this sampling was repeated 1000 times and p-values combined via Fisher’s method, overall resulting in much fewer significant differences across laboratories and, independently, regions.

Author response image 2.. Linear mixed effects model results for two PCs and two activity windows.

For the post-movement window (“Move”), regional influences are significant (red color in plots) for all but one region while only one lab has a significant model coefficient for PC1. For PC2 more labs and three regions have significant coefficients. For the pre-movement window (“Stim”) one region for PC1 or PC2 has significant coefficients. The variance due to session id was smaller than all other effects (“eids Var”). “Intercept” shows the expected value of the response variable (PC1, PC2) before accounting for any fixed or random effects. All p-values were grouped as one hypothesis family and corrected for multiple comparisons via Benjamini-Hochberg.