Anatomically inspired digital twins capture hierarchical object representations in visual cortex

We thank the reviewer for their constructive and encouraging feedback. We are pleased that the biological motivation, structure, and empirical support of our approach were positively received. We appreciate the support for our new approach to modeling neural responses. Below, we address the specific points raised in the review.

1. Contextualizing our approach in prior literature
Previous studies have explored various strategies to understand the visual cortex, including linear-nonlinear (LN) models (Heeger, 1992; Simoncelli et al., 2004), energy models (Adelson & Bergen, 1985), and others. Deep learning has significantly advanced predictive performance, especially through convolutional neural networks (CNNs) trained on image classification tasks (e.g., Nayebi, 2023) or directly on neural data (e.g., Klindt et al., 2017; Wang et al., 2023). For a comprehensive review, see [1]. Importantly, many prior models incorporate biological constraints. For example, Antolík et al. (2016) proposed a structured model of V1 responses inspired by the feedforward visual hierarchy, where populations of V1 neurons share thalamic inputs and are organized into simple and complex-like layers. While such models have been surpassed in performance by end-to-end CNNs, the latter are often still motivated by biological principles. A key innovation, introduced by Klindt et al. [2] and adopted in our work, is the factorization of feature and spatial position into orthogonal components. Traditional approaches often treat each neuron independently, but the visual system exhibits spatial equivariance: similar computations recur across nearby locations. Leveraging this shared structure enables models to learn a common nonlinear feature space, leading to more powerful and generalizable predictions than those based on isolated neurons.
Our introduction of a hierarchical readout continues this biologically inspired modeling tradition. We first recognized that standard models, with a shared readout from the final feature layer, perform well at predicting individual neuron responses, but fail to capture population-level response structure. This motivated us to introduce a hierarchical readout variant, which explicitly models population organization by allowing structured variation in how different neurons access shared features. This approach enables the model to capture higher-order response patterns across the population, going beyond traditional single-neuron metrics.

[1] Turishcheva et al. “The Dynamic Sensorium competition for predicting large-scale mouse visual cortex activity from videos.” arXiv (2023)
[2] Klindt et al. “Neural system identification for large populations separating ‘what’ and ‘where’.” Advances in Neural Information Processing Systems (2017)

2. On prior literature about cortical hierarchies
Regarding the concern about missing citations to earlier literature on cortical hierarchies, we agree this is an important point. For mice, we will add earlier references [1,2] that demonstrate prior knowledge of visual areas beyond V1. For the anatomical hierarchy specifically, [3] was, to our knowledge, the first to systematically characterize interareal connectivity and propose a hierarchical structure based on anatomical data. Please let us know if there are additional references we may have missed. While the literature on primates is indeed more extensive and older, our study focuses on the mouse visual system; for space reasons, we chose to cite a review article [4] to provide broader context.

[1] Wang & Burkhalter. “Area map of mouse visual cortex.” J. Comp. Neurol. (2007)
[2] Marshel et al. “Functional specialization of seven mouse visual cortical areas.” Neuron (2011)
[3] Wang et al. “Network analysis of corticocortical connections reveals ventral and dorsal processing streams in mouse visual cortex.” J. Neurosci. (2012)
[4] DiCarlo et al. “How does the brain solve visual object recognition?” Neuron (2012)

3. Motivation and impact of architectural choices
For the core, we built upon the architecture introduced in the Sensorium 2023 challenge ([1]), which employed a four-layer factorized 3D convolutional network with increasing channel depth (16, 32, 64, 128). This choice was motivated by its suitability for processing dynamic stimuli with spatial and temporal resolution comparable to those used in our study. While previous models in the literature have generally used shallower cores, e.g., [2] employed a four-layer 2D CNN with 64 channels per layer; [3] and [4] used three layers with up to 48 channels, the Sensorium architecture represented the best publicly available option for mouse visual modeling and served as a practical benchmark.
Importantly, our work is among the first to systematically test how increasing core depth affects neural prediction performance. We extended the core from 4 to 8 layers and observed a marked improvement (8.4%) in single-neuron prediction accuracy. However, this modification alone did not enhance alignment with population-level hierarchical organization. This result suggests that deeper networks improve single-neuron fitting but do not resolve functional mapping across cortical areas, highlighting the necessity of architectural biases like our anatomically-inspired hierarchical readout.
For the readout, we adopted the Gaussian spatial readout proposed by [2], which currently represents the best-performing readout framework for mouse visual cortex. The dimensionality of the features used in the readout is determined directly by the channel count of the core layer being accessed, e.g., 128 features per neuron from the final layer in our case. We did not extensively test this dimensionality, but it appeared sufficient for capturing the target responses. As a point of comparison, [5] used 512 features per neuron, yet reported only marginal gains in prediction accuracy, suggesting limited sensitivity to readout dimensionality beyond a certain threshold.
We will revise the manuscript to explicitly include this rationale and emphasize that our work systematically explores design choices that have previously been adopted without extensive justification. This includes evaluating how architectural decisions, such as core depth and readout configuration, affect both single-neuron accuracy and population-level representational structure.

