Reproducibility of predictive networks for mouse visual cortex

Authors: polina  turishcheva
max f burg
fabian h sinz
alexander s ecker
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 17:58:42

Abstract: Deep predictive models of neuronal activity have recently enabled several new discoveries about the selectivity and invariance of neurons in the visual cortex.
These models learn a shared set of nonlinear basis functions, which are linearly combined via a learned weight vector to represent a neuron's function.
Such weight vectors, which can be thought as embeddings of neuronal function, have been proposed to define functional cell types via unsupervised clustering.
However, as deep models are usually highly overparameterized, the learning problem is unlikely to have a unique solution, which raises the question if such embeddings can be used in a meaningful way for downstream analysis.
In this paper, we investigate how stable neuronal embeddings are with respect to changes in model architecture and initialization. 
We find that $L_1$ regularization to be an important ingredient for structured embeddings and develop an adaptive regularization that adjusts the strength of regularization per neuron.  
This regularization improves both predictive performance and how consistently neuronal embeddings cluster across model fits  compared to uniform regularization.
To overcome overparametrization, we propose an iterative feature pruning strategy which reduces the dimensionality of performance-optimized models by half without loss of performance and improves the consistency of neuronal embeddings with respect to clustering neurons.
Our results suggest that to achieve an objective taxonomy of cell types or a compact representation of the functional landscape, we need novel architectures or learning techniques that improve identifiability. 
The code is available https://github.com/pollytur/readout_reproducibility.

Main paper: 
 introduction: 
 One central idea in neuroscience is that neurons cluster into distinct cell types defined by anatomical, genetic, electrophysiological or functional properties of single cells [23, 49, 54, 59, 69, 70] . Defining a unified taxonomy of cell types across these different descriptive properties is an active area of research. The functional properties of neurons, i. e. which computations they implement on the sensory input, constitute an important dimension along which cell types should manifest. However, each neuron's function is a complex object, mapping from a high-dimensional sensory input to the output of the neuron. Classical work in neuroscience has described neuronal function by a few parameters such as tuning to orientation, spatial frequency, or phase invariance [1, 2, 5, 8, 9, 13, 25, 40] . However, when considering neurons from higher visual areas, deeper into the brain's processing network, neuronal functions become more complex and not easily described by few manually picked parameters. In addition, recent work in the mouse visual cortex demonstrated that even in early areas, neuronal functional properties are not necessarily well-described by simple properties such as orientation [47, 71] . One possible avenue to address these challenges is to use functional representations learned by data-driven deep networks [36, 37, 39, 44] trained to predict neuronal responses from sensory inputs. The central design element of these networks is a split into a core -a common feature representation shared across neurons, and a neuron specific readouttypically a linear layer with a final nonlinearity (see Fig. 1 ). These networks are state-of-the-art in predicting neuronal responses on arbitrary visual input, and several studies have used these networks as "digital twins" to identify novel functional properties and experimentally verify them in vivo [37, 47, 57, 58, 67, 68] . Because readout weights determine a compact representation of a neuron's inputoutput mapping, they can serve as an embedding of a neuron's function. Such embeddings could be used to describe the functional landscape of neurons or obtain cell types via unsupervised clustering [59] . However, because the feature representations provided by the core are likely overcomplete, there will be many readout vectors that represent approximately the same neuronal function. Thus, it is not clear how identifiable functional cell types are based on clustering these embeddings. Moreover, since early data-driven networks [22, 28] there have been several advances in readout [48, 51, 65] and network architecture [33, 36, 56, 67] to decrease the number of per-neuron parameters, incorporate additional signals such as behavior and eye movements, or add specific inductive biases, such as rotation equivariance. Currently, there is no study systematically investigating the impact of these architectural choices on the robustness and the consistency of embedding-based clustering.Here we address this question by quantifying how consistent the models are across several fits with different seeds. We measure consistency by (1) adjusted rand index (ARI) of clustering partitions across models, (2) correlations of predicted responses across models, and (3) consistency of tuning indexes describing known nonlinear functional properties of visual neurons. Our contributions are:• We show that L 1 regularization, used in early models [28, 59] , is instrumental in obtaining a structured clustering when using embeddings of newer readout mechanisms [51] .• We introduce a new adaptive regularization scheme, which improves consistency of learned neuron embeddings while maintaining close to the state-of-the-art predictive performance.• We address the identifiability problem by proposing an iterative feature pruning strategy which reduces the dimensionality of performance-optimized models by half without loss of performance and improves the consistency of neuronal embeddings.• We show that even though our innovations improve the consistency of clustering neurons, older readout mechanisms [28] reproduce neuronal tuning properties more consistently across models.Our results suggest that to achieve an objective taxonomy of cell types or a compact representation of the functional landscape, we need novel architectures or learning techniques that improve identifiability while preserving neuronal functional properties. 