[1] Turishcheva et al. “The Dynamic Sensorium competition for predicting large-scale mouse visual cortex activity from videos.” arXiv (2023)
[2] Lurz et al. “Generalization in data-driven models of primary visual cortex.” ICLR (2021)
[3] Sinz et al. "Stimulus domain transfer in recurrent models for large scale cortical population prediction on video." NeurIPS (2018)
[4] Klindt et al. “Neural system identification for large populations separating “what” and “where”.” NeurIPS (2017)
[5] Wang et al. “Foundation model of neural activity predicts response to new stimulus types.” Nature (2025)

4. Clarifications on minor points
We will briefly clarify the use of the ELU+1 activation function in the main text. This ensures strictly positive activations, which is more biologically realistic given that firing rates are inherently non-negative.
Regarding the input video resolution (36 × 64 pixels), this reflects a trade-off between biological plausibility and computational efficiency. The low resolution aligns with the limited visual acuity of mice and is further supported by evidence that lower-resolution training can improve neural response predictivity (e.g., Fig. 3 in [1]). Notably, our model achieves single-neuron prediction performance comparable to that reported in [2], which also modeled responses from the MICrONS dataset but used higher resolution (144×256 pixels), enabled by extensive pretraining across multiple mice. This comparison supports the conclusion that lower resolution is sufficient for accurate modeling. We will clarify this rationale in Supplement A2.

[1] Nayebi et al. “Mouse visual cortex as a limited resource system that self-learns an ecologically-general representation.” PLOS Comp Bio (2023)
[2] Wang et al. “Foundation model of neural activity predicts response to new stimulus types.” Nature (2025)

Once again, we thank the reviewer for their valuable feedback and are confident that the proposed revisions will strengthen the manuscript.

Many thanks to the authors for their comprehensive answers to most of my comments.

I still have one theoretical question that is not completely answered by the authors' rebuttal. From the viewpoint of NN-based learning it seems rather expected that an architecture with sufficient complexity that is trained end-to-end can provide good fits for neural data, while not necessarily providing good matching of the behavior of neural units at intermediate levels of the hierarchy. That this matching becomes better by constraining the network by additional data that constrains also the behavior of units at intermediate levels of the model seems not very surprising. Apparently, in this field people so far chose to focus on end-to-end training, not considering such additional constraints. It would be very interesting to discuss the theoretical motivation why researchers before thought that such end-to-end training is a suitable approach. If this is based on theoretical assumptions about, for example, transformations of feature spaces which fall out generically from end-to-end training of deep architectures that make additional intermediate layer constraints redundant? If this should be the case, the presented work would provide evidence against the validity of these theoretical assumptions, which would be potentially an interesting general theoretical insight.

We thank the reviewer for the thoughtful and thorough review. Below, we respond point-by-point to each of the reviewer’s comments and clarify the steps taken to address the concerns raised.

Major Concerns

1. Rationale for mapping V1 to layer 6, and disentangling hierarchical readout vs. increased expressivity due to depth

Our assignment of V1 to layer 6 was based on empirical validation: performance improved when moving from layer 4 to 6, with no significant gain beyond that (Fig 3). We therefore selected the simplest effective configuration that yielded high performance. It’s also worth noting that while V1 is anatomically early in the cortical hierarchy, its activity results from multiple preceding stages of nonlinear processing. In fact, several recent studies have shown that retinal circuits alone can be well-modeled by multi-layer deep networks [1]. Additionally, V1 responses to thalamic inputs are shaped by complex recurrent dynamics. Given this, it is not surprising that deeper layers in artificial neural networks are required to best model V1 responses.