background and related work: 
 Predictive models for visual cortex. Starting with the work of Antolík et al. [22] , a number of deep predictive models for neuronal population responses to natural images have been proposed. These models capture the stimulus-response function of a population of neurons by learning a nonlinear feature space shared by all neurons (Fig. 1A ; the "core"), typically implemented by a convolutional neural network [26, 28, 51] , sometimes including recurrence [36, 55] . The core can be pretrained and shared across datasets [67] ; it forms a set of basis functions spanning the space of all neuron's input-output function. The second part of these model architectures is the "readout": it linearly combines the core's features using a set of neuron-specific weights.As the dimensionality of the core's features is quite high (height × width × feature channels), several different ways have been proposed to constrain and regularize the readout, accounting to the fact that neurons do not see the whole stimuli but rather focus on a small receptive field (RF). Initially, Klindt et al. [28] proposed a factorization across space (RF location) and features (function) (Fig. 1A ; "Factorized readout"). This approach required relatively strong L 1 regularization to constrain the RF. More recently, Lurz et al. [51] further simplified the problem by learning only an (x, y) coordinate, representing the mean of a Gaussian for the receptive field location, as well as a vector of feature The models used in our study are separated into a core learning a non-linear feature representation shared across neurons, and a neuron specific linear readout. In our case, the core is a rotation-equivariant core, i.e. each learnt feature f is analytically rotated n times by 360/n degrees, resulting in f • n output channels. We apply batch norm on the learnt feature channels f only and do not learn scale and bias parameters in the core's last layer, as this would interfere with the readout regularization. Two commonly used linear readouts are factorized and Gaussian readout. The factorized readout learns spatial weights accounting for neuron's receptive field position and feature weightings -the neural embeddings. The Gaussian readout replaced the spatial weights by sampling receptive field positions from a Gaussian and at inference time samples the Gaussian's mean. B. Readout alignment: To yield rotation invariant neural embedding clusters, we align the orientations of the embedding vectors by cyclically rotating the elements in the embedding vectors to minimize the sum of pairwise distances between the rotated embeddings [46] .C. Our pipeline: For each model architecture and training configuration we train 3 models with different parameter initialization seeds. We orientation-align their embeddings, cluster them, and eventually compute all model pair-wise ARIs. weights (Fig. 1A ; "Gaussian readout"). This removed the necessity for strong L 1 regularization. Current state-of-the-art models use Gaussian readout with weak or no regularization [60, 64, 67] .Functional properties. Previous work verified these model in vivo [37, 47, 58, 67, 68] . There is a long history of work investigating neuronal functional properties experimentally, beginning with Hubel and Wiesel (1962) [1] who showed bar stimuli to cats and found that many neurons in V1 are orientation selective. Since then, a number of studies showed that many V1 neurons are -in a linear model approximation -optimally driven by a Gabor stimulus of specific orientation, size, spatial frequency, and phase [6, 20, 32, 40] . Later, non-linear phenomena were investigated, such as surround suppression and cross-orientation inhibition. Neurons that are surround suppressed can be driven by stimulating their receptive field, while an additional stimulation of their surround reduces neural activity. Interestingly, only stimulating the surround would not elicit a response. Cross-orientation inhibition reduces neural activity of some neurons if a neuron's optimal Gabor stimulus is linearly combined with a 90 degrees rotated version of that Gabor stimulus, forming a plaid [3, 5, 7, 8, 10, 14, 17] . Recently, such experiments were performed in silico with neural predictive models, allowing for large-scale analysis of neuron's tuning properties without experimental limitations [53, 59] .Network pruning. While modern deep neural networks are typically highly overparameterized, as this simplifies optimization [38, 50, 52] , overparametrization poses a problem for reproducibility of neural embeddings: model training might end in various local optima that leads to different neuronal representations even if model outputs and predictive performance are similar. Prior work has shown that learned deep neural networks can be pruned substantially: Frankle and Carbin [34] showed that up to 80% of weights can be removed without sacrificing model performance. Further work investigated various pruning strategies. For instance Li et al. [29] pruned CNN channels based on their L 1 norm, Luo and Wu [30] used entropy as the pruning criteria, and various other criteria were studied [24, 27, 31] . Interestingly, a later study showed that model performance is relatively robust against the specific pruning strategy used [35] .3 Data and model architecture Data. We used the data from the NeurIPS 2022 Sensorium Competition for our study [60] . This dataset contains responses of primary visual cortex (V1) neurons to gray-scale natural images of seven mice, recorded using two-photon calcium imaging. In addition the mice's running speed, pupil center positions, pupil dilation and its time-derivative were recorded.Model architecture. Our model (Fig. 1 ) builds upon the baseline model from the NeurIPS 2022 Sensorium competition [60] . Similar as the baseline model, our model includes the shifter subnetwork [36] to account for eye movements and correct the resulting receptive field displacement. We also concatenate the remaining behavioral variables to the stimuli as the input to the model [56, 60] . For the model core, we use a rotation-equivariant convolutional neural network, roughly following earlier work [33, 46, 59] . Compared to this prior line of work, our model includes an additional convolutional layer (as the dataset is larger). It consists of four layers with filter sizes 13, 5, 5, and 5 pixels. In all layers, we use sixteen channels and eight rotations for each of the channels, which results in 128-dimensional neuronal embeddings. The details of model and training config are provided in Appendix A.1. For the readout, we use both the factorized readout and the Gaussian readout (Fig. 1A ). The factorized readout has been used in prior work on functional cell types [46, 59] . As it does not support the shifter network, this version of the model cannot account for eye movements. 
methods: 
 Training. Following prior work [28, 33, 36, 47, 53, 58, 60, 64, 66] , we trained the model minimizing the Poisson loss L p = N -1 N i=1 (r (i) -r (i) log r(i) ) between predicted r(i) and observed r (i) spike counts for all i = 1, ..., N neurons, as it is oftentimes assumed that neuron's firing rates follow a Poisson process [15] . For factorized readouts we add a L 1 regularization of the readout mask and embeddings to the loss, i.e.L 1 = N i=0 γ ( mask i,j |m ij | + embedding k |w k |)for spatial mask m, embedding w, and regularization coefficient γ, in line with previous work [28, 33, 39, 46, 53, 59] . As the Gaussian readout selects one spatial position, the L 1 penalty term reduces toL 1 = N i=0 γ embedding k |w k |.Thus, for both models the total loss becomes Loss = L p + L 1 . For the factorized readout γ is typically chosen to maximize performance (Fig. 6 ), while the choice of γ for Gaussian readouts is the subject of this paper.Evaluation of model performance. Following earlier work [21, 33, 36, 47, 53, 60, 64, 66] , we report model's predictive performance using Pearson correlation between each measured and predicted neural activity pair across images on the test set.Evaluation of embedding consistency. Assessing how "consistent" neuronal embeddings are across model architectures and initializations is a non-trivial question. They will naturally differ across runs in absolute terms, because core and readout are learned jointly. In the context of identifying cell types, we are interested in the relative organization of the embedding space: Will the same sets of neurons consistently cluster in embedding space? To address this question, we use clustering and evaluate how frequently pairs of neurons end up in the same group by computing the adjusted rand index (ARI, [4] ). The ARI quantifies the similarity of two cluster assignments X and Y :EQUATIONHere, a i is the number of data points in cluster i of partition X, b j is the number of data points in cluster j of partition Y , n ij is the number of data points in clusters i and j of partition X and Y , respectively, and n is the total number of data points. The ARI is invariant under permutations of the cluster labels. It is one if and only if the two partitions match exactly. It is zero when the agreement between the two partitions is as good as expected by chance and, hence, can be negative [61] .For clustering, we follow the rotation-invariant clustering procedure (Fig. 1B ) introduced by Ustyuzhaninov et al. [46] and then cluster the aligned embeddings with k means. As we do not know the true underlying number of clusters, we explore a range of k = 5 up to k = 100 clusters, which covers the typical amount of cell types suggested in various works and modalities [41, 42, 45, 59 ].To visualize the neuronal embeddings, we use t-SNE [16] using the recommendations of Kobak and Berens [43] . We use a perplexity given by 1% of the size of the dataset and set early exaggeration to 15. To avoid visual clutter, we plot a randomly sampled subset of 2,000 neurons from each of the seven animals in the dataset. We use the same neurons across all ARI experiments and t-SNE visualizations.The results do not depend on the specific cluster consistency metric used (see Appendix A.9). We assessed the credibility of t-SNE visualizations in Appendix A.10.Neuronal tuning properties. Unlike the embeddings, which were not yet biologically validated, the predicted activity from our models has been consistently verified in vivo across several studies [47, 58, 59, 68] . By focusing on predicted activity, we leverage the stability and biological relevance of these outputs and relate our findings to real-world scenarios. Specifically, we investigated the predicted neural activity to parametric, artificial stimuli that are known to yield interesting phenomena and were used extensively in biological experiments (see Section 2).To quantify neuron's tuning properties, for phase or orientation tuning, we use the parameters of the optimal Gabor and change only the phase or orientation to create a set of new Gabors. We feed this set of stimuli through the model and fit sinusoidal functions to the stimulated response curve, namely f (α) = A sin(α + ϕ) + B for the phase sweep. For the orientation sweep, rotating the Gabor-based stimuli by 180-degrees will lead the same stimulus, accounted for by the additional factor of two when fitting f (α) = A sin(2α + ϕ) + B to the stimulated responses. We define the corresponding tuning indices as the ratio of the fitted amplitude A to the mean B of the curve.For cross-orientation inhibition we added an optimal Gabor and its orthogonal copy and stimulated the tuning curve by varying the contrast of the orthogonal Gabor. Analogically, for surround suppression we varied the size and the contrast of a grating based on the optimal Gabor's parameters. To quantify both indices, we followed [11, 53] and saved the highest model prediction rmax . Next, we computed the according tuning indices as 1 -rsupp /r max , where rsupp corresponds to the highest suppression/inhibition observed. See Appendix Section A.2 for parameter details.We compute tuning index consistency across three trained models only differing in their parameter initialisation before fitting. Specifically, we compare their similarity by the normalized mean absolute error N M AE(x, y)= neurons i ((|x i -y i |)/(|x i -x|))where x and y are tuning indices for different model fits and x denotes the mean across neurons.Optimal Gabor search. All tuning properties investigated in our paper are based on optimal Gabor parameters. Previous work [33, 53, 59 ] brute-force searched the optimal stimulus. Since their models relied on the factorized readout which could integrate over multiple spatial input positions, they had to show Gabor stimuli of various properties at all spatial input positions. Thus, they needed to iterate through millions of artificially generated Gabor stimuli for each neuron. Here, we introduce a technical improvement upon their work significantly speeding up the optimal Gabor search for Gaussian readouts: as this readout type selects only one spatial position per neuron it is sufficient to limit the presented stimulus canvas size to the model's receptive field size. This significantly decreased the canvas size from 36 × 64 = 2, 304 to 25 × 25 = 625 input pixels, leading to an approximately 3.7-fold speedup. Additionally, our improved approach reduces the number of floating point operations in the model's forward pass, further benefiting the computation time. Another important detail is that for the optimal Gabor search and tuning curves experiments we completely ignored the learnt neuronal position as there was only one position to choose. For all in silico experiments, we removed models' shifter pupil position prediction network [36] , as pupil position only impacted the spacial position selection, which is now limited to a single pixel. All other behavior variables we set to their training set median values, following [56, 67] . The "mask" is a 2D matrix, selecting a "receptive field" from the latent space and the "embedding" is a linear vector representing a learned neuron function. Both are weights learned in the readout. In the top, the mask learnt is much smaller, hence, regularizing a neuronal embedding is more important to reduce the L 1 penalty, while in the bottom we can keep the neuronal embedding less sparse; its impact is smaller as the mask is bigger. I: Adjusted rand index (ARI) for clustering embeddings using k-means. We take embeddings from models trained with different seeds, cluster them and compute ARI between the clusterings. Note that even clustering the same embeddings twice with different clustering initializations will result in ARI < 1 (see Appendix A.4). 