We agree that disentangling the effect of hierarchical readout from increased depth is crucial. This is addressed in Fig 3, which compares models with and without hierarchical readouts at matched depth, and in Fig 4, which controls for depth while varying readout organization. These analyses demonstrate that the performance gains stem from the anatomically inspired hierarchy rather than from depth increase.

[1] Tanaka H, et al. From deep learning to mechanistic understanding in neuroscience: the structure of retinal prediction. Advances in Neural Information Processing Systems (2019)

2. Does the model preserve single-neuron precision while capturing population-level structure?

We agree that strong single-neuron predictive accuracy remains important. To address this, we expanded training to include all 14 MICrONS sessions (rather than only 6 as in our initial submission). With the full dataset (Table 1), our 8-layer hierarchical readout model matches the performance of depth-matched models with non-hierarchical readouts, in which all areas are predicted from the final layer. Note that, while both models use the same number of parameters in the readout (17,396,557), the non-hierarchical model has a substantially larger core (10,629,424 parameters) compared to the hierarchical variant (3,417,392 parameters), resulting in more total parameters. Despite this added capacity, overall performance remains equivalent across models. This performance is comparable to the current state-of-the-art (Wang et al., 2025). When analyzing performance by area (Table 2), we find that the hierarchical model improves accuracy in AL, the top of the functional hierarchy, compared to a depth-matched baseline. Performance for V1, LM, and RL remain unchanged. These findings confirm that hierarchical readout preserves, and in some cases improves, single-neuron prediction accuracy while enabling better alignment with population-level structure.

Table 1. Average single-neuron predictive performance (correlation) per model.
Values are reported as the median normalized correlation coefficient (CC_norm) across all neurons (mean ± SEM across three random seeds), along with the range of per-session medians. Wang et al. (2025) is excluded from the average due to the lack of per-neuron data.

Table 2. Median CC_norm per area for 8-layer models with and without hierarchical readout.
Values are reported as mean ± SEM across three random seeds. Data from Wang et al. (2025) are not included due to the unavailability of area-level measurements.

3. Biological realism and the lack of cascaded processing

In our current architecture, responses in each visual area are computed from the output of a shared convolutional core, with parallel linear readouts from different layers. The core itself is hierarchical and cascaded in the sense that deeper layers (e.g., layer 8 for AL) build on earlier ones (e.g., layer 7 for LM/RL). However, we do not explicitly model inter-area communication as it occurs in biological circuits, where the output of one area directly informs the next. This is a simplification made to accommodate current data-driven modeling approaches.

Explicitly routing activity from one area’s readout to another could in principle be more biologically realistic. However, we do not expect this to substantially alter the results. The transformations between layers in our shared core already approximate such computations and can be learned directly from the data. Moreover, because we observe only a subset of neurons in each area, constraining downstream areas to use only observed upstream activity might limit performance. In contrast, accessing shared internal representations allows the model to capture a more complete and abstract embedding of the population response, improving generalization.

We will clarify this modeling choice in the revised manuscript and acknowledge the potential value of future work exploring modular architectures with explicit inter-area pathways.

4. Exclusion of inhibitory neurons

The MICrONS dataset used in this study contains functional recordings exclusively from excitatory neurons, and our model is trained and evaluated based solely on that activity. This is consistent with the experimental data we aim to replicate, specifically the studies by Hoeller et al. and Froudarakis et al., which also report responses from excitatory neurons only.

Therefore, the functional properties and hierarchical organization we model arise from excitatory populations, which aligns with the scope of both the available data and the experiments we seek to reproduce. While inhibitory interneurons are known to influence receptive fields, gain control, and temporal dynamics, our model is not intended to mechanistically capture those interactions. Instead, it is a functional model that aims to learn statistical regularities in the observed excitatory responses.

Including inhibitory neurons, on the other hand, could potentially improve training by providing additional information about local circuit dynamics. However, because inhibitory neurons are typically less selective, their per-neuron contribution to predictive performance may be limited.

Finally, we emphasize that our goal is not to maximize single-neuron prediction but to reproduce the functional hierarchy observed in population-level data. Our results show that training directly on excitatory responses is sufficient to recover this structure.