results: 
 From factorized to Gaussian readouts. We start with our observation that the structure of the neuronal embeddings depends quite strongly on the type of readout mechanism employed by the model (Fig. 2 ). Earlier studies employing the factorized readout reported clusters in the neuronal embeddings space [46, 59] and suggested that those could represent functional cell types (Fig. 2B ). However, for more recent models with a Gaussian readout [51, 60] this structure is much less clear (Fig. 2C ), although these models exhibit substantially better predictive performance (Fig. 2A ). This vizualization difference was also evident in quantitative terms: Clustering the neuronal embeddings led to significantly more consistent clusters for the factorized readout than for the Gaussian readout (Fig. 2I ; red vs. blue), independent of the number of clusters used.We hypothesized that this difference might be caused by the L 1 regularization in earlier models, which could reveal the structure in neuronal function better and lead to more consistent embeddings. Indeed, the weight distribution for the factorized model was much sparser compared to the Gaussian readout (Fig. 2G ). Increasing the L 1 regularization indeed resulted in more structured embeddingsat least qualitatively (Fig. 2D, E ) and in terms of the distribution of embeddings (Fig. 2G ; pink and green). It also came with a drop in predictive performance comparable to that of the factorized readout (Fig. 2A ). However, somewhat surprisingly, it did not lead to more consistent cluster assignments (Fig. 2I ; pink and green). This result suggests that the structure that emerges for the Gaussian readout, arises somewhat trivially by very aggressively forcing weights to be zero due to the strong L1 penalty, but this does not happen in a consistent way across models, therefore not improving the consistency of embeddings across runs. Note that the factorized readout is always L 1 -regularized -it does not work without it [28] and is very sensitive to the choice of regularization strength (Appendix A.3). This is also why we show only a single gamma value in the table.Adaptive regularization improves consistency of embeddings while maintaining predictive performance. We next asked what could explain the difference in embedding consistency between factorized and Gaussian readout. One important difference is that the factorized readout jointly regularizes the mask and the neuronal weight. As a consequence, the factorized readout could reduce the L1 penalty on the neuron weights by increasing the size of its receptive field mask (Fig. 2H ). As a consequence, the neuronal weights may not be regularized equally.Motivated by this observation, we devised a new adaptive regularization strategy for the Gaussian readout, where the strength of the L 1 regularizer is adjustable for each neuron:EQUATIONHere, γ controls the overall strength of regularization as before. The variables β n are learned coefficients for neuron n, onto which we imposed a log-normal hyperprior p(β n ) ∼ LogNormal(1, σ), achieved by adding the according loss term L β = 1/σ 2 N n=0 log(β n ). This adaptive regularization procedure introduces an additional hyperparameter, σ, that controls how far each neuron's regularization strength is allowed to deviate from the overall average controlled by γ. We found that for relatively large σ = 0.5, the model essentially learns to down-regulate all neurons' weights (Fig. 3B ) and falls back to a very mildly regularized mode countering the effect of increasing γ (Fig. 3A ; blue). When constraining the relative weights more closely around 1 (σ = 0.1), the model indeed learns to redistribute the regularization across neurons. The average weights remain closer to 1 (Fig. 3B ) and changing the overall regularization strength γ does have an effect (Fig. 3A ).The optimal model chosen by this adaptive regularization strategy (γ = 10, σ = 0.1) indeed performs well on all dimensions we explored: First, its predictive performance matches the state-of-the-art model with the Gaussian readout (Fig. 2A ). Second, its learned neuronal embeddings are more consistent that those of both factorized and uniformly regularized Gaussian readout (Fig. 2I ; purple).Pruning models post-hoc improves consistency of neuronal clusters. We saw that the consistency of neuronal clusters could be improved by adaptive regularization of the readout. However, in absolute terms the consistency is still not very high. One potential reason for this could be that deep models are typically highly overparameterized [34] . As a result, several of the feature maps of the core could be redundant, resulting in multiple ways to represent the same neuronal function. In this case, two neurons with the same function could have highly different embeddings. We therefore now investigate whether pruning trained networks could lead to more consistent neuronal embeddings across model fits. We took the following approach to pruning models:1. Train a model as usual.2. In the trained model, mask one core output convolutional channel and compute performance.3. Remove the channel for which masking decreases model performance the least.4. Finetune the model on the reduced core. 5. Repeat steps 2-4 until only one channel is left.In agreement with prior work on pruning neural nets [34] , we observe a slight performance improvement during first pruning iterations and find that 1/3 to 1/2 of the convolutional channels can be removed without affecting performance (Fig. 4A ). We selected the smallest number of channels before the performance starts to decrease. We also observed that for the more regularized models more channels should be kept after pruning (Fig. 4A ), which makes sense as severe L 1 regularization is also a feature selection mechanism. Pruning noticeably improved the consistency of clustering neurons for models with the Gaussian readout (Fig. 4B ), but interestingly only for the regularized models. However, this improvement in consistency still comes at the cost of predictive performance.Overall, the (non-pruned) model with the adaptive readout achieves the best consistency-performance trade-off (Fig. 4C , purple).Neuronal tuning properties. To explore to how improving reproducibility of neural embeddings impact neurons' actual behavior, we examined the predicted neurons' tuning with respect to four known properties of V1 neurons: orientation tuning, phase invariance, surround suppression and cross-orientation inhibition (see Section 4 for details). We chose these properties because they are well-established and widely studied nonlinear phenomena of V1 neurons [12, 17] . We performed in-silico experiments showing classical grating-based stimuli and computed neuronal tuning curves from the predicted model responses. Note that the neurons have not seen any of these stimuli during the in-vivo experiment. Hence, our models have not been explicitly trained on these stimuli. Thus, we cannot compare these predictions against ground truth. However, we can evaluate qualitatively whether the predicted tuning curves reproduce what is known from biology and we can quantify how consistently these phenomena are expressed across predictive models that just differ in their initialization (in an ideal world tuning curves would be identical across seeds).We found that different readouts gave rise to substantially different tuning properties for the same neurons (Fig. 5 , top row; one example neuron's tuning curves). Strong uniform regularization of Gaussian readouts (γ = 50 and γ = 100) made neuron's responses highly biased towards their mean response, while the non-regularized Gaussian model exhibited the biggest responses amplitude, followed by the adaptively regularized Gaussian and the factorized model, which both exhibit moderate shrinkage effects.For each of the four response properties, we computed a tuning index that measures the degree of modulation. As tuning indices measure the effect size, tuning curves with small amplitude and high baseline lead to high phase invariance, and low orientation tuning, surround suppression, and crossorientation inhibition tuning indices (distributions in Fig. 5 bottom). Non-regularized Gaussian and factorized models exhibited similar distributions. This result suggests that even though regularization leads to sparser and less overfitted embeddings, it also removes important functional properties due to reduced expressive power. The adaptively regularized models stayed in between of the nonregularized and highly regularized models, as it puts different regularization weights across various neurons, effectively combining both modes. Thus, it keeps the overall regularization coefficient γ low, offering a good combination of both viable tuning properties and preserved consistent clustering. 
discussion: 
 We performed several important reproducibility checks for the commonly used readout mechanisms of predictive models of the visual cortex. We found that the older factorized readout led to a more structured embedding space, more consistent neuronal clusters and more reproducible in-silico tuning We compared pairwise correlations of response predictions across three model fits (↑: higher is better) and the normalized mean absolute error (NMAE; Sec. 4) across tuning indeces (↓: lower is better).curves than the more recent, performance-optimized Gaussian readout. The structure of the embedding space was primarily related to the L 1 regularization of the weights. By equipping the more recent Gaussian readout with a novel, adaptive L 1 regularization, we could recover a similarly structured embedding space and achieve more consistent neuronal clusters without sacrificing predictive performance. However, all of the models with a regularized Gaussian readout exhibited strongly biased neuronal tuning curves, calling into question how faithfully they can represent neuronal tuning beyond experimentally verified properties such as maximally exciting stimuli [37, 47, 62, 63, 68] . Rigorous, biologically meaningful evaluation remains an open question, but the presence of non-linear phenomena (tuning curves) is crucial. While current methods focus on explained variance/correlation, models with similar scores can differ in explainability [53] . To our knowledge, we conducted state-of-the-art evaluations for biological meaning by testing the tuning curves.An important goal of models is that they should be resistant to small hyperparameters changes and different initial conditions in order to lead to meaningful biological conclusions. The comparably low consistency of state-of-the-art models across seeds suggests that future research is required to improve the consistency of models, potentially by improving identifiability or improved training paradigms that find robust solutions.One limitation of our work is that we focused on one core architecture (rotation-equivariant CNN). We do not have a reason to believe that other core architectures would behave fundamentally differently.One sanity check with a regular CNN core (Appendix A.6) yielded similar results. Whether the same holds for Vision Transformer cores [64] remains to be checked.In conclusion, our work To select the amount of channels, we train same models with 3 different starting seeds and estimate the mean barρ and σ ρ of performance on each pruning stage. accounting for std using the following logic: ρi -σ ρi <= ρi+1 -σ ρi+1 where i is the amount of channels pruned.The motivation behind the rotation-equivariant core is that V1 neurons are orientation selective and neurons with all preferred orientations exist. Thus, any feature is likely to be extracted in multiple rotated versions [33, 46, 59] . To control that our results are not specific to the rotation-equivariant core, we performed a few control experiments with a non-rotation-equivariant CNN core (Fig.The general pattern of results is the same: the adaptive regularization improves consistency of clustering and so does pruning feature maps. 1 but here for the Gaussian readouts the shifter was not used. We can see that with grid mean predictor undefined the shifter actually improves the consistency of responces 1, while pruning still helps, especially for the non-regularized models. However, turning off the shifter causes some displacements in the receptive fields, probably in different directions. This effect is much more present, when grid mean predictor is on as it can compensate for the shifter displacement even more. Optimal Gabors of pruned models 10 we see that for γ = 0 the model became much more consistent across seeds, suggesting that the model was originally overparametrised. Moreover, pruned γ = 0 models now select optimal Gabors, which are close to the non-pruned majority choice across models. While this is not true for the more regularized models, they now disagree more between seeds and tend to choose Gabors with higherfrequencies, which is not a biologically plausible choice. This is probably happening due to too severe regularization. However, this explains why there is a bigger differences between the tuning indexes for more regularized models. The indexes are computed on modifications of the optimal Gabors, and if the Gabors are different the responses and indexes would be different as well.A.9 Other clustering consistency metrics ARI is the most popular metric for pairwise cluster comparisons [19] , which we think is more appropriate than set-based methods. To ensure we are not missing anything, we also considered alternative metrics such as Fowlkes-Mallows, completeness, homogeneity, and v-measure (a particular case of mutual information) and repeated the analysis. The results show the same ordering of methods A.9. Do the metrics reflect biological significance? We reasoned that if neuronal embeddings determine the response function of a neuron in the model, then two neurons that have similar embeddings under one instance of a model, should also have similar embeddings under another instance of this model. Therefore, measuring similarity of clustering is a reasonable proxy. Please note that we used t-SNE only for visualization purposes to show the presence of the density modes in the data. While t-SNE plots are indeed not sufficient to make any conclusions about the data, many papers showed that they are a good tool for exploratory data analysis [72] and we additionally checked trustworthiness metric [18] after performing t-SNE. Trustworthiness shows higher scores for more regularized models with varying numbers of neighbors confirming the presence of the density modes (Fig. 14 ).It turns out that pruning does not further increase ARI for the adaptively regularized model (Fig. 15 ). Thus, regularizing neurons in a non-uniform fashion is already effective enough that additional pruning does not improve it. We speculate that this is because some noisy neurons become sufficiently regularized without pruning, while before such regularization was only achieved with pruning and regularization combined.All of the experiments were performed on a local infrastructure cluster with 8 NVIDIA RTX A5000 GPUs with 24Gb of memory each. 2 models could be trained on a single GPU simultaneously in ≈ 3 hours. Pruning experiment, without parallelization, would take a week. No extensive cpu resources are required. Optimal Gabor search, without parallelization takes around 12hours for Gaussian readout model and 3 days for the factorized readout model. We also pre-saved optimal Gabors in batches to improve the speed, which required ≈ 100Gb of memory.Our work is an important step towards obtaining a functional cell types taxonomy for the primary visual cortex neurons. Such classification can significantly improve our understanding of the brain and help on the way to understand the mechanisms and develop the treatment for different neurodegenerative disease. Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material? Answer: [Yes] Justification: The open source available dataset was used in the experiments. The code will be public upon publication.• The answer NA means that paper does not include experiments requiring code. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted.Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes] Justification: We have 2 section in the appendix specifying hyperparameters. The data splits for train-test-validation are fixed in the open version of the dataset we used. We do not provide the specific ids for 14,000 neurons used for ARI and vizualization because it is supposed to be generalizable and the specific ids are not important.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes] Justification: The 1-sigma error bars are present in the majority of the plots.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.• The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [Yes] Justification: Broader impact section in the appendix.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: The paper poses no such risks Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected? Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?Answer: [NA] Justification: The only new asset is code and it will be available upin publication with extensive documentation.Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
crowdsourcing and research with human subjects: 
 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?Answer: [NA] Justification: the paper does not involve crowdsourcing nor research with human subjects.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 
institutional review board (irb) approvals or equivalent for research with human subjects: 
 Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?Answer: [NA]Justification: The paper does not involve crowdsourcing nor research with human subjects Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.• Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.38th Conference on Neural Information Processing Systems (NeurIPS 2024). 
 
 Appendix:  


Decision: Reviewers were impressed by the paper's contributions in the problem of DNN modeling of biological neurons: the novel methods of iterative feature pruning and adaptive regularization, and well-designed and clearly motivated experiments, as well as the focus on addressing the reproducibility of learned embeddings.  Other strengths of the paper included testing of the biological plausibility of the tuning curves which enhanced the evaluation of the method in addition to more standard ML metrics.  The results on the non-rotationally invariant CNN added during the rebuttal further contribute to the potential generalizability of the results of the paper, which were originally shown only for rotationally invariant networks (due to biological considerations).

Some reviewer comments indicate that the paper could potentially benefit by adding more explanations and clarifying the motivation and technical details, which the authors have agreed to attempt in their revision.

This work has demonstrated a high level of technical sophistication and quality, and is likely to have high impact on of neuronal data modeling.  Given the size of the audience that would find this work to be of interest, I am recommending it for a spotlight presentation.