5. Brain-like invariance and equivariance in DNNs

We agree with the reviewer on the importance of emergent invariance and equivariance in artificial neural networks and appreciate the opportunity to clarify how our study relates to prior findings. Farzmahdi et al. (2024) and Cheon et al. (2022) offer complementary insights, demonstrating that task-driven models, such as AlexNet, can exhibit mirror-symmetric viewpoint tuning, while the same randomly initialized architectures already possess invariance properties. These findings suggest that the emergence of invariance and equivariance is shaped not only by network architecture, but also by the training objective and the dataset.

Although our study does not directly probe invariance or equivariance, it builds upon these insights by comparing model outputs to experimental data from Hoeller et al. and Froudarakis et al. Specifically, we assess population-level representational changes across hierarchical cortical areas using digital twin models. Our objective is to examine how invariance to object transformations, such as viewpoint, scale, and position, is distributed across cortical regions, and how this invariance emerges with anatomical hierarchy.

To address the limitations of standard models, which typically use a shared readout from the final layer and can match single-unit responses but fail to capture hierarchical population geometry, we introduce an anatomically-inspired hierarchical readout. This strategy enables the model to capture higher-order response patterns across brain areas by allowing structured variation in how populations access shared features.

We recognize that a systematic dissection of invariance and equivariance across architectures, training regimes, and biological constraints represents an important future direction. We will revise the discussion section to reflect these connections more explicitly and emphasize how our findings relate to and complement the cited studies.

Minor Points

• RL vs. LM ranking in receptive field size (Figure 1B)
  Thank you for catching this inconsistency. This was a typo in the original manuscript text. The correct experimental ranking is V1 < RL < LM < AL, as shown in the figure. We will correct this in the revised version.

Thank you to the authors for their detailed responses and clarifications. The additional results on single‐neuron–level predictions address my concerns. I also believe that the discussion of emergent invariance and equivalence in hierarchical models, and its connection to prior work, will be valuable for the field. Although I have minor concerns regarding the model design and its biological plausibility, I believe this work offers valuable insights for modeling biological vision and object representations. Therefore, I support the acceptance of this paper.

We thank the reviewer for the thoughtful and thorough review. Below, we respond point by point to each concern and will incorporate clarifications and additional analyses in the revised paper.

1. Complexity bias in hierarchical readout
Our hierarchical readout is composed exclusively of simple linear projections followed by a rectification step to ensure non-negative outputs, consistent with biological firing rate constraints. It does not include any additional nonlinear components such as multilayer perceptrons (MLPs), nor does it increase the parameter count relative to standard readouts used in prior work (e.g., [1,2]). This ensures that the function class of the decoder remains unchanged and maintains the simplicity bias that has historically motivated digital twin architectures. The novelty of our approach lies in the anatomical organization of the readout: rather than pooling all brain areas into a shared readout layer, as is standard, we assign each area its own dedicated readout layer based on its hierarchical position in the cortical structure. Importantly, the term "hierarchical" in our method refers solely to this anatomical alignment and does not imply added computational complexity.

[1] Sinz, Fabian, et al. NeurIPS (2018)
[2] Willeke, Konstantin F., et al. arXiv (2022)

2. Architecture design and applicability
Our choice of a relatively compact convolutional architecture was guided by a balance between computational tractability and the need for sufficient representational power. Importantly, this architecture has been used successfully in prior digital twin work and provides strong baseline performance [1]. To go beyond earlier baselines, we extended the architecture from 4 to 8 layers and incorporated a hierarchical readout. Combined with the full MICrONS dataset (all available sessions), this version achieves performance comparable to the model of Wang et al. (2025), which, to our knowledge, is the current state of the art in modeling mouse visual cortex responses according to standard single neuron metrics. Importantly, our model is significantly smaller (number of parameters: 20,824,610 ours, vs 58,915,491 Wang 2025), is trained on a smaller number of neurons (115k ours, vs 135k Wang 2025) and training set (minutes of video: 692 ours, vs 900 Wang 2025). Despite this, our hierarchical model achieves similar performance on key single-neuron metrics (median CCnorm in session: [0.60–0.72] ours vs [0.58–0.76] Wang 2025).
To directly assess whether hierarchical readouts are essential for success in this context, we adapted the Wang et al. model to our population-level evaluation framework and re-ran the core experiments from Hoeller et al. and Froudarakis et al. In both experiments, the Wang model performs significantly worse than our hierarchical model (see Table 1 and 2 below). Specifically, while it captures the broad difference between V1 and higher areas (i.e., lower discriminability in V1), it fails to reproduce distinctions across higher order areas (i.e., AL vs. RL/LM), predicting nearly identical performance. This is inconsistent with empirical findings that show a clear AL advantage. In contrast, our hierarchical model replicates both the V1-vs-higher-area difference and correctly ranks AL above RL and LM, aligning closely with experimental data.
These findings support our central claim: hierarchical readouts provide a useful inductive bias that helps even small models capture key aspects of the visual hierarchy. Increasing model size or training data alone does not recover these features. While larger models may boost overall performance, our results suggest they still benefit from anatomically structured readouts.