Reviews: review: paper_summary: paper_summary: The paper "Reproducibility of predictive networks for mouse visual cortex," explores the reproducibility of neuronal embeddings in the mouse visual cortex using deep predictive models. By introducing adaptive regularization and iterative feature pruning, the authors address key issues related to model overparameterization and provide a robust framework for achieving consistent functional representations of neurons. The work lays the groundwork for future research aimed at developing reliable and interpretable models for understanding neuronal function. 

The primary goal is to investigate the stability and reproducibility of neuronal function embeddings derived from deep predictive models. These models aim to predict neuronal responses to sensory inputs and have been proposed to define functional cell types via unsupervised clustering. The paper addresses the concern that deep models are often highly overparameterized, leading to multiple solutions that can represent the same neuronal function, thereby questioning the reliability of embeddings for downstream analysis.

The paper demonstrates that L1 regularization, which was used in early models, is crucial for obtaining structured and consistent neuronal embeddings when newer readout mechanisms are used. A novel adaptive regularization scheme is introduced, which adjusts the strength of regularization for each neuron. This method improves the consistency of neuronal embeddings across different model fits while maintaining predictive performance. The paper proposes an iterative feature pruning strategy to reduce the dimensionality of performance-optimized models by half without losing predictive performance. This pruning improves the consistency of neuronal embeddings concerning clustering neurons.  The consistency of neuronal embeddings is evaluated using the Adjusted Rand Index (ARI) of clustering partitions across models, correlations of predicted responses across models, and the consistency of tuning indexes describing known nonlinear functional properties of visual neurons.
strength_weakness: strength_weakness: - Bias from Regularization and Pruning: While the adaptive regularization and pruning strategies improve consistency, they may introduce biases that affect the biological validity of the neuronal embeddings. Over-regularization, for instance, can reduce the model’s expressive power and lead to less biologically plausible representations.

- The study focuses primarily on models with a rotation-equivariant convolutional neural network (CNN) core. The authors acknowledge that other core architectures, such as regular CNNs or Vision Transformers, were not evaluated. This limitation means the findings may not generalize across different model architectures.

- Despite the improvements in clustering consistency, there is a trade-off between consistency and predictive performance. The pruning and regularization strategies, while improving consistency, sometimes result in a drop in predictive performance, which is not ideal.

- The introduction of adaptive regularization adds another layer of hyperparameters (e.g., the log-normal hyperprior parameter) that need to be carefully tuned. This increases the complexity of the model training process and may require significant computational resources.

- The focus on achieving high consistency in clustering and embeddings may overshadow other important factors, such as the interpretability and biological relevance of the model outputs. Balancing consistency with these factors is crucial for developing useful predictive models in neuroscience.
questions: questions: - Is ARI a metric you developed or did it exist before? If it existed before, it requires referencing. 

- Have you considered evaluating other core architectures, such as Vision Transformers or regular CNNs, to see if the findings generalize across different model types? What specific reasons led you to choose a rotation-equivariant CNN core for this study?

- Could you elaborate on the potential biases introduced by the adaptive regularization and pruning strategies? How do you mitigate these biases? Have you explored other regularization techniques or pruning strategies that might offer a better balance between consistency and predictive performance?

- Are there specific aspects of model performance that are most critical for maintaining biological plausibility?

- The adaptive regularization introduces additional hyperparameters. How do you approach hyperparameter tuning to ensure optimal model performance? What are the computational costs associated with this tuning process? 

- Why did you choose the Adjusted Rand Index (ARI) and other metrics for evaluating consistency? Are there alternative metrics that might provide additional insights? How do you ensure that the chosen metrics accurately reflect the biological significance of the neuronal embeddings?
limitations: limitations: NA
paper_summary: paper_summary: This work studies models trained to predict the responses of neurons in visual cortex. The model has a shared multilayer network core followed by a final layer which maps the core features into individual neuron responses. The key question that this paper asks is, How reproducible are the individual neuron properties and "cell-type" clusterings inferred by such techniques? The main conclusion is that sparsity-inducing regularization is important for finding consistent cell properties. Pruning certain channels in the core is also shown to improve consistency.
strength_weakness: strength_weakness: * I found the comic sans font and cartoon diagram style distracting in Figs 1-4.
* It isn't clear that much is gained from the t-SNE plots in Fig 2. Some evidence of different clusters, yes, but it's known that these kind of plots can be deceiving.
* How do you pick the optimal regularization parameters? By what metric/cross-validation procedure? (Line 242 gives $\gamma=10, \sigma=0.1$) Validation performance is reported but I'm unclear of the splitting procedure and whether there was a train/valid/test 3-way split or not.
* Minor points are made in "questions" section
questions: questions: * Line 35 "well described" -> "well-described" typo
* Line 62 "We address the identifiability problem, by proposing" typo, remove comma
* Fig 1 caption "sampling receptive field positions form a Gaussian" typo
* Is the clustering shown in Fig 1B truly "rotation invariant" or related to cyclic permutations of the vector indices? It is unclear from the text. The figure does not depict a rotation (i.e. action by an orthogonal matrix), since the orbit is elliptical rather than circular and it is not centered at the origin.
* Fig 3A has $\sigma=01$ which should be $\sigma=0.1$
* Fig 3B axis "amount of such weights" sounds strange
limitations: limitations: The authors mention that their work is limited to a particular type of core architecture. It's also limited to being applied to just one dataset, whereas open data in other animals (or collected by other research groups, e.g. the Allen Institute's work) are available. I am not suggesting the authors do this in their revision but that they mention it as a limitation.
paper_summary: paper_summary: The authors present a systematic investigation of the use of deep neural network fits to biological neurons as the basis for neuron cell type classification. The authors explore how various factors such as regularization and model pruning can influence both the predictive model fit and the consistency of the neural clustering.
strength_weakness: strength_weakness: Overall, the biggest weakness I think is some lack of explanation. As someone who is very familiar with work in computational neuroscience but not a neuroscientist myself, there were a few things that felt under-explained and it wasn't clear if it was because I'm not the correct audience or because the authors did not provide enough context in the text. The authors do a good job of providing clear motivation for some of the fundamentals (L24-53) but then leave out context when things get a little more technical specialized (i.e. comparisons between readout mechanisms, significance of biological properties and tuning indices, etc) as a result it feels unclear what the reader is expected to know before reading.