Table 1. Difference in classification accuracy between areas (Hoeller)

Table 2. Difference in object discriminability between areas (Froudarakis)

[1] Turishcheva et al., arXiv (2024)

3. Network depth and parameter count
Varying the depth of the convolutional core affects the number of internal parameters, but the number of readout parameters is kept constant across all model variants (including those with hierarchical and non-hierarchical readouts). This allows us to isolate the effect of core architecture depth without conflating it with decoder complexity.

4. Methodological and design choices
In Supp. Section A we provide a rationale for our architectural and design choices, referencing prior work that supports each of them. Our goal was to follow established practices while introducing the hierarchical readout as the key novel component. We would greatly appreciate it if the reviewer could point out any specific design choices they feel were insufficiently justified or that significantly alter the narrative, so that we can address them more directly and clarify our intentions in the revised manuscript.

5. Bypassing subcortical regions
Our current focus is on explaining the increase in invariance from V1 to higher areas like AL, as shown in prior work (e.g., Froudarakis et al., Hoeller et al.). While the retina and LGN may contribute to early visual processing, digital twins without a hierarchical readout fail to capture the functional gradient across cortical areas, even if some invariance arises earlier.

6. RL and LM readout depth
Our original hierarchical readout design was informed by both performance optimization and anatomical constraints. V1 performance was best captured using a readout from the 6th layer of our convolutional core. To model downstream areas, we extended the core and placed readouts at progressively deeper layers based on anatomical hierarchy scores from [1], as shown in Fig. 1. This hierarchy follows the order V1 < LM < RL < AL. Since the separation between LM and RL is small, we approximated their depth as equivalent and used a shared readout at layer 7. We acknowledge this as a simplification.
To assess the validity of this approximation and directly address the reviewer’s concern, we trained an additional model variant that explicitly separates LM and RL readout paths. Inspired by the dorsal/ventral stream distinction in mouse visual cortex [2], this architecture splits into two branches after layer 6. The dorsal stream targets RL and AL with readouts at layers 7 and 8, respectively. The ventral stream targets LM with a separate readout at layer 7. Due to time constraints, we evaluated this variant only in the original six sessions, not the full set of 14 used in the new analysis referenced elsewhere in the response. The model yields improved performance across all areas (Table 1), and preserves the global structure of inter-area relationships (Table 2, 3). This confirms that the hierarchical readout framework is not sensitive to the specific treatment of LM and RL depth, and that our findings hold under a more anatomically refined design.

Table 1. CCnorm per area for 8-layer models

Table 2. Classification accuracy and Kendall's Tau (Hoeller)

Table 3. Mutual information and Kendall's Tau (Froudarakis)

[1] Harris et al., Nature (2019)
[2] D’Souza et al., Nat. Commun. (2022)

7. Neuron types
Our decision to treat all neurons identically reflects the limitations of the MICrONS dataset, which contains only excitatory neuron recordings. Including inhibitory neurons could enhance model training by providing additional information about local circuit dynamics. While prior work shows that model performance improves with more neurons [1], inhibitory neurons are typically less selective, so their per-neuron contribution may be smaller. Importantly, we show that single-neuron performance does not reliably predict a model’s ability to capture the functional hierarchy. Moreover, the experimental data we aimed to reproduce also focused solely on excitatory neurons. Thus, this simplification does not compromise our main findings. Extending the approach to include both excitatory and inhibitory populations is a valuable future direction but beyond the scope of this study.