- The font in the figures is a little unprofessional. Comic sans? In this economy?
- Figure 2G is hard to parse. I recommend splitting up the histograms
- The intuitive explanation of the factored vs gaussian vs adaptive readouts is not super clear to me. It might be worthwhile to spend more time on it, especially to motivate your proposed adaptive mechanism.
- minor typos ("initialiation" in L183)
- L216 is nonsensical. I think some words are out of order
- Confidence intervals in Table 1 would be nice to gauge significance. Also, the table format makes it a little hard to visualize. Perhaps a graphical representation would be better.
questions: questions: - Why is the adaptive model not pruned in the experiments shown in fig 4? I feel like it makes sense to include it in the pruning experiments if possible for consistent comparison across experiments.

- I think more explanation of table 1 is warranted. First, what is the ideal case? It's not clear to me that the optimal computational model has neurons with precise tuning curves for these chosen functional biological properties. If this is obvious to the authors, then perhaps there should be explicit motivation in the text. Why these biological properties and not others? Do biological neurons in HVC demonstrate each of these properties robustly (and not others)? Is it fair to expect that neurons in silica should each individually exhibit these properties? Second, the difference in NMAE of tuning indices between the factorized and the gaussian readouts seems very large. Large enough that I would expect more discussion of it in the text. How big of a difference in score is this really, is there something that can put it in context? If the biological property fits are a big deal, then is this a nail in the coffin for gaussian readout? Is gaussian readout a hack for improving overall fit while sacrificing something more important? Or is the functional biological properties just a cool trick that would be nice to have but not essential.
limitations: limitations: Limitations are clear.
paper_summary: paper_summary: This work demonstrates that overparameterized neural network models, which have many non-unique solutions, can lead to inconsistencies in representing the mouse visual cortex. It suggests a novel approach called "adaptive regularization," where the regularization parameters in the loss term are learnable rather than fixed. The study also examines other approaches, such as normal regularization and pruning, and finds that these methods are beneficial for improving representation consistency. The proposed method significantly enhances consistency. This contribution is significant to the computational neuroscience community.
strength_weakness: strength_weakness: 1. The study is limited to only one type of model, the Rotational Equivariance CNN.

2. In the Methods section, the terms "Embedding" and "Mask" are not defined. Moreover, Lp and L1 are described as functions of model parameters, but no parameters are defined in the text. This lack of clarity poses an issue for understanding. For instance, Lp could be interpreted in two different ways: as a function of the core but not the readout, or as a function of both the core and the readout. The same ambiguity applies to L1.

3. There are two types of computational neuroscience models regarding their outputs. The first type is task-optimized models, whose outputs are task-related, such as object class, object representation, or action. These models are trained to perform downstream tasks, such as supervised object classification, reconstruction, or playing games [1,2,3]. The second type is neural response fitting models, whose outputs are predicted neural responses, and they are trained to predict these responses directly. However, the reasons why the authors chose response fitting models over task-optimized models are not mentioned.

Reference

[1] Performance-optimized hierarchical models predict neural responses in higher visual cortex

[2] Unsupervised neural network models of the ventral visual stream

[3] Using deep reinforcement learning to reveal how the brain encodes abstract state-space representations in high-dimensional environments
questions: questions: 1. Oother models that, despite lacking rotational equivariance, still demonstrate high capability in predicting neural responses measured by electrodes [1]. Why did the authors choose the Rotational Equivariance CNN? What was the motivation behind this choice?

2. Figure 2A shows the models' performances. However, what about the factorized readout with L1 regularization? Does the predictive performance decrease?

3. Why did the authors choose a log-normal prior for the learnable coefficient?

4. Why did the authors not examine adaptive regularization for factorized readout models? The implementation should be straightforward, similar to what was done for Gaussian.

5. Why does the factorized model have one γ in Table 1?


Reference

[1] https://www.brain-score.org/vision/
limitations: limitations: This study investigates the consistency of neuronal properties and the prediction performance of regularized models. However, the models used by the authors are trained to predict neural responses. There is another type of computational model that is trained to perform downstream tasks, with neural-like representations emerging from the training. Future work could consider these latter models.

Comments: title: Gentle reminder to respond to rebuttal
Gentle reminder to respond to rebuttal
Response
comment: Dear Reviewer iat5,

As the reviewer-author discussion closes soon (Nov 13 11:59 pm AoE), please let us know if the author rebuttal addressed your questions and whether your maintain or modify your original score.

Gratefully,

AC
Dear Reviewer RvHn,

As the author-reviewer discussion period is about to close, I would like to know your thoughts on the author rebuttal.  Especially since the other reviewers all currently leans towards acceptance, it would be extremely informative for me to know if you still maintain your original score following the rebuttal.  I would very much appreciate if you could reply to the authors before the close of the discussion (Nov 13 11:59 pm AoE).

Gratefully,

AC
Thank you for the authors response. I will update my review accordingly.
Thanks for your response. I hope that you will include the references and clarifications that you've provided in the final form of the paper.