[1] Lurz et al., ICLR (2021)

The authors have answered all my questions and addressed my concerns, so I am increasing the rating to 5.

We thank the reviewer for their thoughtful and encouraging comments. We appreciate the recognition of the importance of moving beyond single-unit correlations to study population-level representations.

Novelty and Relation to Prior Work

Regarding the concern about novelty, we would like to clarify a key distinction between findings in primates and mice that underpins the contribution of our work. In primates, Yamins et al. [1] demonstrated that performance-optimized deep networks trained on object recognition exhibit a clear hierarchical correspondence with the ventral visual stream: early layers best predict V1, middle layers align with V4, and deeper layers correspond to IT cortex. This provided strong evidence that task-driven models can approximate the functional hierarchy of primate vision.

However, when the same approach was applied to the mouse visual system, Cadena et al. [2] found no evidence of such hierarchical alignment. While CNN features outperformed classical subunit energy models, randomly initialized networks performed comparably to ImageNet-trained models in predicting mouse neural responses. Critically, all visual areas studied (including V1, LM, AL, and RL, which we also examine in our work) mapped similarly across layers, showing no consistent progression from lower to higher areas.

More recently, Nayebi et al. [3] extended this analysis by systematically varying model architectures, training objectives, and input statistics. They found that certain unsupervised objectives, particularly contrastive learning, outperformed both random and supervised models in predicting mouse neural activity. While these models consistently mapped V1 to early layers and higher visual areas such as LM, AL, and RL to intermediate layers, they failed to differentiate among these higher-order areas. As a result, they cannot account for the empirical findings (e.g., Froudarakis et al., Hoeller et al., both examined in this manuscript) that reveal clear differences in object recognition performance across these regions (see Fig. 1 of our manuscript).

In this context, we believe our findings are both surprising and novel. Instead of relying on predefined tasks or training objectives, we use a data-driven approach, training a model directly on neural responses. The resulting digital twin successfully recovers meaningful distinctions between visual areas, particularly in their object recognition performance. By allowing each area to be read out from an individually matched model layer, we reveal a functional hierarchy that aligns with observed behavioral and representational differences, something that task-driven models have not achieved. This suggests that models trained directly on neural data may be better suited to capturing the computational organization of the mouse visual cortex than those relying on external tasks or objectives.

Following the reviewer’s suggestion, we will add a brief section to clarify the distinction between task-driven and data-driven approaches and to acknowledge the broader relevance of such comparisons. We hope this helps reinforce the novelty and significance of our findings.

Clarification of Anatomical Hierarchy Score

In Figure 2B, the x-axis shows the anatomical hierarchy score introduced by Harris et al. [4], which assigns a continuous value to each brain area reflecting its relative position in a cortical–thalamic hierarchy. This score is computed by analyzing the laminar termination patterns of axonal projections between areas. Each projection is classified as either feedforward or feedback based on its laminar pattern, and an optimization procedure finds the mapping that yields the most self-consistent global hierarchy.

The initial score of each area is based on the average directionality of its incoming and outgoing projections, and is refined via iterative updates that account for the positions of source and target areas. This method captures the idea that, for example, feedback projections from higher to lower areas tend to terminate in deep layers, while feedforward projections from lower to higher areas terminate more superficially.

In our manuscript, we introduce this hierarchy score briefly in Figure 1A and use it throughout the text, including as the independent variable in Figure 2B, where we relate it to receptive field size, object decoding accuracy, and discriminability. We will expand the description of this score in the revised manuscript.

Minor Comments

Typos on lines 82 and 247 will be corrected in the revised manuscript.

References

[1] Yamins, D. and DiCarlo, J. Using goal-driven deep learning models to understand sensory cortex. Nature Neuroscience (2016)
[2] Cadena, S.A., et al. How well do deep neural networks trained on object recognition characterize the mouse visual system? NeurIPS (2019)
[3] Nayebi, A., et al. Mouse visual cortex as a limited resource system that self-learns an ecologically-general representation. PLOS Computational Biology (2023).
[4]: Harris et al. “Hierarchical organization of cortical and thalamic connectivity”. Nature (2019)