Jacobian-Based Interpretation of Nonlinear Neural Encoding Model
Authors: Xiaohui Gao
Haoran Yang
Yue Cheng
Mengfei Zuo
Yiheng Liu
Peiyang Li
Xintao Hu
Venue: NeurIPS 2025 spotlight ()
Keywords: fMRI
neural encoding model
nonlinearity
Jacobian matrix
Primary area: neuroscience_and_cognitive_science

Abstract: In recent years, the alignment between artificial neural network (ANN) embeddings and blood oxygenation level dependent (BOLD) responses in functional magnetic resonance imaging (fMRI) via neural encoding models has significantly advanced research on neural representation mechanisms and interpretability in the brain. However, these approaches remain limited in characterizing the brain’s inherently nonlinear response properties. To address this, we propose the Jacobian-based Nonlinearity Evaluation (JNE), an interpretability metric for nonlinear neural encoding models. JNE quantifies nonlinearity by statistically measuring the dispersion of local linear mappings (Jacobians) from model representations to predicted BOLD responses, thereby approximating the nonlinearity of BOLD signals. Centered on proposing JNE as a novel interpretability metric, we validated its effectiveness through controlled simulation experiments on various activation functions and network architectures, and further verified it on real fMRI data, demonstrating a hierarchical progression of nonlinear characteristics from primary to higher-order visual cortices, consistent with established cortical organization. We further extended JNE with Sample-Specificity (JNE-SS), revealing stimulus-selective nonlinear response patterns in functionally specialized brain regions. As the first interpretability metric for quantifying nonlinear responses, JNE provides new insights into brain information processing. Code available at https://github.com/Gaitxh/JNE.

Jacobian-Based Interpretation of Nonlinear Neural Encoding Model
Xiaohui Gao!," Haoran Yang i Yue Cheng Mengfei Zuo', Yiheng Liu', Peiyang Li
Xintao Hu
Northwestern Polytechnical University Beejjig Jiaotong University
*Chongqing University Posts and Telccommunications
gaitxh@foxmail com, yhr06028mail nwpu. edu. en, yuecheng0bjtu. edu. cn, zuomengfei mwpu. edu. cn, Liuyiheng123@ma .wpu edu, cn,
mmp.
Lipeiyangccqupt.oo edu. cn, xhu@nwpu. edu. cn Equual contribution *Corresponding author
Abstract
In recent years, the alignment between artificial neural network (ANN) embed- dings and blood oxygenation level dependent (BOLD) responses in functional magnetic resonance imaging (IMRI) via neural encoding models has significantly advanced research on neural representation mechanisms and interpretability in the brain. However, these approaches remain limited in characterizing the brain's ss inherently nonlinear response properties. To address this, we propose the Jacobian- based Nonlinearity Evaluation (JÑE), an interpretability metric for nonlinear neu- ral encoding models. JNE quantifies nonlinearity by statistically measuring the dispersion of local linear mappings (Jacobians) from model representations to predicted BOLD responses, thereby approximating the nonlinearity of BOLD signals. Centered on proposing JNE as a novel interpretability metric, we vali- dated its effectiveness through controlled simulation experiments various tivation functions and network architectures, and further verified it on real fMRI data, demonstrating hierarchical progression of nonlinear characteristies from primary to higher-order visual cortices, consistent with established cortical orga- nization. We further extended JNE with Sample-Specifice stimulus-selective nonlinear response patterns in functionally specialized brain regions. As the first interpretability metric for quantifying nonlinear responses, JNE provides new insights into brain information processing. Code available at https://gg /gi thub. com/co////ooooc. Comm/aaitai GaaithhJJE...
(JNE-SS), revealing
Introduction
In recent years, researchers have made significant progress in constructing interpretable computational models of the brain by aligning artificial neural network (ANN) representations with functional magnetic resonance imaging (fMRI) signals [1, 2, 4. 5, 6]. These studies have typically been based on the neural encoding model framework, wherein hierarchical embeddings from ANNs are used input features predict voxel-level neural responses (i.e., blood oxygenation level dependent responses in this study) using either linear models (e.g., ridge regression) 7, 8, 9] or nonlinear regression models (c.g., deep neural 110 11. 5. 12], and the prediction performance is evaluated by metrics such as the coefficient of determination (fi²) or Pearson correlation coefficient (r)[11 13, 14] (see Fig. la). Based this framework, substantial body of work has shown that ANN-derived embeddings correspond well with the spatial distribution of functional brain regions in language [3, 4, 15. 16], vision [2. 17], and multimodal tasks [8, 18, 19]. With the advancement of neuroscience research, it has been observed that the brain's 's responses to external stimuli exhibit
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Feature space
Stimull
Linear mapping
Activity space
Encoder X companison
Brain nonincar soore (Ours)
Non-linear
feature mapping
or Non-linear mapping
Vini
JNC
Figure Comparison of brain encoding frameworks and nonlinearity analysis approaches. Neural encoding aligns ANN representations with fMRI BOLD responses via linear/nonlinear models to predict voxel responses. Nonlinearity estimation via R²² ratios between lincar/nonllneaa models fails when deep features embed nonlinearities, yielding comparable R²² and poor regional discriminability. c. JNE characterizes BOLD nonlinearity by statistically measuring Jacobian dispersion across inputs, showing higher discriminability across regions.
marked nonlinear characteristics. Numerous studies have demonstrated the nonlinear nature of neural responses [10. 11, 20, 21, 22, 23, 24, 25, 26]. However, key questions such as how to quantitatively characterize the degree of blood oxygenation level dependent (BOLD) response nonlinearity and how it is distributed across the cortex remain insufficiently addressed by systematic studies and theoretical frameworks.
Current neural encoding models face fundamental limitations in terms of nonlinear interpretability: linear neural encoding models can only establish simple mappings and fail to capture higher-order nonlinear features, while nonlinear neural encoding models are difficult to interpret due to their inherent "black-box" nature [11]. Therefore, quantitatively interpreting the nonlinear characteristics of neural encoding models has become a critical challenge. A common strategy is to infer the degree of nonlinearity in a given brain region by comparing the prediction performance of linear and nonlinear neural encoding models in that region [27, 28, 29]. However, this strategy essentially reflects only performance differences between models and fails to accurately reveal the intrinsic nature of BOLD response nonlinearity (see Fig.1b l and Section 3.2).
To address the above issues, we propose a novel interpretability metric for nonlinearity- Jacobian- based Nonlincarity Evaluation (JNE) to quantify the nonlinear characteristics of neural encoding models, thereby approximating voxel-level BOLD nonlinearity. The core idea is to calculate the local linear mappings (i.e., local derivatives, also known as Jacobian matrices) of the neural encoding model under different input stimuli, and to evaluate the level of nonlinearity by statistically measuring their variability across inputs (theoretical description in Section 1) Intuitively, an ideal linear neural encoding model should maintain a constant mapping across different inputs, resulting in a JNE value of 0; in contrast, if the model exhibits significant nonlinearity, the Jacobian matrices will vary across inputs, yielding a higher JNE value. Centering on JNE, we validate via simulations for nonlinearity guantification, then apply to fMRI for approximating visual cortex BOLD nonlinearity. JNE serves as a conditional interpretability metric for approximating voxel-level BOLD nonlincarity, rather than direct neural measure, due to BOLD confounds. In summary, the main contributions of this study include:
We propose JNE as novel interpretability metric for nonlinear neural encoding models, which quantifies nonlinearity by statistically measuring the dispersion of input-output Jacobian matrices (Section 2.4), and theoretical derivation indicates that JNE can serve an approximate measure of BOLD response nonlinearity (Section A.8). Controlled simulation experiments validate the effectiveness of JNE in quantifying the nonlinearity of output responses (Sections? 3.1 and A.7);
Our results demonstrate that inferring nonlinear properties by comparing R² values be- tween linear and nonlinear neural encoding models is insufficient to reveal the brain's 's nonlinear responses, especially when deep visual representations already embed nonlinear transformations (Section 3.2);
We find that the primary visual cortex tends to exhibit more linear responses, while higher- order visual cortices show stronger nonlinear characteristics under natural visual stimulation (Section 3.3). Overall, a hierarchical structure emerges with increasing nonlinearity from the primary to intermediate to higher-order visual cortex (Section 3.4);
We further define JNE-SS (Section A.10), which characterizes the nonlinear properties of individual samples at specific voxels. Our results indicate that BOLD response nonlinearity exhibits sample-selectiv preferences (Section 3.5).

---

2 Methods
2.1 IMRI Dataset Description and Feature Extraction
In this study, we used the Natural Scenes Dataset (NSD) [30], which provides whole-brain 7T fMRI responses from multiple subjects while viewing natural scene images. We focused four subjects (S1. S2, S5, S7) who completed the full experimental protocol as the primary participants in our analysis. Each image and its corresponding averaged beta map were treated dataset was split into training, validation, and testing sets with a ratio of 8:1:1 (8000:1000:1000 samples), ensuring
single sample. The
subject overlap between the training and validation sets, while the test set included repeated samples across subjects. We employed the pre-trained CLIP-ViT [31]' image encoder obtaín computational representations of the visual stimuli (Fig. 2a). Model performance was evaluated using N², and statistical significance was assessed through bootstrap resampling and FDR correction identify statistically significant activated voxels (see Section A.2 for details).
2.2 Linear and Nonlinear Neural Encoding Models
The general neural encoding model can be formulated as follows:
ÿ = fxx
(I)
where x denotes the feature space of external stimuli, ÿ is the predicted brain activity space (i.e., BOLD responses in this study), and f.) denotes the function of neural encoding models. Specifically, the computational representation of visual images spans the feature space, and the preprocessed beta maps serve as the brain activity space (y). We constructed both linear and nonlinear neural encoding models in simplified form. Each model consisted of four fully connected layers with skip connections, where the nonlinear neural encoding model additionally incorporated ReLU activation functions between layers (as shown in Fig.2b, and more details in Section A.3).
2.3 Computation of the Jacobian Matrix
The Jacobian matrix characterizes the local linear mapping between the representations (x) and the predicted neural responses (ý)) reflects the sensitivity of the neural encoding model to variations in the input. Upon completion of model training, we fixed the parameters of the neural encoding model and input the test set data to obtain the corresponding predictions. Given the i-th sample 2, R in the test set and the prediction of the nonlinear neural neural encoding model ÿ f(x,) RM, the Jacobian matrix JM, can be calculated by taking the derivative of the model's 's output ÿs with respect to its input v; as follows:
000L
80..
JM;- = Oij Ox,
T
D0,,1
E R
(2)
89LL
89..
D2.D
We denote the collection of Jacobian matrices of testing samples as JM E RNNDDM where D = 768 is the dimensionality of the representation in CLIP-VIT, N = 1000 is the sample size of the test set, and M is the number of voxels the brain. Notably, we adopt a forward-mode analytical differentiation approach (Section A.4) to compute Jacobians, accelerating the process by over 99.9% compared to traditional backward-mode automatic differentiation.
2.4 Jacobian-based Nonlinearity Evaluation Index
To systematically quantify the nonlinearity characteristics, this study proposed the Jacobian based Nonlinearity Evaluation index (JNE). This metric established deviation measure from the linear superposition principle in neural response systems by analyzing inter-sample statistical properties of Jacobian matrices. The computational workflow is defined as follows:
Step 1: Jacobian matrix centralization
JM = J.JM] N
JM, E E PXM
(3)
'tttps:/hhuggingface.com http://hhggggaaooeoo///////////s patch16
3

---

Regresentaatinss extracsin
Only 0093 during model branirg
N imeges
Inputs
Encoding models (neural network (T)
Outpuls
IMRI ditta
Linear neural encoding model
Residual connection
CLIP-VIT i
N-M
ANN structure
Oulputs
N-M
a
F
JM Bÿ /xx
Nonlinear neural encoding model
JM.
JNE
CRN
Inputs
connection
6dJM
AJM - JM JM J.
i
FC
PC2
RaLU V
FC4
Outouls
JMras
R.JIM
RaLU
Rai
Figure 2: Framework for representation extraction, neural encoding model, and JNE Compu- tation. a. Extract ANN representations from images using the CLIP image encoder. Use the extracted representations in (linear/nonlinear) neural encoding model to predict brain responses to images. c. Compute the Jacobian matrix to represent the mapping relationship between inputs and outputs of the neural encoding model. Further, calculate the mean, sum, and standard deviation of the Jacobian matrix to obtain the JNE metric.
This operation constructed baseline linear mapping matrix by eliminating stimulus-specific idiosyn- crasies through sample averaging.
Step 2: Sample-specific deviation computation
sAM, = JM, - JMmeanlll R.M
Manhattan norm (I- contraction along the feature dimension (D-axis) yields an LI deviation vector in voxel space for each sample. This metric quantified the absolute deviation between local linear approximations (per-sample Jacobians) and the global averaged mapping. Empirical tests in Section A.9.2 show that both LI and L2 norms yield highly consistent JNE spatial patterns, but adopt L1 for its robustness in nigh-dimensional settings [32, 33, 34).
Step 3: Mean absolute deviation estimation
1
Ap E[sAJM
Σ sAJM, RRM N
(5)
This calculated the cross-sample averaged absolute deviation, establishing baseline reference for fluctuation intensity of mapping relationships at the voxel level.
Step 4: Nonlinear dispersion quantification
JNEm asAJMM m)
(sJJM,, Apm) Vm (1,...,
(6)
By computing voxel-wise standard deviations (a(-),, this precisely characterized the dispersion of mapping relationships across stimuli. This statistic intrinsically reflects eccond-moment features of Jacobian matrices in sample space, with higher-order nonlinear systems exhibiting significant o-value growth. Thus, JNE quantifies nonlinearityvia Jacobian dispersion statistics, approximating voxel-level BOLD nonlinearity. The complete flow of JNE calculation is shown in Fig.2c.
Mathematically, the JNE metric rigorously corresponds to the quantitative evaluation of the homo- geneity and superposition principles in linear systems:
an ideal linear neural encoding model, where JM, JMmean JNEm-
Nonlinear mechanisms induce cally increasing with nonlinearity strength.
Jacobian shifts, with dispersion JNE,, monotoni-
Critically, both input features and output beta maps were z-score normalized prior model training and JNE computation. This standardization ensures that the input space, output space, and mapping space operate under consistent scale assumptions, enabling meaningful comparison Jacobian variability across voxels and brain regions. Therefore, the proposed metric inherently circumvents

---

Falu
Relli
GELU
Input
+ Oulput
Activation function
inl Suaeal and behwetee
ReLU
Lesky ReLU
64
GELU
Swish
a
Input X
Leaky
cornaction
015
2.17
A12
514
214
0.2
Loc1
1
Imput
Output
Loc2
0
0
Combination
2
50
0.0
C
Combination
Figure Validation of the JNE metric in controlled simulation experiments. a. Simulated neuron model; b. Comparison between the normalized JNE (min-max normalization) and the response curves of the normalized second-order derivatives (absolute value and min-max normalization) of various activation functions(Respong curves, 1st and 2nd derivatives of the four activation functions A.5); e. Schematic of the ANN-based simulation framework, where the model consists three-layer fully connected residual network. Locl and Loc2 represent pluggable activation function modules that serve as controllable sources nonlinearity; d. Binary configurations of activation function placement (0 indicates absence, indicates presence); e. Comparison of JNE values across four activation functions under different activation placement configurations (meantstd).
cross-regional scale normalization- magnitude of standard deviation directly encodes biological sensitivity differences in neural responses. Stimulus-sensitive brain regions (c.g., carly visual cortex) naturally exhibit greater Jacobian variability across inputs, which mathematically aligns with the sample-dependent characteristics of nonlinear mechanisms. This approach preserves neuro- physiological interpretability while providing a unified measurement benchmark for cross regionalgional nonlinearity comparisons.
3 Results
3.1 Simulation-Based Validation of JNE
Single-neuron-level validation of activation functions. To validate the ability of JNE to accurately characterize system-level nonlinearity, we first constructed simple and controllable simulation framework. Specifically, focused on four commonly used activation functions-ReLL, Leaky ReLU, GELU, and Swish-and designed a minimal neuron model consisting solely of input, nonlinear transformation, and output (Fig.3a). JNE essentially measures the "statistical dispersion" of the local linear mappings (i.e., Jacobians) of neuron's 's response to input perturbations across samples. Given an activation function o(x). its first derivative ø' (x(x) corresponds to the Jacobian matrix. If o' (x) varies significantly within a local neighborhood [x - 6, x + 6]. it indicates that the neuron's 's response exhibits non-un unifforo shifts with input changes thus violating linear assumptions. When this rate of change can be described by the second derivative (z)). Therefore, in theory, the second derivative of (xx) can serve as a structural reference for JNE: larger ('(x) implies stronger local nonlinearity, which should be reflected by higher JNE value. In our experiment, inputs were defined over the range [-10, 10), and a sliding window [a: 6, + 8) (with 0.001 and step size of 0.02) was used. For each central point v, we sampled 1001 perturbed inputs within the local window and computed the sequence of first derivatives of the activation function to calculate the JNE value. This process repeated over the entire input range to produce JNE curves, which were then compared to the second derivative curves of the corresponding activation functions.
As shown in Fbg.3b, the JNE curves closely resemble the structural patterns of the second derivatives for all activation functions. Smooth functions like GELU and Swish exhibited high JNE values over broad input ranges, whereas ReLU only produced peak around its discontinuity at 0, with a value of zero elsewhere. Additionally, JNE was strictly zero linear regions of the activation functions, validating its theoretical soundness. However, it is worth noting that JNE measures the dispersion of local linear mappings across input samples, which correlates with the magnitude of second derivatives in theory, but is not identical due to its statistical nature.
5

---

Linear encoder
Nonlinear encoder
Mecoi perfere Boroangee UPO
Proddan por PAEE nereted
Figure 4: Comparison of prediction performance between linear and nonlinear neural encoding model. a. Prediction performance (R²) the linear neural encoding model. b. R of the nonlinear encoder. c. 2D histogram of linear and nonlinear neural encoding model performance relative to the noise ceiling (calculation [30]) of significantly activated voxels and the 80% noise ceiling. For visualization purposes, only voxels with prediction performance significantly higher than chance (P<0.05, FDR corrected, one oiddsided t-test) are displayed in the brain maps.
0.82
0.85
0.93
0.95
0.99
Nontr werr enrodder per performnno (Re)
Brain, (=0.98
PPA
RSC
OPA
FFA-1
FFA-2
0.99
0.99
0.98
r= 0.99
1.00
3.
36 1.
13
03
81
10
Linear enooder perfomaane (84)
Figure 5: Scatter plot comparison between linear and nonlinear neural encoding models. Scatter plots comparing the R2 of linear and nonlinear neural encoding models across the whole brain and 10 regions of interest (ROIs) defined by the HCP-MMP atlas [35]. The Pearson correlation coefficient between the R² sequences of the two models is also reported.
Network-level validation via ANN-based simulation. To further assess JNE's practical utility in multilayer network structures, we developed an ANN simulation framework (Fig.3c). The model consisted a three-layer fully connected residual architecture and was not subject to any training was used solely for forward propagation to obtain the Jacobian matrices between inputs and outputs. Nonlinearity was introduced via two plug-in locations (Loc1 and Loc2), where activation functions could be flexibly inserted in fully controlled manner (Fig.3c and d). This setup allowed to systematically evaluate JNE's sensitivity to different nonlinear configurations without altering the overall network architecture. To simulate the real neural encoding models used in this study, we set the input dimension to 768 (mimicking CLIP- VIT image representations) and generated an input matrix of size 1000 768, with the corresponding output being simulated response vector of size 1000 1. For each input sample, we computed the Jacobian of the output with respect to the input and evaluated the overall nonlinearity using the standard JNE pipeline. We conducted comparative analyses for layers).
four activation functions under different Loc 1/cl/2 configurations (none/shallow/deepholl Each configuration was sampled 200 times ensure stable statistical estimates (see Section A.6 for details).
As illustrated in Fig.3e, the JNE value was under the ideal linear system (no activation function), consistent with thcoretical expectations. Upon the introduction activation functions, JNE values increased significantly. Interestingly, we observed that the position of the nonlinearity affected the results: configurations with activations in the deeper layer (Loc2) generally yielded higher JNE values compared to those in the shallow layer (Locl), 1), and the two-layer configuration led to the strongest nonlinearity. This trend was consistent across all activation functions, revealing an pattern.
These results demonstrate that JNE not only reconstructs the intrinsic nonlinear structures of activation functions at the theoretical level but also reliably detects and quantifies diverse sources of nonlinearity

---

0.1
JNE
0.5
a
(All colored vorels 0.05, FDR corrected)
0.30
-0.35
PVC
IVC
HVC
0.25
0.20
JNE
0.55
6.10
Brain
V2
19
V4
EBA
PPA
RSC
OPA
FRA-2
Urain
V1
V2
V3
V4
EBA
PPA
RSC
OPA
FAA-1
ITA-2
b
Figure 6: Nonlinear distribution in the brain. a. Activation map of JNE across the brain, illustrating the spatial distribution of model-predicted BOLD response nonlinearity. b. Mean JNE values across the whole brain and 10 visual ROIs (mean=std). c. JNE distribution across the whole brain and 10 visual ROIs, with the x-axis indicating density.
in complex wwwworks-validatinn its effectiveness as a nonlinearity metric, paving the way for approximating BOLD nonlinearity in fMRI.
3.2 Nonlinearity Analysis Based on Neural Encoding Models Comparison
As response the introduction, this paper begins by addressing central question: Can infer the nonlinearity individual voxels by comparing the R2 of linear versus nonlinear neural encoding models? As first step, we extracted embedding features from the final layer of CLIP-VIT and aligned them to fMRI responses using both a linear and nonlinear neural encoding model, respectively. Fig.4 presents the relevant results for Subject (all results in the main text refer to Subject 1; other subjects in Section B).
As shown in Figs.4a and 4b, both types of models exhibit highly consistent predictive performance across the whole brain, with activations primarily concentrating the primary, intermediate, and hidher-order visual cortices, well as parts of the prefrontal cortex. Furthermore, we observed that both models achieved predictive performance close to the noise ceiling for many voxels (Fig.4c), with some voxels even exceeding 80% of the noise ceiling. This indicates that both linear and nonlinear neural encoding models can effectively account for voxel responses in the held-out test data (explaining up to 74% of the response variance). To further assess the consistency of the two models in predictive performance, we conducted a correlation analysis between the R sequences of the linear and nonlinear neural encoding models across the whole brain and within ten predefined vision-related ROIs (Fig.5). Results show strong positive correlations between the R² sequences of the models in all ROIs, with the lowest Pearson correlation observed in VI (min 0.82), indicating a strong relationship between the activation patterns of the two neural encoding models.
Together, these findings support the notion that when the deep features of visual encoder have already sufficiently captured the nonlinear transformations of visual information, a linear neural encoding model can still effectively fit neural response patterns (Fig. lb, demonstrating the high similarity between linear and nonlinear neural encoding models). In other words, the difference in predictive performance between linear and nonlinear neural encoding models becomes marginal, making unreliable to infer nonlinearity solely based on model performance comparisons. To address this limitation, we proposed JNE, concise and principled metric for characterizing nonlinearity. In the following section, we apply INE to fMRI data as real-world validation.
3.3 Spatial Distribution of BOLD Response Nonlinearity in Visual Cortex
Fig.6n 6a illustrates the spatial distribution of JNE across the cortical surface. We observed that the nonlinearity of neural responses was not uniformly distributed across the brain but instead

---

PVC
IVC
PVC
PVD
NC
HWC
<<<< <
b
Figure 7: Hierarchical analysis of nonlinear characteristics in visual cortex. a. Bar plot com- parison of mean JNE values across primary, intermediate, and higher-order visual cortices; Distributional differences in JNE across the three cortical levels; Significance matrix of pairwise JNE differences among primary, intermediate, and higher-order visual cortices (one-sided test, P 0.05), where entry (i. indicates that cortex i exhibits significantly higher JNE than cortex d. Colors
Region-level visualization hierarchical nonlinearity based on the significance matrix represent cortical levels (primary = black, intermediate = red, higher-order yellow), not voxel-wise specificity. This visualization validates the consistent hierarchical trend across the visual cortex.
concentrated in higher-order cortical parietal-occipital junction (TPOJ), and parts of the prefrontal cortex. Figs.6b-c and Fig. 17 further quantify the statistical characteristics of JNE across the whole brain and within ten predefined ROIs. The Primary Visual Cortex (PVC-V1) exhibited the lowest mean JNE values, with distribution tightly concentrated in the low-value range, suggesting that neural responses this region are more consistent with linear characteristics. In Intermediate Visual Cortex (IVC-V2, V3, V4), JNE values showed a moderate increase and broader spread, indicating the presence of moderate nonlinear response components. In contrast, Higher-order Visual Cortex (HVC-EBA, PPA, RSC, OPA, FFA-1 and FFA-2) demonstrated the highest JNE values, with distributions skewed toward higher values. This pattern suggests that these regions exhibit more pronounced nonlinear response properties, approximating hierarchical BOLD patterns. Collectively, these results indicate that brain regions exhibiting stronger JNE values, approximating hierarchical BOLD nonlinearity higher-order cortices, which are typically involved in advanced information integration and abstract semantic processing.
higher-order visual cortices, the temporo-
3.4 Hierarchical Progression of Nonlinearity Across the Visual Cortex
Furthermore, we observed a hierarchical progression of INE-approximated BOLD nonlinearity across visual cortex. As shown in Fag.7a, bar plots comparing JNE across primary, intermediate, and higher-order visual cortices reveal clear trend: the number of voxels exhibiting high JNE values increases with cortical hierarchy. Fig.7h 7b further illustrates the distribution JNE values within each visual cortical tier. In PVC, JNE values are predominantly concentrated in the lower range. In IVC, the distribution broadens with slight shift toward higher values. In contrast, HVC displays significantly right-skewed distribution with pronounced long-tail pattern, indicating stronger nonlinear response characteristics. Figs.7c 7c and 7d present heatmaps and statistical comparisons that highlight the significant differences in JNE across the three cortical levels. These results reveal robust, layer-wise enhancement of nonlinearity across the visual hierarchy-ize i.e., primary < intermediate < higher-order visual cortex. These findings suggest that JNE, as applied to the neural encoding model, uncovers progressive nonlinear hierarchy in the predicted BOLD responses across the visual system, beginning in the primary cortex, transitioning through intermediate areas, and peaking higher-level regions (d)..7d). This hierarchical organization of BOLD response nonlinearity aligns well with previous neurophysiological evidence supporting functional stratification across the visual cortex [36, 37, 38).
3.5 Sample Selectivity of BOLD Response Nonlinearity
To characterize the sample-specific nonlinear responses of different brain regions to visual stimuli, we further introduced the JNE with Sample-Specifiche (JNE-SSE NR theoretical details in Section A.10). In brief, JNE-SS quantifies the degree of nonlincarity in a voxel's response to individual image

---

Lab4
28
PPA
Rec
OPA
ITA-1
Figure 8: Sample-specific analysis of brain nonlinear encoding. a. The test set of 1,000 images was divided into 43 categories and sorted by category; the top images in each category for display (see detailed results in Fig. 12). b. MoM was computed across the whole brain and within 10 predefined ROIs. c. The results in (b) were normalized within each ROI using min-max normalization to reflect category preferences of each ROI. Note that the heatmap in Fig. 8c applies min-max normalization within each ROI, which may attenuate visual differences in absolute MoM values, For quantitative comparisons across regions, we refer to Table 5 in Section A.12, which reports MoM values without applying min-max normalization.
selected
stimuli: a lower JNE-SS value for voxel-image pair indicates that the voxel exhibits a more linear response preference for that specific stimulus, whereas a higher value implies stronger nonlinear tuning. For practical analysis, we first applied t-SNE for 2D embedding of all image stimuli and then performed K K-means clustering to categorize them into 43 clusters. This clustering structure preserved a meaningful continuum along the dimension of semantic complexity (Fig, Sa; more details in Section A.1 Subsequently, we computed the Mean of JNE-SS (MoJ E R43x(M for each voxel across image clusters, and further derived the Mean of MoJ (MoM E 1133111 at the whole-brain and ROI levels (including the whole brain and 10 ROIS). This metric reflects the variation in BOLD response nonlinearity elicited by different image categories across brain regions.
The result reveals that MoM values are generally lower in primary and intermediate visual cortices, while notably higher in higher-order visual areas--articuulalyy the EBA (Fig. 8b, Fig. 21a), cor- roborating our previous findings in Fig. 6. Additionally, EBA maintains consistently high JNE-SS values across multiple image categories, suggesting potentially central role of EBA in mediating nonlinear visual processing. Further examination of regional nonlinearity preferences for specific image types revealed additional insights (Fig. 8c. Fig. 21b). Images depicting outdoor scenes or complex backgrounds evoked stronger nonlinear responses in PPA and RSC, while stimuli involving people or faces induced pronounced nonlinear fluctuations with the known functional specificity of these regions [2, 30, 39, 40, 41]. Interestingly, we also observed broad nonlinear responses in primary and intermediate visual areas across a wide range of image categories (Fig.21b), contrasting with the more selective tuning seen in higher-order visual regions.
FFA and EBA. These findings align
4 Conclusion and Discussion
In this work, we introduce JNE, a novel interpretability metric for nonlinear models via Jacobian dispersion, approximating voxel BOLD nonlinearity. Through rigorous simulation experiments and its application to iMRI data from natural scene viewing, we validate JNE's effectiveness in discerning hicrarchical gradients and stimulus-selective patterns of nonlinearity across the visual cortex. Importantly. JNE evaluates input-dependent variability in the model's 's local linear mappings to BOLD signals, serving as a diagnostic tool for encoding model interpretability rather than direct assay of intrinsic neural nonlinearity.
The metric deliberately retains scale information in its lacobian-based standard deviation, as variations predicted response amplitudes may encode physiologically relevant differences in neural gain. This preserves neuroscientific utility, despite the inherent coupling of gain and nonlinearity. Although JNE could principle confound response gain, nonlinearity degree, and input dimensionality, theoretical derivations (Section A.9.1) and empirical assessments real data (Section A.9.2) demonstrate that

---

these confounds are minimal under the constrained architectures and stimulus conditions typical of neural encoding frameworks.
In real data experiments, our study demonstrates that the nonlinear characteristics of neural responses exhibit significant differences across brain regions (Section 3.3). Specifically, JNE values are generally lower in primary cortices and significantly higher in higher-order cortices. The elevated JNE values in higher-order regions suggest that the encoding model must employ more nonlinear mappings to predict BOLD responses in these areas, which is consistent with the hypothesis that these areas engage in more nonlinear processing of visual information. However, JNE captures the net nonlinearity of the input-to-BOle pathway, which integrates both neural and hemodynamic factors.
In addition, we observed hierarchical trend of BOLD response nonlinearity within the visual cortex (Section 3.4). This result reveals a gradient pattern consistent with the hierarchical organization of the visual processing pathways [36, 37, 38]. Specifically, PSV (V1) is primarily responsible for processing low-level perceptual features such as edges, orientation, and spatial frequency [30, 35). The low JNE values in this region indicate that it may perform only simple linear processing of visual information. Subsequently, the processed information is transmitted along the dorsal and ventral visual pathways IVC and HVC, during which JNE values progressively increase. The observed monotonic increase in JNE from primary to higher-order visual cortices supports a hierarchical progression of information processing along the ventral visual stream, consistent with established models of visual cortical organization. [2. 35, 42, 43, 44]. We note that fMRI responses integrate both feedforward and feedback signals, and thus JNE captures the net nonlinearity of the system under natural viewing conditions. Additionally, it suggests that the propagation of visual information in the brain may follow a progressive pattern from now-noniinearity (low JNE) regions to high-nonlinearity (high JNE) regions, gradually transitioning from linear to nonlinear processing.
Section 3.5 further supports this view, as reflected in the relatively low MoM values in PVC, contrasted with the higher MoM values in higher visual cortices (such as PPA, RSC, FFA), which also exhibit stronger selectivity different image categories. Moreover, we found that PPA and RSC showed stronger nonlinear responses to outdoor scenes, whereas FFA displayed highly significant nonlinear responses to faces and human figures. These results corroborate the functional selectivity of these regions to different stimuli [2. 30, 45], and also suggest that the higher nonlinearity observed in these areas may underlie their roles in supporting semantic generalization, contextual integration, and complex perceptual discrimination core aspects of high-level visual cognition. Notably, we also found that even in V1-V4, MoM values increased when images contained rich backgrounds, multiple objects, or higher visual complexity, revealing a limited yet genuine capacity for nonlinear modulation in early visual processing stages.
Generally, the primary contribution of this work is the proposal and validation of the JNE metric. While we demonstrate its application in revealing nonlinearity in visual cortex, these results should be seen as preliminary validations of the method's ' feasibility.
Limitations and Prospects
Despite its contributions, this study has several limitations that suggest avenues for future work. First, INE approximately quantifies nonlinearity in mappings to BOLD signals rather than directly to neural activity. BOLD is indirect hemodynamic measure confounded by vascular factors (e.g.. blood flow refractory effects and HRF nonlinearity), which cannot be fully separated from neural signals in standard fMRI. Thus, JNE reflects the composite nonlinearity of the Input-to-BOI. pathway, limiting direct neuronal inferences. While simulations used a controllable ANN framework, we did not exhaustively assess how activation function count, placement, or type influences overall nonlinearity, nor track its evolution during training. Furthermore, the computation of JNE essentially depends on the specific encoding model architecture and stimulus distribution used. Future studies should systematically evaluate the consistency of JNE across different model architectures (e.g., CNNs, Transformers), training objectives (supervised,
superviedd, unsupervised), and stimulus sets, thereby validating its robustness as comparative tool. Additionally, no absolute JNE threshold was defined to classify voxels as "linear" or "nonlinear," " restricting analyses to relative comparisons and voxel-level conclusions. Moreover, our sample-seeciiic analysis relied on informal t-SNE ordering for image categories, lacking rigorous standards. Finer-graines feature/patch-leee investigations are needed to pinpoint what drives region-specific nonlinear responses.
10

---

6 Acknowledgements
This work was supported by the National Natural Science Foundation of China under Grants 62076205 and 62576065.
References
[1] Haibao Wang, Lijie Huang, Changde Du, Dan Li, Bo Wang, and Huiguang He. Neural encoding for human visual cortex with deep neural networks learning "what" and "where". IEEE Transactions on Cognitive and Developmental Systems, 13(4):827-840, 2020.
[2] Aria Y Wang, Kendrick Kay, Thomas Naselaris, Michael J Tarr, and Leila Wehbe. Better models of human high-level visual cortex emerge from natural language supervision with a large and diverse dataset. Nature Machine Intelligence, 1415):1415- 1426, 2023.
[3] Xu Liu, Mengyue Zhou, Gaosheng Shi. Yu Du, Lin Zhao, Zihao Wu, David Liu, Tianming Liu, and Xintao Hu. Coupling artificial neurons in bert and biological neurons in the human brain. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 37, pages 8888-8896, 2023.
[4] Mengyue Zhou, Xu Liu, David Liu, Zihao Wu, Zhengliang Liu, Lin Zhao, Dajiang Zhu, Lei Guo, Junwei Han, Tianming Liu, et al. Fine-grained artificial neurons in traff-rransformess for disentangling neural auditory encoding. In Findings of the Association for Computational Linguistics: ACL 2023, pages 7943-7956, 2023.
[5] Xiaohui Gao, Yue Cheng, Peiyang Li, Yijie Niu, Yifan Ren, Yiheng Liu, Haiyang Sun, Zhuoyi Li, Weiwei Xing, and Xintao Hu. Linbridge: A learnable framework for interpreting nonlinear neural encoding models. arXiv preprint arXiv:2410.20053, 2024.
[6] Zitong Lu. Visualizing the mind's eye: a future perspective applications of image reconstruc- tion from brain signals to psychiatry. Psychoradiology, 3:kkad022, 2023.
[7] Alexander G Huth, Wendy A De Heer, Thomas L Griffiths, Frédéric E Theunissen, and Jack L Gallant. Natural speech reveals the semantic maps that tile human cerebral cortex. Nature, 532(7600):453 458, 2016.
[8] Jerry Tang, Meng Du, Vy A Vo, Vasudev Lal, and Alexander G Huth. Brain encoding models based multimodal transformers transfer across language and vision. arXiv preprint 2305v /2244 2023.
[9] Jerry Tang, Amanda LeBel, Shailee Jain, and Alexander G Huth. Semantic reconstruction of continuous language from non-invasive brain recordings. Nature Neuroscience, 26(5):858-866, 2023.
[10] Danny Dongyeop Han, Yunju Cho, Jiook Cha, and Jay-Yoon Lee. Mind the gap: Aligning the brain with language models requires nonlinear and multimodal approach. arXiv preprint 1rXiv:2502.12771, 2025.
[11] Thomas Naselaris, Kendrick Kay, Shinji Nishimoto, and Jack Gallant. Encoding and decoding in fmri Neuroimage, 56(2):400 410, 2011.
[12] Yuanning Li, Huzheng Yang. and Shi Gu. Enhancing neural encoding models for naturalistic perception with multi-level integration of deep neural networks and cortical networks. Science Bulletin, 17311))7738-111
2024.
[13] Subba Reddy Oota, Zijiao Chen, Manish Gupta, Raju S Bapi, Gael Jobard, Frédéric Alexandre, and Xavier Hinaut. Deep neural networks and brain alignment: Brain encoding and decoding (survey). arXiv preprint arXiv:2307. 2023.
[14] Jian Gao. R-squared (r2)-how much variation explained? Research Methods in Medicine & Health Sciences, 1(4):104-109 2024. 2024.
11

---

[15] Sreejan Kumar, Theodore R Sumers, Takateru Yamakoshi, Ariel Goldstein, Uri Hasson, Ken- neth Norman, Thomas L Griffiths, Robert D Hawkins, and Samuel A Nastase. Shared functional specialization in transformer-based language models and the human brain. Nature communications, 15(1)55533, 2024.
[16] Ariel Goldstein, Zaid Zada, Eliav Buchnik, Mariano Schain, Amy Price, Bobbi Aubrey, Samuel A Nastase, Amir Feder, Dotan Emanuel, Alon Cohen, et al. Shared computational principles for language processing in humans and deep language models. Nature neuroscience, 25(3):369-380, 2022.
[17] Oliver Contier, Chris I Baker, and Martin N Hebart. Distributed representations of behaviour- derived object dimensions in the human visual system. Nature Human Behaviour, 8(11):2179- 2193, 2024.
[18] Sara F Popham, Alexander G Huth, Natalia Y Bilenko, Fatma Deniz, James S Gao, Anwar O Nunez-Elizalde, and Jack L Gallant. Visual and linguistic semantic representations are aligned at the border of human visual cortex. Nature neuroscience, 24(1 1):1628-1636, 2021.
[19] Jing Sui, Dongmei Zhi, and Vince D Calhoun. Data-driven multimodal fusion: approaches and applications in psychiatric research. Psychoradiology, 3:kkad026, 2023.
[20] Gottfried Kress. Non-linear mechanisms in the brain. Zeitschrift für Natwrforschung C, 53(7-8):677-685, 1998.
[21] D Vignesh, Shaobo He, and Santo Banerjee. review the complexities of brain activity: Insights from nonlinear dynamics in neuroscience. Nonlinear Dynamics, 113(5):451--552 2025.
[22] Betty Tuller, Noël Nguyen, Leonardo Lancia, and Gautam K Vallabha. Nonlinear dynamics in speech perception. Nonlinear Dynamics in Human Behavior, pages 135-150, 2011.
[23] Shailee Jain and Alexander Huth. Incorporating context into language encoding models for fmri. Advances in neural information processing systems, 31, 2018.
[24] Yibo Cui, Kai Qiao, Chi Zhang. Linyuan Wang, Bin Yan, and Li Tong. Gabornet visual encoding: A lightweight region-based visual encoding model with good expressiveness and biological interpretability. Frontiers in Neuroscience, 15:614182, 2021.
[25] Nanyin Zhang, Xiao-Hong Zhu, and Wei Chen. Investigating the source of bold nonlinearity in human visual cortex in response to paired visual stimuli. Neuroimage, 43(2):204-21,, 2008.
[26] Xiaohong Wan, Jorge Riera, Kazuki Iwata, Makoto Takahashi, Toshio Wakabayashi, and Ryuta Kawashima. The neural basis of the hemodynamic response nonlinearity in human primary visual cortex: Implications for neurovascular coupling mechanism. Neuroimage, 32(2):616-625, 2006.
[27] Chi Zhang, Kai Qiao, Linyuan Wang, Li Tong, Guoen Hu, Ru-Yuan Zhang, and Bin Yan. visual encoding model based on deep neural networks and transfer learning for brain activ- ity measured by functional magnetic resonance imaging. Journal of neuroscience methods, 325:108318, 2019.
[28] Jingwei Li, Panpan Chen, Chi Zhang, Linyuan Wang, Bin Yan, and Li Tong. An improved gabornet visual encoding model with nonlinear mapping and data augmentation. In 2022 7th International Conference on Computational Intelligence and Applications (ICCIA), pages 180- 184. IEEE, 2022.
[29] Yibo Cui, Chi Zhang, Kai Qiao, Linyuan Wang, Bin Yan, and Li Tong. Study on representation invariances of cnns and human visual information processing based on data augmentation. Brain Sciences, 1009):602, 2020.
[30] Emily J. Allen, Ghislain St-Y St-Yves, Yihan Wu, Jesse L. Breedlove, Jacob Prince, Logan Dowdle, Matthias Nau, Bradley Caron, Franco Pestilli, Ian Charest, J. Benjamin Hutchinson, Thomas Naselaris, and Kendrick Norris Kay. A massive 7t fmri dataset to bridge cognitive neuroscience and artificial intelligence. Nature Neuroscience, 25:116 126, 2021.
12

---

[31] Alec Radford, Jong Wook Kim, Chris Hallacy, Aditya Ramesh, Gabriel Goh, Sandhini Agarwal, Girish Sastry, Amanda Askell, Pamela Mishkin, Jack Clark, et al. Learning transferable visual models from natural language supervision. In International conference machine learning, PMLR, 2021.
pages
[32] Peiyang Li, Xiaohui Gao, Cunbo Li, Chanlin Yi, Weijie Huang, Yajing Si, Fali Li, Zehong Cao, Yin Tian, and Peng Xu. Granger causal inference based dual laplacian distribution and its application to mi-bei classification. IEEE Transactions on Neural Networks and Learning Systems, 2023.
[33] Xiaohui Gao, Zhengya Ma, Yue Cheng, Peiyang Li, Yifan Ren, Pengcheng Zhu, Xiaoxu Wang, and Xintao Hu. A novel residual-based bayesian pnpectaaioo-maxxmizatioo adaptive kalman filter with inaccurate and varyiig noise covariances. Measurement, 235:1 14937, 2024.
[34] Peiyang Li, Cunbo Li, Joyce Chelangat Bore, Yajing Si, Fali Li, Zehong Cao, Yangsong Zhang, Gang Wang, Zhijun Zhang, Dezhong Yao, et al. L1-norm based time-varying brain neural network and its application to dynamic analysis for motor imagery. Journal of Neural Engineering. 19(2):22019, 2022.
[35] Matthew F Glasser, Timothy S Coalson, Emma C Robinson, Carl D Hacker, John Harwell, Essa Yacoub, Kamil Ugurbil, Jesper Andersson, Christian F Beckmann, Mark Jenkinson, et al. A multi-modal parcellation human cerebral cortex. Nature, 536(7615):111-175 2016.
[36] Daniel L. K. Yamins, Ha Hong, Charles F. Cadieu, Ethan A. Solomon, Darren Seibert, and James DiCarlo. Performance-optimee hierarchical models predict neural responses in higher visual cortex. Proceedings of the National Academy of Sciences, 11 (11(())8619--88) 9--624, 2014.
[37] Kalanit Grill-Spector and Rafael Malach. The human visual cortex. Annu. Rev. Neurosci., 27(():649-677) 677, 2004.
[38] Victor AF Lamme, Hans Super, and Henk Spekreijse. Feedforward, horizontal, and feedback processing in the visual cortex. Current opinion in neurobiology, 8(4):529-5333 1998.
[39] Peter Donaldson, Nicole J Rinchart, and Peter G Enticott. Noninvasive stimulation of the temporoparietal junction: systematic review. Neuroscience & Biobehavioral Reviews, 04:547-572, 2015.
[40] Kathleen A Hansen, Carlton Chu, Annelise Dickinson, Brandon Pye, J Patrick Weller, and Leslie G Ungerleider. Spatial selectivity the temporoparietal junction, inferior frontal sulcus, and inferior parietal lobule. Journal of Vision, 15(13): 15, 2015.
[41] Kevin A Pelphrey, Juliana Lopez, and James P Morris. Developmental continuity and change in responses to social and nonsocial categories in human extrastriate visual cortex. Frontiers in human neuroscience, 3:885, 2009.
[42] Nadine Dijkstra, Peter Zeidman, Sasha Ondobaka, Marcel AJ van Gerven, and Karl Friston. Distinct top-down and bottom-up brain connectivity during visual perception and imagery. Scientific reports, 7(1):5677, 2017.
[43] Hristos S Courellis, Juri Minxha, Araceli R Cardenas, Daniel L. Kimmel, Chrystal M Reed, Taufik A Valiante, C Daniel Salzman, Adam N Mamelak, Stefano Fusi, and Ueli Rutishauser. Abstract representations emerge in human hippocampal neurons during inference. Nature, 632(8026):841-88))) 2024.
[44] Lijuan Wang, Ling Huang, Mengsha Li, Xiaotong Wang, Shiyu Wang, Yuefa Lin, and Xilin Zhang. An awareness-dependee mapping of saliency in the human visual system. Neurolmage, 247:1 18864, 2022.
[45] Libi Kliger and Galit Yovel. The functional organization of high-level visual cortex determines the representation of complex visual stimuli. Journal of Neuroscience, 40439):7545-7558 2020.
[46] Paul Richard Halmos. Fimete-dimensional vector spaces. Springer Science Business Media, 2012.
13

---

[47] Christopher M Bishop and Nasser M Nasrabadi. Pattern recognition and machine learning. volume 4. Springer, 2006.
[48] Tsung-Yi Yi Lin, Michael Maire, Serge Belongie, James Hays, Pietro Perona, Deva Ramanan, Piotr Dollár, and C Lawrence Zitnick. Microsoft coco: Common objects in context. In Computer Vision-ECCV 2014: 13th European Conference, Zurich, Switzerland, September 6-12, 2014, Proceedings, Part V 13, pages 740-755. Springer, 2014.
[49] Aria Y. Wang, Kendrick Kay, Thomas Naselaris, Michael Tarr, and Leila Wehbe. Incorporating natural language into vision models improves prediction and understanding of higher visual cortex. bioRxiv, 2022.
[50] Vighnesh Subramaniam, Colin Conwell, Christopher Wang, Gabriel Kreiman, Boris Katz, Ignacio Cases, and Andrei Barbu. Revealing vision-language integration in the brain with multimodal networks. arXiv preprint arXiv:2406. /4481, 2024.
[51] Heng Huang, Xintao Hu, Yu Zhao, Milad Makkie, Qinglin Dong, Shijie Zhao, Lei Guo, and Tiamming Liu. Modeling task fini data via deep convolutional autoencoder. IEEE transactions on medical imaging, 37(7):1551-1561 2017.
[52] Zijun Zhang. Improved adam optimizer for deep neural networks. In 2018 IEEE/ACM 26th international symposium quality of service (IWQoS), pages 1 1-2. leee, 2018.
[53] Laurens Van der Maaten and Geoffrey Hinton, Visualizing data using t-sne. Journal of machine learning research, 9(1 1), 2008.
14

---

A Appendix A
A.1 Theoretical Basis of Nonlinearity Measurement
In a liiteedimensional space, an ideal linear system (here the system refers to neural encoding model) should satisfy two fundamental principles: additivity (f(æ1 + x2) = f(x1) ) + f(x2)) and homogeneity (f(ar) = a.((x) [46, 47]. These two properties must hold simultaneously for the system to be formalized as = Wx, where W a mapping matrix that remains invariant across samples. A nonlinear system longer adheres this constant mapping. Its input-output relation f(x) can instead be locally approximated as yi = W, where W; denotes the local derivative (i.... the Jacobian matrix) corresponding to the i-th input sample. Since the system's 's response to input perturbations depends on the current input, W; varies across samples, thereby violating the input-invariance property of a linear system.
Based on this analysis, the present study constructs a quantitative measure of deviation from linearity by evaluating the statistical variation of Jacobian matrices across samples. Specifically, if W, remains consistent across samples, the system can be regarded as linear; conversely, greater variability of W; indicates stronger deviation from linearity and hence higher degree of nonlinearity in the neural response at that region or voxel.
A.2 [MRI Data Description, Preprocessing, and Statistical Significance Evaluation
This study uses the Natural Scenes Dataset (NSD) [30], which contains IMRI recordings from eight subjects passively viewing 73,000 colored natural images over period of more than 40 hours. The images are sourced from the MS-COCO dataset [48], with each image presented for 3 seconds and repeated three times across 50-40 scanning sessions. Each subject completed approximately 22,000 to 30,000 trials. The fMRI data were acquired using a whole-brain 7T gradient-echo EPI sequence with spatial resolution of 1.8 mm and TR of 1.6 seconds. Single-trial beta maps were estimated using customized general linear model (GI.M) and released alongside the raw fMRI data [49]. Following previous studies [49]. the beta maps were normalized (zero mean and unit variance) within each scan run and averaged across repeated presentations of the same image, serving as functional measures of brain activity.
In our experiments. each image and its corresponding averaged beta map are considered a sample, with an 8:1:1 split into training, validation, and test sets (8000:1000) 1000 samples). The training and validation sets involve non-on overlappiin subjects, while the test set includes repeated samples across subjects. Notably, S1, S2, S5, and S7 completed the full experimental protocol, and their fMRI data are used for analysis in this study.
We use the pre-trained CLIP-ViT model [31] to extract deep visual representations of the stimuli and construct both linear and nonlinear neural encoding models to predict fMRI responses.
In addition, we use the coefficient determination (R²)) to assess model performance on the test set. To determine the statistical significance of the predictions, we adopt the approach described in [2], performing 200 bootstrap resampling iterations on the test data and computing FDR-corrected P-value thresholds for multiple evaluation metrics [2, 50]. As result, all activated voxels reported in this study exhibit statistically significant activation after FDR correction, confirming task- specific neural responses under rigorous statistical criteria.
A.3 Training Strategy
The mean squared error (MSE) [51] is used as the loss function in the encoder, and the training process of the encoder utilizes the Adam optimizer [52] for parameter optimization. To prevent overfitting, early stopping is employed, ceasing training if the validation loss fails to improve over 8 consecutive epochs. Additionally, the model parameters that demonstrate the optimal performance on the validation set are recorded and preserved.
A.4 Efficient Jacobian Computation via Forward-Mode Analytical Differentiation
To address computational costs, traditional backward computation methods (e.g.. using torch.autograd.grut grad in PyTorch with pretrained model and corresponding inputs/outputs) require
15

---

ReLU
Leoky ReLU
GELU
Swish
10-
-10
-10
-5
10
-10
-5
10
-10
1.0
1.0
1.0
0.0
0.0
0.5
0.5
0.6
0.0
0.4
0.0
110
-10
-5
-10
-3
-10
0.5
0.04
0.04
0.00
2nd
0.2
0.9
0.0
-5
Input X
Figure 9: Visualization of activation function response and derivative characteristics. We present the response curves (first row), first-order derivative curves (second row), and second-order derivative curves (third row) of four activation functions-Redd Leaky ReLU, GELU, and Swish-over the input range [0, 10).
approximately several days (Intel Xeon E5-2620 v4 CPU in this study). We optimized this process through explicit forward-mode analytical differentiation. For the nonlinear model employed in our study, the forward structure as follows, where (k) denotes the hidden activations layer k. by represents the bias vectors, WeCk are the weight matrices, and Φ1, Ø2, Φ3 are the activation functions:
(0) = Φ1 (W FCI + b1),
(1) Φ2 (22WWrr
b2) (0)
h (2) = (WF3hh (1) + 03) + (1),
0: WFcch +b
Since commonly used activation functions (e.g., ReLU, GELU) possess analytical first-order deriva- tives, record the sample-specific derivative matrices W during forward propagation and leverage them to compute the full Jacobian:
J3 WFCI W. W.,,WWFFWWW.oo W WFw.W WFos + WFCW Wo.. + FF W & WFC2 W WFr + WFC1 Wo,.WFes.
This approach reduces the computation time to 40-40 seconds, achieving acceleration of over 99.9%.
A.5 Response Curves and Derivative Characteristics of Activation Functions
This section presents the response curves, first-order derivatives, and second-order derivatives of four commonly used activation functions-RRed Leaky ReLU, GELU, and Swish-over the input range E [-10, 10j. By visualizing their local responses and higher-order derivative characteristics, these plots provide intuitive references for analyzing nonlinear structures (Fig.9).
A.6 Details of the Simulation Procedure
To ensure that the evaluation of JNE in the simulation environment both statistically robust and structurally controlled, designed the following experimental procedure:
16

---

Stepl-Model architecture setup. We constructed three-layer feedforward artificial neural network (ANN) (see Fig. 3c), where layers are fully connected. The input and hidden layers are both to 768 dimensions, and the output layer has single unit, simulating voxel-level prediction in neural encoding models. All network weights are randomly initialized prior each experimental run to ensure generalizability and representativeness of the results.
Step2-Activale ctivation function insertion strategy. To introduce nonlinearity, activation functions are inserted between network layers. Four configurations of activation function placement are constructed (see Fig. 3d), ranging from fully linear structure (no activation functions) to configurations where activations are applied between all layers. Identical weights are maintained across all configurations ensure fair comparison.
Step3-Input data generation. The input feature matrix is defined as R00000665 with each row representing high-dimensional feature vector of simulated sample. All elements sampled from a uniform distribution within the fixed interval [-10, 10), ensuring sufficient randommess and variability the input space.
Step4-Jacobian matrix and JNE computation. For each activation function configuration, the Jacobian matrix of the ANN is computed for each sample using automatic differentiation tools implemented in Python. The Jacobian has a shape 1000 x 768 x 1, and is used to calculate the JNE metric. One JNE value is obtained per configuration, resulting in total of four values.
Steps-Repeated experiments for statistical robustness. To enhance statistical reliability, Steps through 4 are repeated 200 times with independently initialized weights. This results in 200 x 4 JNE response matrix, allowing for a systematic characterization of how nonlinearity varies with different activation placements in the ANN.
Step6-Comparison of activation function types. To comprehensively evaluate the nonlinear characteristics of different activation functions, we introduce four commonly used types: ReLU, Leaky ReLU, GELU, and Swish. For each activation function, Steps through 5 are repeated, yielding a complete JNE response matrix with dimensions 4 x 200 x 4 (4 activation types x 200 repetitions placement configurations). which enables detailed analysis of ietiivatoo--ppcifii trends and distinctions under the JNE metric.
A.7 Robustness to Input Sample Size, Distribution and Scale
To address the sensitivity of the JNE metric to the choice of input set and sample size, we conducted controlled simulation experiments. These simulations systematically evaluated the stability and robustness of JNE under variations in sample size, input perturbations, distribution types, and numerical scales. The results indicate that JNE exhibits strong robustness input perturbations and remains stable and reliable even with relatively small sample sizes.
A.7.1 Experiment on Sample Size Variation
The estimation of JNE can be influenced by sample size and input perturbations. To quantify this impact and assess whether moderate sample sizes suffice for stable performance, we conducted a simulation evaluating JNE's ability track known modulation pattern under varying N. The details of simulation experiments are as follows:
Sample size variation: N {10, 20, 30, 50, 100, 200, 300, 500, 1000).
Signal perturbation factors: We generated 1001 evenly spaced points in the interval
x E [0, 2nr], defined a modulating signal sin (x), and used |c] as the input perturbation factor.
Jacobian simulation: For each N, we generated random tensors E IR Nx768x1 A/ 1), then modulated them pointwise using |c] to obtain 1001 Jacobian sequences under varying input amplitudes.
JNE computation: We computed JNE values for each modulated Jacobian sequence using the standard pipeline.
Evaluation metric: We measured the Pearson correlation coefficient r between the resulting
JNE curve and the target trend sin(())] to assess tracking accuracy.
17

---

Generally, from Table I, we find that even with N = 10. the estimated JNE shows strong correlation with the target modulation
0.87). Performance improves with increasing N, achieving near- perfect alignment (r 2 0.99) for N 2 100. These findings indicate that JNE is robust to input perturbations and stabilizes effectively with reasonable sample sizes.
Table 1: Correlation coefficients (r) between JNE and [sin(x)] across sample sizes.
Sample size
10
20
30
50
100
200
300
500
1000
Correlation
0.87
0.94
0.96
0.97
0.99
1.00
1.00
1.00
1.00
A.7.2 Experiment on Input Distribution and Scale
Evaluating sensitivity to input distributions essential for verifying JNE's stability and generalizabil- ity. We systematically varied distribution type and scale while maintaining the network architecture from Network-level validation via ANN-based simulation in Section 3.1. Generally, the simulation experiments are defined as follows:
Distribution type: Inputs were sampled from either uniform distribution U(a, b) or normal distribution N M(u o²).
Numerical scale: Three amplitude ranges were tested small [-100, 100], medium 10000 1000), and large [- 10000, 10000).
Network configurations: Under each input setting, we evaluated JNE across four activa- tion functions (ReLU, Leaky ReLU, GELU, Swish) and four activation location settings (Loc Combination E {1, 3, 4}. see Fig. 3d for details).
Input size and output size: Input fixed at 1000 768, with output as 1000 x 1. Each configuration was sampled 200 times for statistical estimates.
Table 2: Mean JNE values under uniform distributions with varying scales.
Distribution ()
100,
a
1000,
1000
10000,
10000
Loc
Loc
Loc
Loc
Loc
Loc 2
Loc 4
Loc
Loc
Loc 3
0.14
0.16
0.00
0.12
0.16
0.00
0.12
0.14
Leaky ReLU
0.00
0.11
0.13
0.16
0.00
0.12
0.14
0.16
0.00
0.11
0.13
0.16
GELU
0.00
0.12
0.14
0.16
0.00
0.12
0.14
0.16
0.00
0.12
0.14
0.16
0.00
0.12
0.14
0.16
0.00
0.12
0.16
0.00
0.12
0.14
Mean JNE Values under Uniform Input Distributions. Results are highly consistent across scales, indicating insensitivity to input amplitude changes (see Table 2).
Table 3: Mean JNE values under normal distributions with varying variances.
Distribution. N()o2
Jt = 0, o2
100
J == 0,
1000
p = 0, 10000
Loc
Loc
Loc 4
Loc
Loc 2
Loc 4
Loc
Loc
Loc
0.12
0.16
0.12
0.16
0.00
0.14
Leaky ReLU
0.00
0.12
0.13
0.16
0.00
0.11
0.13
0.16
0.00
0.11
0.13
0.16
GELU
0.00
0.12
0.14
0.16
0.00
0.12
0.14
0.16
0.00
0.12
0.14
0.16
Swish
0.12
0.16
0.12
0.14
0.16
0.00
0.14
Mean JNE Values under Normal Input Distributions. Under varying variance configurations of \i o²), JNE remains nearly unchanged, confirming robustness to distribution types (see Table 3).
A.8 Theoretical Justification of JNE as a Proxy for BOLD Response Nonlinearity
The JNE metric quantifies a model's sensitivity to input perturbations, defined as the standard deviation of the Jacobian norm computed across different input samples (i.e., predicted fMRI responses). By the chain rule, the complete onput-output mapping in neural encoding framework can be expressed as:
ml-) 1
(9)
18

---

Here, I denotes the image stimuli, m(.) denotes a pretrained ANN that extracts high-level features from the input I, and f(-) linear or nonlinear neural encoding model that maps the features into the BOLD response space ý. Therefore, applying the chain rule gives:
09 0ū dr 0f(x) dr
(10)
01
dz
01
dz
01
Since the feature extractor m(.) is fixed, the term 0z/01 is constant across all voxels. Therefore, voxel-wise differences in onput-output sensitivity are entirely determined by (x(x)/Ox. From a
physiological perspective, the actual brain response pathway can be approximated
g(-) Neural Activity
x(-)
y
(11)
where g(-) denotes the human brain system, 2(-) denotes the fMRI signal acquisition equipment, and denotes the BOLD responses. If the ANN representation x m( partially approximates the neural state (a widely supported assumption in prior neural encoding work), then 0f(z)/0a can be
viewed as first-order approximation of (z(g( That is:
0f(m(()) 0I
0f(x) aÿ 0z(g(1)) dy dr
(12)
dx aT
aI
Thus, variation in the JNE score reflects the degree of nonlinear sensitivity in the composite system =(g(-). JNE only increases when this system exhibits substantial apput-deeendent nonlinearity, providing theoretical justification for interpreting JNE as approximate proxy for BOLD response nonlinearity.
A.9 Theoretical and Empirical Analysis of JNE Sensitivity
As described in Section 2.4, the JNE metric is designed to quantify nonlinearity in neural encoding models by measuring the dispersion of sample-specific Jacobian deviations from their global mean. While this approach provides principled interpretability tool, it inherently incorporates elements of response gain (magnitude of neural responses) alongside structural nonlinearity and dependencies across input dimensions. This integration is intentional, it allows JNE to reflect physiologically relevant sensitivity differences across voxels or brain regions under standardized input and output spaces (e.g., 7-scored features and beta maps). However, this design choice raises potential concerns about interpretability. such as whether JNE could yield equivalent values for functionally distinct mappings due scaling or dimensional effects. To address these issues systematically, this appendix presents a theoretical derivation of the conditions for JNE equivalence and an empirical evaluation on real fMRI data, demonstrating that such confounds have limited practical impact in neural encoding applications.
A.9.1 Theoretical Derivation
To analyze conditions under which two neurons yield equivalent JNE values, consider a simplified representation of the nonlinear encoding model. For the i-th sample, the predicted responses for two voxels are approximated as linear projections over transformed input features:
U1. = h wT.
f2,1 = h W2:
(13)
where h; [M( A(ID,)) represents the nonlinear transformations of the D-dimensional input x,, and W1, E IR are the learned output projection weights of the two neurons in the nonlinear encoding model.
The Jacobians (gradients) between inputs and outputs are then:
J = Vx,01.5 g. WT, J2 Vx, = g, W2.
(14)
where g: = [0h(z1,) ,)) 00000,,,........ ah(a .)
/)rD.)) [ is the vector of local derivatives.
Assuming zero-mean gradients across N samples (a simplification for analysis), JNE for each voxel is the standard deviation of the absolute Jacobian norms:
JNE; v Var Var (Jr.b-... (JIN))),
JNEg V Var VaJ J2,N
(15)
19

---

Let Σ E RDxD be the covariance matrix of Under the assumption that the entries of g follow Gaussian distribution (approximating the central limit theorem in high dimensions), the absolute norm g w| follows folded normal distribution, yielding:
2
JNE: V 1- 11--- W EWi: JNE2- V
W Ew2-
(16)
Thus, JNE; JNE₂ holds if and only if Ew, = W Ew2 However, in real-world neural encoding models, the fully connected architecture typically induces complex dependencies across input dimensions. As result, Σ is generally non-identity, non-diagonal matrix that includes substantial off-diagonal coupling terms. In such scenarios, even if ww12 o a/w2|2 it does not guarantee a linear scaling relation JNE, = a JNE,, because the contribution of each input direction to JNE is modulated by the covariance structure.
More critically, functional heterogeneity across neurons leads divergent encoding goals, input sensitivities, and nonlinear response characteristics. These structural constraints imply that the output weights w1 and W2- We-lead through training- are unlikely to align direction magnitude, making the assumption W11 00l2 almost never satisfied in practice. Thus, while JNE conflates response gain with nonlinearity in theory, such confounds have limited practical impact due to architectural and functional constraints.
A.9.2 Empirical Evaluation on Real Data
To quantify the practical extent of these theoretical sensitivities, implemented four variants of JNE and applied them to the NSD dataset [30] for S1, computing voxel-wise nonlinearity across significant visual cortex voxels (as defined in Section 3.3). These variants isolate the effects of norm choice and scale normalization, allowing us to assess stability under different formulations:
JNE-L1: Original (centered Jacobians with L1 norm contraction along the feature dimen- sion).
JNE-L2: Centered Jacobians with L2 norm, to evaluate stability across different norm types.
JNE-zscore-L1 Z-score normalization of Jacobian entries across input dimensions before L1 norm, which theoretically mitigates confounding from response gain.
JNE-zscore-L2, Z-score normalization before L2 norm, combining scale invariance with geometric symmetry.
Table 4: Pairwise Pearson correlations between JNE variants NSD fMRI data.
Variant JNE
INE-L
INE-12
JNE-
INF-2score-1.2
1.00
0.95
0.92
INF-L2
0.94
0.92
INE-zscure zscure-LL
0.99
JNF-2score-112
Pairwise Pearson correlation coefficients between the resulting JNE maps shown in Table 4. All correlations exceed 0.9, indicating high consistency despite variations in norm and scaling. For example, the original JNE-LI correlates (e.g., LI vs. L2 in rotationally invariant tunings) do not substantially alter regional patterns. Z-score variants, which mitigate pure gain effects, still align closely (r > 0.90) with unnormalized ones, confirming that scale sensitivity has minimal impact on relative comparisons in this dataset. These findings align with the theoretical analysis, underscoring JNE's robustness for interpretability in neural encoding applications.
1.00 with JNE-L2, suggesting that directional biases
A.10 Theoretical Description JNE-SS
JNE fundamentally characterizes linear nonlinear properties by quantifying the standard deviation of sAJM across different samples at each voxel, which reflects the variability of its mapping relationships across different samples. Further functional analysis of the standard deviation term (sAJM,
pA.)) reveals that
(sAJMi. - Apm)) Es [sAJM Es[sAJM, ...JM
20

---

is essentially second-order moment operator over the sample space S C RY quantifying the covariance structure of the Jacobian matrix within the sample distribution. Its non-zero value directly indicates deviations from the principle of linear superposition.
When the neural system is strictly linear, the Jacobian matrix JM, of the mapping function f: H remains invariant across samples (JM; = JMmea)), causing the sample subspace (sAJM; } to collapse to the origin, leading to J, O. However, nonlinear mechanisms disrupt this invariance, as the local linear approximation (JM,) of a nonlinear system varies with the position of input x; on the manifold, introducing
deviations AsJM,
In this context, the squared deviation term (sAJMi,n pm)) 2 serves microscopic measure of Neural Compliance, capturing the variability of neural encoding at each sample level:
Low compliance (small deviation): Voxel adheres to the linear superposition principle in its response to stimulus changes, exhibiting high stability in neural response.
High compliance (large deviation): Voxel m is governed by higher-order nonlinear mechanisms, with its information transmission exhibiting critical fluctuations.
Thus, (sAJM,,m Apm)) 2 quantifies the sample-dependent dispersion of encoding at individual voxel level, providing insights into the variability of neural responses across samples. Based on this, we define novel metric for quantifying sample-specific nonlinear properties at each voxel, termed JNE with Sample-Specificity (JNE-SS), formally expressed as
JNE-SS (sAJM Aμ)² R
A.11 Dimensionality Reduction, Ordering, and Classification Based on t-SNE
To systematically characterize the trend of voxel-wise nonlinear response properties across stim- ulus samples-ther that is, to identify which samples elicit more linear versus more nonlinear neural responses-we we first applied t-SNE [53] to reduce the dimensionality of CLIP-ViT ViT final-layer repre- sentations for all 73,000 stimulus images. Based on the resulting ID t-SNE embedding, the image samples were then sorted in descending order, and the test set samples were selected accordingly. Since each subject's test set consisted of the same 1,000 images shared across the NSD dataset, this image ordering remained consistent across participants. Fig. 10 presents the sorted images, revealing smooth transition and the clustering of semantically similar samples.
Subsequently, we performed K-means clustering on the 2D t-SNE representation (with the optimal number of clusters selected via silhouette analysis over the range 2-100; see Fig. 11). The 1,000 test images were grouped into 43 clusters, which were then organized according to the 1D t-SNE order. Notably, we observed that these 43 clusters also exhibited smooth transitional structure (Fig. 12).
A.12 Quantitative Comparisons across Regions
Table 5 provides the unnormalized MoM values used in the sample-specifee analysis (Section 3.5), enabling quantitative comparisons across regions without the visual bias introduced by within-ROP normalization in Fig. 8c.
21

---

1:25
26-50
76-100
101-125
126-150
351-175
176-200
201-225
226-250
251-275
276-300
326-350
351-375
401-225
426-400
476-500
501-525
551-575
576-600
626 650
651-675
701-725
726 750
751-775
801-825
826-850
851-875
876-900
901-925
851-975
976-1000
Figure 10: Sample ranking based on t-SNE (ID dimensionality reduction, sorted in ascending order) and visualization of the ranked test set (consistent across subjects).
22

---

Silhouette Coefficlent
0.44
0.42
Dimension
0.40
0.38
-100
10
20
30
40
60
70
80
90
100
120
-30
190
a
Number of clusters in K-means
Dimension
Figure 11: Optimal cluster selection for K-Means based on the silhouette coefficient method (range: 2-100) and 2D dimensionality reduction visualization using t-SNE (with K-Means clustering and cluster number=43).
Labe
30
31
32
34
20
36
41
Figure 12: Image category ranking based on K-Means clustering (number of clusters = 43).
23

---

Table 5: Raw Mean of MoJ (MoM) values across brain regions and image categories.
Categories
EBA
PPA
RSC
OPA
HEA-1
FEA-2
C1 0.016
0.012
0.012
0.014
0.018
0.019
0.020
0.024
0.020
0.015
0.015
C2
0.017
0.016
0.016
0.016
0.023
0.021
0.017
0.020
0.021
0.017
0.18
C3
0.017
0.013
0.013
0.013
0.017
0.021
0.020
0.030
0.017
0.016
0.017
C4
0.015
0,014
0.016
0.019
0.032
0.017
0,019
0.018
0.031
0,040
CS
0.018
0.019
0.018
0.017
0.021
0.019
0.025
0,034
0.020
0.016
0.018
0.014
0.012
0.012
0.013
0.013
0.018
0.020
0.015
0.019
0.014
0.017
C7
0,022
0.016
0,016
0.016
0.015
0.028
0.018
0.020
0.017
0.023
0.030
CR
0.014
0.011
0.011
0.012
0.014
0.022
0.019
0.013
0.016
0.015
0.023
C9
0.023
0.188
0.018
0.018
0.019
0.035
0.022
0.025
0.026
0.032
0.043
C10
0.015
0.012
0.012
0.012
0.014
0.020
0.019
0.016
0.016
0.014
0.019
C11
0,015
0,012
0.012
0.012
0.014
0.020
0.017
0,016
0.016
0.014
0,017
C12
0.023
0.021
0.020
0.020
0.025
0.030
0,023
0,024
0.024
0.021
0.028
C13
0,016
0.013
0.013
0.015
0.018
0.021
0.019
0.016
0.018
0.016
0.021
C14
0,488
0.018
0.018
0.021
0.020
0.067
0.017
0,020
0.016
0,062
0.095
C15
0.021
0.020
0.018
0.021
0.026
0.032
0.018
0.014
0.020
0.023
0.032
C16
0.019
0.011
0.011
0.013
0.013
0.035
0.016
0,015
0.016
0.027
0,038
C17
0.020
0.013
0.012
0.014
0.014
0.032
0.018
0.017
0.016
0.024
0.027
0,017
0.012
0.012
0.013
0.015
0.025
0.019
0.018
0.018
0.018
0.023
C19
0.018
0.012
0.012
0.014
0.015
0.026
0
0,015
0.017
0.019
0.024
C20
0.015
0.012
0.012
0.013
0.014
0.023
0.018
0,018
0.017
0.014
0.018
C21
0,017
0.015
0.014
0,016
0.020
0.021
0.022
0.024
0.024
0.015
0.016
C22
0.017
0.011
0.012
0.013
0.015
0.025
0.023
0.019
0.019
0.015
0.020
C23
0.012
0.012
0.013
0.015
0.022
0.024
0.025
0.022
0.015
0.020
C24
0.021
0.018
0.016
0.016
0.017
0.032
0.018
0.015
0.017
0.022
0.030
0,021
0.015
0.016
0.018
0.021
0.030
0.025
0.023
0.023
0.020
0.025
C26
0.021
0.017
0.016
0.018
0.023
0.034
0.025
0,026
0.026
0.017
0.023
C27
0.015
0.011
0.011
0.012
0.014
0.021
0.020
0.019
0.019
0.015
0.022
C28
0.015
0.012
0.012
0,013
0.014
0.030
0.016
0.015
0.015
0.014
0.021
C29
0.015
0.012
0.012
0.013
0.015
0.018
0.019
0.016
0.016
0.014
0.018
0.011
0.012
0.012
0.013
0.021
0.018
0,016
0.016
0.016
0.022
C31
0.019
0.013
0.013
0.014
0.015
0.031
0,021
0,020
0.020
0.027
0.034
0,023
0.014
0.014
0.016
0.017
0.041
0.022
0.020
0.020
0,029
0.040
C35
0.018
0.011
0.012
0.013
0.014
0.027
0.019
0.018
0.018
0.021
0.027
C34
0.040
0.015
0.016
0.023
0.023
0.097
0.028
0.025
0.025
0.050
0.091
C35
0.021
0.011
0.011
0,013
0.013
0.041
0.015
0.016
0.016
0.022
0.043
C36
0.015
0.011
0.011
0.012
0.014
0.026
0.016
0.016
0.016
0.018
0.026
0.013
0.013
0.016
0.016
0.065
0.016
0.015
0.015
0.048
C38
0.020
0.011
0.011
0.013
0.014
0.035
0.017
0,016
0.016
0.023
0.032
0.023
0.010
0.011
0.013
0.014
0.057
0.018
0.016
0.016
0.027
0.052
C40
0.017
0.012
0.013
0.014
0.016
0.027
0.019
0.018
0.018
0.018
0.021
C41
0.020
0.014
0.014
0.016
0.016
0.037
0,018
0.017
0.018
0.026
0.036
C42
0.020
0.013
0.013
0.015
0.016
0.031
0.023
0.019
0.019
0.022
0.029
C45
0.022
0.015
0.015
0.017
0.018
0.037
0.020
0.020
0.021
0.028
0.037
B Appendix B
B.1
Supplementary Results on Comparison Between Linear and Nonlinear Neural Encoding Models
This section presents results for other subjects (S2, S5, S7). Consistent with the findings reported in the main text, we observed in Fig. 13 that both linear and nonlinear neural encoding models exhibited similar spatial patterns of predictive performance across the whole brain, with peak accuracy localized in visual cortices and several higher-order association areas. These results suggest that, regardless of inter-individual variability, both models reliably capture neural response patterns in the visual cortex.
To further quantify the spatial differences in predictive performance between the two models, we projected the R² scores of the linear and nonlinear models onto a common cortical surface and computed voxel-wise difference maps (Fig. 14). In these maps, white regions indicate comparable performance, red indicates higher accuracy for the nonlinear model, and blue indicates the opposite. The predominance of white regions suggests that both models achieve highly similar predictive accuracy at the voxel level across the brain, with only subtle differences localized in higher-order cortical areas, potentially reflecting nonlinear processing of abstract or semantic information in those regions.
We also performed Pearson correlation analyses on the R2 sequences from the linear and nonlinear models, both across the whole brain and within ten visual-related ROIs. As shown in Fig. 15, all correlations were significantly positive across subjects and ROIs, with the lowest correlation being r = 0.86 in VI of S5 (P < 0.001). These results further confirm the strong consistency in predictive
24

---

Linear
Nonlinearr encoder
S2
encodar
Voninear encodar
S5
Linear
Nontinaar ancoder
S7
Figure 13: Supplementary results comparing the predictive performance of linear and nonlinear neural encoding models. This figure presents the prediction performance of both linear and nonlinear neural encoding models for S2, S5, and S7, The colorbar ranges in all supplementary figures are consistent with those used in the main text.
capacity between the two models and indicate that, when driven by deep visual features extracted from CLIP-ViT, both models characterize IMRI responses in highly similar manner.
From a neuro-computational perspective, these findings suggest that linear neural encoding models are already sufficient capture the neural response structure in the visual cortex when using CLIP-VI VIT features. While nonlinear models offer greater representational flexibility, they do not yield systematic gains in predictive performance under these conditions. Therefore, relying solely differences in R2 scores between linear and nonlinear models is not a reliable strategy for inferring the presence of nonlinear neural response characteristic--paae particclarry when using features that already integrate high-level semantic information.
In conclusion, the cross-subject analysis supports main-text finding: the difference in R between linear and nonlinear models cannot robustly reflect voxel-level nonlincar properties. This provides a theoretical motivation for the development of our JNE metric, which is designed to reveal fine-gu fine-grained comparisons can offer.
structure in neural responses beyond what conventional prediction-
B.2 Supplementary Results Cortical Distribution of Nonlinearity and Hierarchical Effects in Visual Cortex
Fig. 16 presents the cortical distribution BOLD response nonlinearity in additional subjects (S2, S5, S7), further validating the generalizability and consistency of this property across individuals. Consistent with S1, high JNE values were predominantly localized higher-order association cortices, particularly at the temporoparietallocciiiii prefrontal cortex. These regions have long been implicated mul-timodal semantic integration, higher-order cognition, and social reasoning, and are known for their high selectivity to abstract features such as semantics, identity, and context [2. 30, 35, 39, 40, 41]. The observed higher variability in neural responses across different stimulus conditions in these regions likely reflects stronger nonlinear characteristics.
junction (TPOJ), ventral temporal lobe, and medial
Fig. 17 further presents the average JNE values across different hierarchical levels of the visual cortex. Across four subjects, higher-level visual areas consistently exhibited the highest JNE values, while primary visual areas showed significantly lower values, with intermediate areas falling in between. A plot comparison Fig. 18 illustrates this pattern, highlighting that higher-level visual cortices not only have higher mean JNE values but also greater proportion of voxels with high JNE. This pattern
25

---

S1 S2
S5 S7
Figure 14: Mapping of prediction performance for linear and nonlinear neural encoding models on the same cortical atlas. White regions indicate voxels where both models perform similarly; blue regions indicate voxels where the linear neural encoding model outperforms the nonlinear model; red regions indicate the opposite.
Dram
122
PTA
OTA
F1A-2
CSA
3.50
S2
S5
S7
Figure 15: Bar plots of Pearson correlation coefficients of R² values across the whole brain and 10 ROIs for S2, S5, and S7.
supports the presence of a hierarchical gradient in BOLD response nonlinearity within the visual system: the level of information processing increases, neural response variability and nonlinearity also increase.
Fig. 19 further analyzes the distribution of JNE values, presenting the density curves across primary, intermediate, and higher-level visual areas. We observed that JNE values primary areas are concentrated at the lower end, while intermediate areas show broader and slightly right-skewed distribution, and higher-level areas exhibit a clear right-skew with a long-tailed structure. Collectively, Figs. 17, 18, 19, and 20 consistently demonstrate a hierarchical increase in nonlinearity across the visual cortical hierarchy.
Notably, in Fig. 20, the intermediate visual cortex of S2 did not show markedly higher JNE levels compared to the PVC. To clarify this, we further examined JNE values in VI, V2, V3, and V4 of S5, and found that V4 exhibited significantly higher JNE values than V1-V3 (P 0.001). This finding still supports our general conclusion that BOLD response nonlinearity increases with cortical hierarchy and suggests that even in the presence of individual variability, hierarchical patterns of nonlinearity remain evident within specific visual subregions.
In summary, cross- suosect analyses reinforce the conclusion that BOLD response nonlinearity is not uniformly distributed across cortex. Instead, shows pronounced spatial clustering and functional specificity, particularly in higher-order association areas involved in mul-timode integration and semantic abstraction. Furthermore, nonlinearity within the visual pathway follows clear hierarchical organization, with systematic increases from primary to higher-level cortices.
26

---

S2
s5
S7
Figure 16: Spatial distribution of BOLD response nonlinearity in S2, S5, and S7.
S1
S2
S5
S7
0.225
0,200
0.175
0.150
0.100
IVC
HVC
PVG
IVC
HVC
PVC
IVC
100
HVC
Figure 17: Comparison of JNE statistical metrics (mean, standard deviation) across primary visual cortex (PVC), intermediate visual cortex (IVC), and higher-order visual cortex (HVC) for S1, S2, S5, and S7. Primary visual cortex (PVC) VI; Intermediate visual cortex (IVC) V2. V3, V4; Higher-order visual cortex (HVC) EBA, PPA, RSC, OPA, FFA-1, FFA-2.
S2
PvC
S5
S7
30/
HVC
Zer
030
0.3
020
025
0.30
0.00
013
0.2
0.39
0.3
030
0.5
0.00
0.5
Figure 18: The bar plot comparison of JNE primary, intermediate, and higher-order visual cortex for S2, 85, and S7.
27

---

S2
PVC
S5
S7
Densyy
6.
0.5
0.3
3.5
0.5
0.2
0.1
0.2
.NE
Figure 19: The comparison of JNE distribution maps across primary visual cortex, intermediate visual cortex, and higher-order visual cortex for S2, S5, and S7.
S5
S7
PVC
IMC
0
PVC
DVC
HVG
PVC
DVC
INC
PVC
DVC
Figure 20: Significance heatmap of the primary visual cortex, intermediate visual cortex, and higher-order visual cortex for pairs of regions (S2, S5, S7). In the heatmap, value of position (i. j) indicates that the JNE of brain region is significantly greater than that of brain region while value of 0 indicates no significant difference.
S1
S1
S2
S5
S7
Mean (cross subjects)
Brain
3005
Mean of MoM across image calegories
2.00
2005
3:43
$2
V2
2003
V3
s.ns
EBA
PPA
RSC
S7
FFA-1
2003
FFA-2
1005
HNMST
Brain
EBA
PPA
RSC
OPA
FPM-1
FFA-2
b
Label
Figure 21: Statistical (mean) display of MoM values across the whole brain and 10 ROIs, along with the nonlinear category preferences of S1, S2, S5, and S7. Shows the MoM values across the whole brain and 10 ROIs for the four subjects, reflecting the spatial distribution of MoM across different brain regions; Displays the max-min normalized MoM (NMoM) statistical results for the nonlinear response category preferences across the four subjects, with the average values of the four subjects shown characterize statistically significant selective patterns for different categories in nonlinear encoding.
28

---

B.3 Supplementary Results on Sample-Specific Preferences of Nonlinear Neural Responses
We further analyzed the sample-specife preferences of nonlinear neural response characteristics across different subjects, as shown in Fig. 21. Specifically, Fig. 21a illustrates the mean MoM values across all image categories within each visual cortical region for the four subjects. This result further supports the existence of significant spatial inhomogeneity in nonlinear response characteristics across the cortex. An overall trend emerged, suggesting that nonlinearity increases progressively along the visual hierarchy, i.e., PVC < IVC HVC. This hierarchical enhancement of nonlinearity aligns well with classical models of visual processing [36, 37, 38), indicating that neural responses tend to become increasingly nonlinear as the complexity of information processing increases.
Although the overall pattern is largely consistent, also observed individual variability across subjects. For instance, in S7, the mean JNE value of FFA-2. a higher-order visual area, was unexpectedly lower than that of the primary and intermediate visual cortices (Fig.2la). This may reflect subject-specific structural or functional differences in facial stimulus processing. Interestingly. the same subject, we observed notably high JNE level FFA-1, suggesting that within the same functional module, different subregions may show complementary reconfigurable selectivity to specific stimulus categories. Such pattern may arise from the heterogeneous microcircuitry within cortical areas and the diverse semantic feature dimensions being encoded.
Fig. 21b further reveals the category-level selectivity of nonlinear responses. We found that higher- order visual regions exhibited peak nonlinear responses for certain image categories, such as portraits, buildings, and vehiules-stimlli that often contain rich structural and semantic content and thus require more complex perceptual and integrative processes. In contrast, the primary visual cortex showed relatively high nonlinear responses across a broader range of image categories, suggesting that although these early regions are generally considered more linear in processing low-level visual features (e.g., edges, textures, local contrast), the richness of natural stimuli can still elicit complex and possibly nonlinear response patterns due to contextual combinations or stimulus-specife configurations. Overall, these findings complement the results in Fig.
29

---

NeurIPS Paper Checklist
Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: see Section 1.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results be expected to generalize other settings.
is fine include aspirational goals as motivation as long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: see Section 5.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification. asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low images are taken in low lighting. Or speech-to-text system might be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach address problems of privacy and faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question:
For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
30

---

Justification: the theoretical description of the method can be found in Section 2 and Section A.I. Moreover, the analysis of the results is presented in their respective sections.
Guidelines:
The answer NA means that the paper does not include theoretical results. All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: to avoid the disclosure of personal information and the compromise of the double-blind review process, the code and corresponding results will be made open-source after the paper accepted.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer this question will not perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
specific model and empirical evaluation, may
While NeurIPS does not require releasing code, the conference does require submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm.
(b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way access this model for reproducing the results way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reprodueibility may be tricky in some cases, in which authors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path reproducing verifying the results.
5. Open access to data and code
31

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification: the data is from public dataset provided by NSD [30]. The code will be open-sourced after the paper accepted.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips. cc/
cetaanissssooooooooooooooooooooe oeriionolliie Policy) for more details.
While encourage the release of code and data, understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https: / //nips. cc/ wwc/pbllc///aeeecceeeeceecceelee http//uddee//eeeee Ceaees//ooSSubmieec cessinnoolccyy for more details.
The authors should provide instructions data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: Section and Section A.2. Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper 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, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: see Section and their respective sections of results.
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, confi-
dence intervals, statistical significance tests, 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 with given experimental conditions).
32

---

The method for calculating the error bars should be explained (closed form formula, call 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.
OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time execution) needed to reproduce the experiments?
Question:
Answer: [No]
Justification: the code runs on a device with an Intel(R) Xeon(R) CPU E5-2620 v4 @ 2.10GHz and 256 GB of memory. Depending the number of voxels for each subject, the training time per epoch for each subject's 's neural encoding model approximately 80 seconds.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https: http:////wwrrps...pppciiilccceee httpp//pcccctthhccctteeeeeeeeeee
Answer:
[Yes]
Justification: in compliance with ethical regulations.
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 deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive socictal impacts and negative socictal impacts of the work performed?
Answer: [Yes]
Justification: see Section 5.
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 why the paper does not address societal impact.
33

---

Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), faireess 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 direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that 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 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 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 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 monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper deseribe 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)?
Question:
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 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 best faith effort.
12. Licenses for existing assets
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?
Question:
Answer: [Yes]
Justification: we complied with this regulation.
Guidelines:
The answer NA means that the paper does use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include 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 use in package should be provided. For popular datasets, paperswi wwwode.hcoe...c. com datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
34

---

For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer:
[No]
Justification:
the relevant content will be released after the paper is accepted.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moddll 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 used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include anonymized zip file.
14. 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: [No]
Justification: we used publicly available dataset.
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 contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Bthies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. 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 equivalent approval/review based 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 NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, not include any information that would break anonymity (if applicable), such the institution conducting the review.
35

---

16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology. scientific rigorousness, originality of the research, declaration is not required.
Answer: [NA]
Justification: the core method development in this research does not involve LLMs as any important, original, or non-standard components.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, or non-standard components.
Please refer to our LLM policy (https: / /ttp///ersp////// theurips.cc/ooffeeee eeenceee////// /LLMM for what should or should be described.
36


Reviews:
[{"review_id": "tv4Aho6DD2", "reviewer": "Reviewer_m4zX", "summary": "The key contribution of this paper is to provide a new metric for quantifying non-linearity in representations. In particular, they refer to this metric as Jacobian nonlinearity evaluation index (JNE). It measures how the Jacobian of the mapping from inputs to outputs changes over the input over—for a linear model this should be 0, while for a non-linear model we should see higher values for higher sensitivity to the space of inputs being considered. They apply this metric to study the amount of non-linearity in visual representations in different parts of the human visual cortex, by mapping embeddings of images to their voxel responses. \n\nAligning with previous results, they find that high-order visual orders have the most non-linearity, and that nonlinearity is hierarchically organized in the visual cortex. Separately, they also compute sample selectivity of representations in different regions and again consistent with literature they find that regions like FFA show pronounced non-linearity in response to faces (which they are known to be encoding).", "questions": "Following up on the weaknesses listed above:\n\n1. I would like to know computational costs for these experiments to get a sense of how feasible they are to be used across domains.\n2. While for the real experiments, it makes sense to stick to the data collected, it would be nice to know how sensitive the metric is to input samples in the simulations in order to strongly trust the real data results. \n3. I am also curious about the choice of L1 norm here to measure deviation from the mean, and would like more clarification on that.", "limitations": "yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n\n1. I absolutely love the figures and how the paper is written, everything is very clean and well-organized.\n2. The metric is interesting, it is a novel way of quantifying non-linear between ANNs and brain encoding that goes beyond prediction performance.\n3. Such mapping is widely studied in the human neuroscience community across modalities: audio, vision, etc and so this metric could be influential in extracting similar information about encodings in other parts of the brain.\n\nWeaknesses:\n1. This is minor but the limitation section should really go into the main text, and I'd urge the authors to do so if the paper gets accepted (in which case they will have an extra page).\n2. Computing the Jacobian for a large model for various values of the input is a computationally expensive operation (as the authors detail in the limitations). The computation really depends on the input-output dimensionality, and size of sample space, which may restrict applicability. \n3. A crucial part of computing this metric is sampling the right set and number of inputs over which the Jacobian is computed. This hasn't really been explored in the paper as of yet.", "was_revised": true, "final_justification": "I am keeping my original score of 5. The paper was strong to start with, and the authors have thoughtfully answered my questions during the rebuttal phase.", "rebuttal": "We would like to thank the reviewer for the constructive comments and suggestions, as well as the appreciation of our study. Our point-by-point responses are detailed as follows.\n\n**W1:** *This is minor but the limitation section should really go into the main text...*\n\n**Response:** We will follow the reviewer's suggestion and move the \"Limitations\" section to the main text.\n\n**W2:** *Computing the Jacobian for a large model is computationally expensive...*\n\n**Response:**  Estimating the Jacobian of the models used in the current study via sample-wise backward computation, as traditionally done, may require 50 hours on a machine equipped with an Intel Xeon E5-2620 v4 CPU . To reduce this burden, our implementation employs an explicit forward-mode inference based on the structural formulation of the model. Specifically, we derive the closed-form Jacobian by analytically differentiating the linear transformations and activation functions at each network layer. This precedure relies solely on intermediate values and derivatives that are saved during the forward pass. As a result, this strategy reduces the computational time to 20–40 seconds, achieving an acceleration factor of over 99.9%.  We will provide implementation details and open-source the code. Below are the detailed formulation of the implimentation.\n\nFor the nonlinear model used in our study, the forward structure is as follows:\n\n```math\nh^{(0)} = \\phi_1(W_1 x_i + b_1) \\\\\nh^{(1)} = \\phi_2(W_2 h^{(0)} + b_2) + h^{(0)} \\\\\nh^{(2)} = \\phi_3(W_3 h^{(1)} + b_3) + h^{(1)} \\\\\n\\hat{y}_i = W_4 h^{(2)} + b_4\n```\n\nSince commonly used activation functions (e.g., ReLU, GELU) have analytical first-order derivatives, we can record the sample-specific derivative matrices $W_{\\phi_j, i}$ during forward propagation, and use them to calculate  the full Jacobian:\n\n```math\nJ_i = W_1 W_{\\phi_1,i} W_2 W_{\\phi_2,i} W_3 W_{\\phi_3,i} W_4 \\\\\n+ W_1 W_{\\phi_1,i} W_3 W_{\\phi_3,i} W_4 \\\\\n+ W_1 W_{\\phi_1,i} W_2 W_{\\phi_2,i} W_4 \\\\\n+ W_1 W_{\\phi_1,i} W_4\n```\n\n**W3:** *A crucial part of computing this metric is sampling the right set and number of inputs...*\n\n**Response:** The choice of input set and the number of samples indeed have a significant impact on the esitmation of JNE. To alleviate this problem, in the original submission, we applied multiple random sampling with averaging to reduce randomness and enhance statistical robustness in simulations. In real-data analysis, the training, validation, and test sets were randomly partitioned.\n\nTo directly address the concerns such as \"*Is the number of samples sufficient?*\" and \"*Is JNE sensitive to input perturbations?*\", we inaddition designed a controlled simulation to quantitatively assess the impact of sample size. Our experimental results show that JNE is robust to input perturbations. When supplied with a reasonable number of samples (N>50), JNE can stably works well.\n\nThe details of simulation experiments are as follows:\n\n1. **Sample size variation**: We tested a range of sample sizes $N \\in [10, 20, 30, 50, 100, 200, 300, 500, 1000]$;\n2. **Signal perturbation factors**: We generated 1001 evenly spaced points in the interval $[0, 2\\pi]$, defined a modulating signal $c = \\sin(x)$, and used $|c|$ as an input perturbation factor;\n3. **Jacobian simulation**: For each sample size $N$, we generated random tensors $J \\in \\mathbb{R}^{N \\times 768 \\times 1} \\sim \\mathcal{N}(0, 1)$, then modulated them pointwise using $c$ to obtain 1001 Jacobian sequences under varying input amplitudes;\n4. **JNE computation**: We computed the JNE values corresponding to each modulated Jacobian sequence;\n5. **Evaluation metric**: We measured the Pearson correlation coefficient between the resulting JNE curve and the target trend $|\\sin(x)|$, to assess JNE's ability to track the input modulation.\n\nThe results are summarized as follows:\n\n| Sample size$N$ | 10   | 20   | 30   | 50   | 100  | 200  | 300  | 500  | 1000 |\n| ---------------- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- |\n| Correlation$r$ | 0.87 | 0.94 | 0.96 | 0.97 | 0.99 | 1.00 | 1.00 | 1.00 | 1.00 |\n\nEven with as few as $N = 20$, the estimated JNE shows a strong and significant correlation with the target modulation pattern ($r = 0.87$). When $N \\geq 100$, JNE almost perfectly fits the target curve.\n\n**Q1:** *I would like to know computational costs...*\n\n**Response:** Please see our response to W1.\n\n**Q2:** *It would be nice to know how sensitive the metric is to input samples...*\n\n**Response:** We aggree with the reviewer that evaluating the sensitivity of the JNE metric to input distributions is a critical step toward verifying its stability and generalizability. To this end, we extended the original simulation setup (Section 3.1) by systematically varying the input configuration along two dimensions—distribution type and numerical scale—to assess the robustness of JNE. Our experimental results demonstrate the high robustness of JNE.\n\nSimulation experiments are detailed as follows.\n\n1. **Distribution type**: Inputs were sampled from either a uniform distribution $\\mathcal{U}(a, b)$ or a normal distribution $\\mathcal{N}(\\mu, \\sigma^2)$;\n2. **Numerical scale**: Three amplitude ranges were tested—small $[-100, 100]$, medium $[-1000, 1000]$, and large $[-10000, 10000]$.\n\nUnder each configuration, we evaluated JNE across four commonly used activation functions (ReLU, Leaky ReLU, GELU, and Swish), and four activation location settings (Loc Combination ∈ {1, 2, 3, 4}).\n\n**Mean JNE Values under *Uniform* Input Distributions**\n\n| Distribution Range   | $[-100, 100]$`<br>`Loc 1 | Loc 2 | Loc 3 | Loc 4 | $[-1000, 1000]$`<br>`Loc 1 | Loc 2 | Loc 3 | Loc 4 | $[-10000, 10000]$`<br>`Loc 1 | Loc 2 | Loc 3 | Loc 4 |\n| -------------------- | ---------------------------- | ----- | ----- | ----- | ------------------------------ | ----- | ----- | ----- | -------------------------------- | ----- | ----- | ----- |\n| **ReLU**       | 0.00                         | 0.12  | 0.14  | 0.16  | 0.00                           | 0.12  | 0.13  | 0.16  | 0.00                             | 0.12  | 0.14  | 0.16  |\n| **Leaky ReLU** | 0.00                         | 0.11  | 0.13  | 0.16  | 0.00                           | 0.12  | 0.14  | 0.16  | 0.00                             | 0.11  | 0.13  | 0.16  |\n| **GELU**       | 0.00                         | 0.12  | 0.14  | 0.16  | 0.00                           | 0.12  | 0.14  | 0.16  | 0.00                             | 0.12  | 0.14  | 0.16  |\n| **Swish**      | 0.00                         | 0.12  | 0.14  | 0.16  | 0.00                           | 0.12  | 0.14  | 0.16  | 0.00                             | 0.12  | 0.14  | 0.16  |\n\nThese results are highly consistent across all tested ranges, indicating that JNE is insensitive to changes in input scale.\n\n**Mean JNE Values under *Normal* Input Distributions**\n\n| Distribution$\\mathcal{N}(0, \\sigma^2)$ | $\\sigma^2 = 100$`<br>`Loc 1 | Loc 2 | Loc 3 | Loc 4 | $\\sigma^2 = 1000$`<br>`Loc 1 | Loc 2 | Loc 3 | Loc 4 | $\\sigma^2 = 10000$`<br>`Loc 1 | Loc 2 | Loc 3 | Loc 4 |\n| ---------------------------------------- | ------------------------------- | ----- | ----- | ----- | -------------------------------- | ----- | ----- | ----- | --------------------------------- | ----- | ----- | ----- |\n| **ReLU**                           | 0.00                            | 0.12  | 0.14  | 0.16  | 0.00                             | 0.12  | 0.14  | 0.16  | 0.00                              | 0.12  | 0.14  | 0.16  |\n| **Leaky ReLU**                     | 0.00                            | 0.12  | 0.13  | 0.16  | 0.00                             | 0.11  | 0.13  | 0.16  | 0.00                              | 0.11  | 0.13  | 0.16  |\n| **GELU**                           | 0.00                            | 0.12  | 0.14  | 0.16  | 0.00                             | 0.12  | 0.14  | 0.16  | 0.00                              | 0.12  | 0.14  | 0.16  |\n| **Swish**                          | 0.00                            | 0.12  | 0.14  | 0.16  | 0.00                             | 0.12  | 0.14  | 0.16  | 0.00                              | 0.12  | 0.14  | 0.16  |\n\nUnder different variance configurations of $\\mathcal{N}(0, \\sigma^2)$, JNE remains almost unchanged, demonstrating the robustness of JNE to various distributions of inputs. These results suggest that JNE can reliably capture structural nonlinearity under a wide range of input settings.\n\n***Q3:** The choice of L1 norm*\n\n**Response:** We adopted the L1 norm primarily due to its robustness and stability in measuring deviations from the mean in high-dimensional spaces. Compared to L2 norm, L1 norm computes absolute deviations, which mitigates the influence of extreme values. This makes it more effective in capturing the total magnitude of nonlinear perturbations, without being disproportionately affected by outliers."}, {"review_id": "iPV9p3Dn5Y", "reviewer": "Reviewer_gFGX", "summary": "This paper introduces a measure (called JNE) to functionally compare output dimensions that are nonlinear functions of the input, focusing on nonlinear encoding models of biological sensory neurons. Using simple simulations, they show that JNE increases as model nonlinearities increase, and they also show JNE computed on models of visual cortical areas reflects hierarchical (and increasingly nonlinear) computations along visual areas.", "questions": "### Below are the questions for which I am happy to increase the rating given that they are convincingly addressed:\n\n1\\. How is the scale-dependence of the method justified? The current formulation conflates response gain with degree of nonlinearity, which is problematic: due to its scale dependence, the metric does not accurately reflect the intrinsic dimensionality or structure of curvature. This can yield identical scores for functionally distinct mappings through arbitrary rescaling—conflating response gain with actual nonlinearity. To illustrate this, consider two functions with inputs $x_1, x_2 \\sim \\mathcal{N}(0, 1)$:\n\n* $f_1(x_1, x_2) = a x_1^2 \\Rightarrow \\nabla f_1 = [2a x_1,\\ 0]$\n* $f_2(x_1, x_2) = x_1^2 + x_2^2 \\Rightarrow \\nabla f_2 = [2x_1,\\ 2x_2]$\n\nThen the JNE (defined as the standard deviation of the L1 norm of the Jacobian across inputs) becomes:\n\n* $\\text{JNE}(f_1) = \\mathrm{std}(|2a x_1|) = 2a \\cdot \\mathrm{std}(|x_1|)$\n* $\\text{JNE}(f_2) = \\mathrm{std}(|2x_1| + |2x_2|) = 2 \\cdot \\mathrm{std}(|x_1| + |x_2|)$\n\nMatching the JNEs yields: $ 2a \\cdot \\mathrm{std}(|x_1|) = 2 \\cdot \\mathrm{std}(|x_1| + |x_2|) \\Rightarrow a = \\frac{\\mathrm{std}(|x_1| + |x_2|)}{\\mathrm{std}(|x_1|)}$\n\nThis shows that **a 1D nonlinear function** can be made to appear **equally nonlinear** as a 2D one by simply amplifying its output. I am not convinced that response gain should be entangled with nonlinearity in a single metric, as they capture distinct computational aspects. This could be mitigated by standardizing Jacobian entries across inputs (e.g., subtracting mean and dividing by std), rather than just centering them. Why did the authors choose not to do this?\n\n2\\. Why use L1 norm and not L2? The proposed Jacobian-based metric aims to quantify a neuron's nonlinearity by aggregating input-specific gradients. However, for neurons with symmetric and rotationally invariant tuning profiles (i.e., invariance under rotations of the input space), the use of the L1 norm introduces a directional bias that breaks isoresponse symmetry. For example, consider a neuron with a 2D Gaussian tuning function: $f(x_1, x_2) = \\exp\\left(-\\frac{x_1^2 + x_2^2}{2\\sigma^2}\\right)$. The gradient is $\\nabla f = -\\frac{f(x)}{\\sigma^2} [x_1, x_2]$, and the **L2 norm** becomes: $\\lVert \\nabla f \\lVert_2 = \\frac{f(x)}{\\sigma^2} \\sqrt{x_1^2 + x_2^2}$. But since $x_1^2 + x_2^2 = -2\\sigma^2 \\log f(x)$, this yields: $\\lVert \\nabla f \\lVert_2 = \\sqrt{\\frac{-2 f(x)^2 \\log f(x)}{\\sigma^2}}$, which depends only on the function value $f(x)$ — and thus is constant along isoresponse contours. In contrast, the **L1 norm** is: $\\lVert\\nabla f\\lVert_1 = \\frac{f(x)}{\\sigma^2} (|x_1| + |x_2|)$, which varies across directions even when $f(x)$ is constant (i.e., on the same contour). In other words, equally nonlinear points are scored differently simply due to orientation. Using L2 could address this specific issue, and I wonder why the authors chose L1 norm.\n\n3\\. Why is the method \"scale-normalized\" as mentioned in the paper? which operation/step is related to this scale normalization? and why does that mean “the magnitude of standard deviation directly encodes biological sensitivity differences in neural responses”?\n\n4\\. I like that this approach is architecture-agnostic and can be applied on any model as long as the Jacobian matrix can be computed, but its dependence on the input introduces other challenges/questions: How and to what degree does JNE depend on the input distribution? How many inputs are needed? how to select the inputs? how does the number of inputs (and their selection) change with different number of output dimensions, and their degree of nonlinearity?\n\n5\\. While the method recovers the expected hierarchical trend in visual areas—consistent with prior findings—it is not immediately clear to me what new insights it enables beyond this. What does the metric reveal beyond \"this area/neuron is more/less nonlinear than another\"? The sample-specific analysis seems to hint at deeper insights, but ultimately still confirms known properties (e.g., early areas encode simple features). Could the authors clarify the kinds of new questions or discoveries this metric enables?\n\n6\\. Regarding the sample-specific analysis, I am not sure how much one can rely on or interpret these results. For instance, you state that “stimuli involving people or faces induce pronounced nonlinear fluctuations in FFA and EBA.” Looking at cluster 14 in Figure 12, this seems to hold—there are faces in the images, and FFA shows high MoM in the heatmap (Figure 8c). However, cluster 9 also contains many face stimuli, yet FFA does not exhibit similarly high MoM. Could the authors clarify this discrepancy or comment on the robustness of these observations?\n\n---\n\n### Additional questions\n1\\. The interpretation of JNE in terms of “stimulus sensitivity” was unclear to me. The authors write that \"stimulus-sensitive regions (e.g., early visual cortex) exhibit greater parameter fluctuations, which mathematically aligns with the sample-dependent characteristics of nonlinear mechanisms\" seemingly equating higher Jacobian variability (across input dimensions) and increased JNE. Could the authors clarify what does “parameter fluctuation” mean, in this context? and how does that “mathematically align with the sample-dependent characteristics of nonlinear mechanisms”?\n\n2\\. In figure 3d,e combination 2 and 3 have different JNE values, and you do mention it as an interesting observation but I could not find an explanation of why this happens. Intuitively in both cases you have linear → nonlinear → linear, and it is not clear to me why they should have different degrees of nonlinearity.\n\n3\\. In figure 7c, the caption says: where entry $(i, j)=1$ that means cortex i exhibits significantly lower JNE than cortex j. This means for instance HVC ($i=3$) exhibits significantly lower JNE than PVC ($j=1$). Shouldn’t it be the other way around? i.e. $(i, j)=1$ indicates statistically significantly “higher” as opposed to “lower”?\n\n4\\. What do the colors in figure 7d correspond to? Unfortunately, I did not get that from the caption. Is this a voxel-specific coloring or the voxels are based on the area label and the color is simply coming form the statistical test of difference between the areas? If it is the latter I am not sure what is the contribution of this figure exactly.\n\n5\\. I found the appendix defining JNE-SS a bit of an overkill. Isn’t JNE simply the square-root of JNE-SS averaged across stimuli? Changing the order to define JNE-SS first and then getting to JNE seems more intuitive to me, simply because JNE is an aggregate measure of JNE-SS, and maybe it could help simplify the definition/explanation.", "limitations": "Some limitations are acknowledged by the authors and could serve as directions for future work. However, key limitations/questions—such as the method’s sensitivity to scale, and dependence on the input distribution were not discussed. These aspects are important for interpreting the utility and generalizability of the proposed metric and it is important to communicate them.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "### Strength\n\n- Technically sound derivation; tractable and well-motivated; reasonable evaluation on simulations and hierarchy\n- Generally well-written, logical flow, readable\n- Addresses an important topic—quantifying single-unit nonlinearity with a scalar, model-agnostic metric\n- The proposed approach is novel in that it aims to provide a direct measure of nonlinearity that can be applied across different models, filling a gap not directly addressed by model comparison or geometry-based approaches.\n\n### Weaknesses\n\n- There are some conceptual issues (scale sensitivity, isoresponse inconsistency); some design choices are not clear (e.g. why L1 norm?)\n- Some unclear phrasing (\"stimulus sensitivity\", “parameter fluctuations”); some definitions and the relation of different terms could be improved (e.g., JNE vs JNE-SS).\n- Practical utility depends on resolving open issues (input distribution characteristics, scale sensitivity). The impact/significance of the method depends on how these are addressed.\n- The interpretability of the sample-specific MoM analysis is unclear; some stimulus clusters with similar content (e.g., faces) elicit inconsistent responses in expected regions (e.g., FFA), raising questions about the robustness or specificity of the metric.\n- Overlap with existing conceptual directions (intrinsic dimensionality, representational geometry); novelty could be better articulated relative to prior work.", "was_revised": true, "final_justification": "The authors have provided a clear and well-argued response to my concerns, including a general derivation clarifying under what conditions JNE values can be equivalent across neurons, and why such equivalences are unlikely in practice. They also provided additional empirical validation and improved clarity around the metric’s design trade-offs. Overall, the authors' handling of these issues and the practical stability of the metric lead me to consider the contribution valuable and above the acceptance bar. I am therefore increasing my rating.", "rebuttal": "We would like to thank the reviewer for the constructive comments and suggestions, as well as the appreciation of our study. Our point-by-point responses are detailed as follows.\n\n### **Response to Questions**\n\n**Q1:** *How is the scale-dependence of the method justified?*\n\n**Response:** We appreciate the reviewer's insightful example, which demonstrates that JNE is sensitive to output scale. However, it is important to clarify that JNE was not designed to quantify curvature in a purely mathematical sense, but rather to serve the interpretability of nonlinear neural encoding models. In this context, we intentionally preserve the magnitude information in the Jacobian standard deviation, based on the assumption that the *amplitude* of model-predicted responses may itself carry neurophysiological significance. For example, the overall responsiveness of a brain region or neuron to a stimulus may reflect biologically meaningful differences in sensitivity.\n\nGrounded in this motivation, we did not remove scale effects from the JNE formulation. Instead, we allow the absolute size of the Jacobian variability to reflect such sensitivity differences. Since the \"structural space\" of the input-output mapping is fixed during training, the response amplitude becomes an intrinsic part of the learned mapping. In other words, the model learns not only the shape of the curvature or representational manifold, but also its corresponding magnitude. Therefore, we argue that retaining response gain in the metric facilitates more direct links between the degree of nonlinearity and biological properties of the system, thereby enhancing physiological interpretability. This trade-off was a deliberate design choice in JNE: while we may sacrifice some theoretical purity in terms of scale-invariant curvature measurement, we gain a more direct characterization of real neural response behaviors. To further illustrate this, we note the following structural comparison between two functions:\n\n$$\nf_d(x) = \\sum_{i=1}^d x_i^2, \\quad f_{1D}(x) = a x_1^2,\n$$\n\nUnder $x \\sim \\mathcal{N}(0, I_d)$, we have:\n\n$$\n\\text{JNE}(f_d) = 2 \\sqrt{d} \\cdot \\sigma_{|x|}, \\quad \\text{JNE}(f_{1D}) = 2a \\cdot \\sigma_{|x|} \\Rightarrow a = \\sqrt{d}.\n$$\n\nThis result shows that to match the JNE value of a $d$-dimensional quadratic function, one must scale a single-dimension function by $a \\sim \\sqrt{d}$, e.g., $a \\approx 28$ when $d = 768$. Such extreme amplification in a single input dimension is practically implausible in real neural encoding models, due to constraints from regularization, input-output normalization, parameter initialization, and activation saturation. These mechanisms jointly make it nearly impossible for a model to concentrate nonlinearity in a single direction at such scale.\n\nIn the revised manuscript, we will explicitly state that the scale sensitivity of JNE is intentional, as it supports the interpretability of nonlinear responses in biological systems rather than undermining theoretical soundness.\n\n**Q2:** *Why use L1 norm?*\n\n**Response:** Our choice of the L1 norm was primarily motivated by its intuitive interpretability: it captures the overall response magnitude of a voxel to input perturbations via the sum of absolute gradient values across dimensions. This makes it more amenable to potential associations with neural sensitivity. In contrast, L2 norm, while geometrically symmetric, may obscure directional differences in local responses. It is important to emphasize that the goal of JNE is not to characterize the shape of iso-response profiles, but rather to compare the amplitude of nonlinear fluctuations across voxels (brain regions) under the same input perturbation distribution. Accordingly, we prioritize perturbation sensitivity, not geometric invariance. In addition, we tested the L2 norm and found that the overall spatial distribution pattern of JNE on the cortical surface was identical with that from L1, despite neumerical difference. We will clarify this design choice in the revised manuscript.\n\n**Q3:** *Why is the method \"scale-normalized\" ?...*\n\n**Response:** In both encoder training and JNE computation, we applied z-score normalization to the input features (i.e., CLIP representations) so that each input perturbation dimension had a standardized scale (0 mean, unit deviation). The target outputs of the encoder (fMRI beta values) were also z-scored prior to training, and both the model architecture and training protocol were kept consistent across voxels. Under this setting, the input space, output space, and mapping space (i.e., nonlinear neural encoding models) are all standardized. The normalized test inputs are likewise constrained within this normalized feature space, ensuring that both the training and evaluation phases operate under the same scale assumptions and statistical priors. As a result, the variability of the Jacobian can be directly interpreted as the relative sensitivity of each voxel to input perturbations, supporting comparisons across voxels and brain regions.\n\n**Q4:** *I like that this approach is architecture-agnostic ...*\n\n**Response:** We thank the reviewer for the appreciation of JNE. As a statistical measure, JNE is indeed influenced by the input distribution and sampling range. Sparse or skewed sampling may lead to an underestimation of the overall nonlinearity.\nTo systematically assess its robustness, we conducted additional simulation experiments by varying both the sample size,  input scale and input distribution. Results show that even with as few as 30 samples, JNE can reliably works well under various input scales and distributions. Moreover, expanding the input scale from $[-100, 100]$ to $[-10^4, 10^4]$ produced highly stable JNE values, suggesting that the metric is primarily driven by model structure, rather than input scale alone.\nIn the case of real fMRI data, we recognize that the stimulus set is often constrained by the experimental task, leading to non-ideal input distributions. Similar to $R^2$, we emphasize that JNE is a conditional metric—its value reflects nonlinear response characteristics under a specific model–data configuration, and are best interpreted in a relative sense across brain regions or models, rather than as an absolute measure of nonlinearity.\n\nIn the revised manuscript, we will incorporate these results and discuss the limitation of JNE in interpreting real fMRI data.\n\n**Q5:** *While the method recovers the expected hierarchical ...*\n\n**Response:** The primary contribution of this work lies in the proposal and validation of the JNE metric for interpreting nonlinear encoding models.  While we did observe some interesting patterns (e.g., category-specific nonlinearity preferences), these were intended solely as preliminary demonstrations of the feasiablity of the proposed method. The validation of those findings requires additional experiments involving multiple vision-language models, multiple fMRI datasets and more sophisticated statiscal comparisons. In the revised manuscript, we will revise the wording to better emphasize the methodological contribution of our work.\n\n**Q6:** *Regarding the sample-specific analysis...*\n\n**Response:** In Cluster 9, the MoM value for FFA is higher than both the whole-brain average and early visual areas such as V1–V4 (see table below). However, since the figure applies max–min normalization within each region, this difference may be visually attenuated. In the revised manuscript, we will explicitly clarify this normalization procedure to avoid misinterpretation during cross-region comparisons.\n\n| Region | Brain | V1    | V2    | V3    | V4    | EBA   | PPA   | RSC   | OPA   | FFA-1 | FFA-2 |\n| ------ | ----- | ----- | ----- | ----- | ----- | ----- | ----- | ----- | ----- | ----- | ----- |\n| MoM    | 0.023 | 0.018 | 0.018 | 0.018 | 0.018 | 0.035 | 0.022 | 0.025 | 0.026 | 0.032 | 0.043 |\n\n### **Response to Additional Questions**\n\n**AQ1:** *The interpretation of JNE in terms of...*\n\n**Response:** In the manuscript, the term \"stimulus sensitivity\" refers to the degree to which a brain region's Jacobian mapping depends on the input sample, i.e., the extent of variation across different visual stimuli. Similarly, \"parameter fluctuation\" denotes the variation of the Jacobian vector across input samples. These two terms essentially describe the same phenomenon: the variability of the model's Jacobian mapping in response to input perturbations. We will make these terms consistent in the revised manuscript.\n\n**AQ2:** *In figure 3d,e combination 2 and 3...*\n\n**Response:** Unfortunately, despite our best efforts to derive a formal explanation, the underlying mechanism remains unclear at this stage. We hypothesize that when nonlinearity is positioned in the later layers of a network, the longer gradient propagation path may increase the proportion of extreme Jacobian values, thereby elevating the JNE score. This pattern remains consistent across different activation functions and network architectures, suggesting that the position of nonlinearity may systematically affect its distribution. We plan to further investigate the underlying mechanism in future work and will design controlled experiments to validate this effect.\n\n**AQ3:** *In figure 7c, the caption says...*\n\n**Response:** Thank you for pointing this out and we will correct it accordingly.\n\n**AQ4:** *What do the colors in figure 7d correspond to?...*\n\n**Response:** You are absolutely correct: the colors in Figure 7d do not reflect voxel-level specificity, but are instead based on the JNE significance results aggregated at the region level. The purpose of this figure is to validate the effectiveness of the JNE metric, rather than to localize voxel-wise effects.\n\n**AQ5:** *I found the appendix defining JNE-SS...*\n\n**Response:** This is a great suggestion that help us to improve the clarity and readability of the manuscript. We will follow this suggestion to revise the manuscript."}, {"review_id": "D3UW8BpB4m", "reviewer": "Reviewer_VQtJ", "summary": "The authors propose a non-linearity metric (Jacobian Non-linearity Evaluation – JNE) to characterize the non-linear response properties of the brain. The JNE for mapping Clip-VIT output activations to brain responses (from the Natural Scenes Dataset) showed increasing JNE scores from primary visual cortex, to intermediate (V2-V4), to higher level visual cortex (e.g., EBA, FFA, PPA, etc.). Moreover, JNE scores were higher for preferred categories in higher visual cortex (e.g., images of outdoor scenes evoked stronger JNE in PPA and RSC; images of people and faces evoked stronger JNE in FFA and EBA). The authors conclude that JNE characterizes the degree of non-linearity in nonlinear neural encoding models, and describes the degree of nonlinearity across different brain regions.", "questions": "1)\tThe JNE score measures the degree of non-linearity in the “mapping” from model-space to voxel-space. Is it the case that high JNE scores can be observed if-an-only-if the target space exhibits high non-linearity? If this is not a provable constraint on JNE scores, we should be cautious about interpreting the mapping-JNE scores as a measure of the encoding-target-space-nonlinearity. \n2)\tOn this same topic, I’m not sure what we can/should conclude about the brain regions from these JNE scores. To illustrate my concern, let’s pretend that some Clip-ViT output units respond more to “profile-faces” and others to “frontal views of faces”, and we are predicting FFA responses (face-selective cortex). Given that feature-response preference, I assume the JNE score would be higher if you presented a diverse image set (mixture of profile vs. frontal views), and lower if you presented a homogenous set (frontal views only). The brain region itself has “the same non-linearity” no matter what image set you present, but the JNE scores would vary. \n3)\tOne more question on interpretation. The fact that the linear and non-linear encoders show equivalent performance mapping between CLIP-output and brain responses, it really seems like JNE scores are capturing the non-linearity of the mapping (which itself seems like it would depend on the variability of the dataset along dimensions captured by the encoder; see previous question), rather than non-linearity of the brain-response per se. \n4)\tI found myself wondering whether/how the JNE scores would vary as a function of model-layer. Typically earlier model layers show a stronger alignment with primary and intermediate visual cortex, whereas later model layers show stronger alignment with higher visual cortex. To what extent do JNE scores depend on the model-layer that feeds into the encoding model? Is the low JNE score for primary visual cortex due to the mismatch between model-space (CLIP-output activations)?\n5) In theory, with a powerful enough encoder, it should be possible to map even an early model layer to higher-visual cortex, presumably requiring greater non-linearity in the mapping. Does this mean higher JNE scores? If so, given that the target brain regions aren’t changing, that variability would imply that JNE is a measure of the non-linearity in the mapping, not the target brain region.\n6)\tIt would be helpful if the authors could situate this work more fully in the context of what is known about the biological visual processing hierarchy, it’s known increase in complexity (particular selectivity profiles + invariances), and non-linearity in visual processing (providing some clarity on how this is measured and quantified in biological systems). It would also be helpful to relate the current approach/findings to what has been demonstrated in prior work, particularly work showing a hierarchical correspondence between model activations and brain responses (early model layers predicting early brain regions better, later model layers predicting later model layers better), which has also been proposed to reflect comparable non-linear stages in artificial and biological visual systems. What new insights do we gain from the current approach?\n7)\tIf the goal is to characterize properties of the brain, independent of the encoding model, it seems important to present results for encoding models mapping multiple model layers (from early to late to output layers) and multiple models (various architectures, visual diets, and objectives) to ensure that the conclusions about the brain are general/consistent across layers and models. My guess is that there will be significant variability based on layer (and possibly model depending on the suite drawn from), and that more nuanced conclusions regarding non-linearities in the mapping from model-space to brain-space are required.\n8)\tIn the Discussion the authors propose that the present results support a bottom-up processing mechanism in the brain’s visual system. I’m not sure I understand exactly what the authors mean here, and what the alternative is. The biological visual system is known to have ubiquitous feed-forward, feed-back, and skip connections, and fMRI responses are too sluggish to isolate purely feed-forward processing, so the signal being modeled here is necessarily a result of both feed-forward and feed-back processing. If the authors mean that the input signal (image) is the primary driver of responses along the ventral processing stream, I would agree with that, but I’m not sure the present results speak to that much. \n\nI want to emphasize that I think this is very solid work and that the results are intriguing enough for me to think the authors are onto something here. The hierarchical progression of JNE scores with brain-depth, and the category-dependent JNE scores (faces in FFA; scenes in PPA) are very interesting results which suggest the measure has potential to provide some new insight. I’m just not sure we can safely interpret the JNE values as a measure of non-linear processing in the target brain region. If not, if they are a measure of the non-linearity in mapping model-layer-A to brain-region-B, then I’m not sure exactly what they tell us about biological visual processing. In that case, I think we need more results to learn something new, some interesting pattern across model layers or across model architectures or models trained on different objectives. And then we’d need a framework for thinking about what “non-linearity in the mapping between spaces” means. Does higher JNE mean two spaces are “less natively aligned” and thus less similar? Does the perfect “model-layer for brain region B” show low JNE or high JNE when mapped to that region? (maybe that could be answered by aligning human-to-human?).", "limitations": "yes", "rating": 4, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strengths:  \n- Well-written, clear presentation. \n- Well-motivated, strong methodological approach. \n- Intriguing results. \n\nWeaknesses:  \n-  The authors aim to understand non-linearities in biological neural processing, but the methods measure non-linearities in the mapping between spaces (model-output-space and voxel-space), and it’s not clear the measure can speak directly to non-linearities in brain responses.  \n- The work is not fully situated in the visual neuroscience literature, and what we already know about non-linear hierarchical processing in visual cortex, and hierarchical correspondence between models and brains. Authors need to clarify what have we learned about biological vision that we didn’t already know from prior work…It can’t be that non-linearities increase as you progress in the visual processing stream (already established in prior neurophysiology work; already inferred in hierarchical correspondence between models and brains). The authors can argue that JNE provides the first direct measurement of this putative model/brain correspondence, but then it needs to be more clearly demonstrated that JNE scores are directly attributable to properties of the target-space (brain responses). If the JNE scores are model-dependent, or model-layer-dependent, then the JNE scores can only be understood as a measure of the non-linearity in the mapping from model-layer-A to brain-region-B (a property of the mapping between two spaces).", "was_revised": true, "final_justification": "The quality and rigor in the methods seems clear to me. I thought that caution/clarity on the interpretation was needed, and deeper integration with related neuroscience research would increase impact of the work. The authors seem to have a solid plan for both.", "rebuttal": "We would like to thank the reviewer for the constructive comments and suggestions, as well as the appreciation of our study. Our point-by-point responses are detailed as follows.\n\n### **Response to Weaknesses**\n\n**W1:** *It's not clear the measure can speak directly to non-linearities in brain responses.*\n\n**Response:** We agree that JNE does not directly reflect neural-level nonlinearity, as fMRI BOLD signals are shaped by both neural activity and vascular dynamics. Therefore, JNE reflects a composite nonlinearity from input to BOLD response. The estimated nonlinearity captured by JNE may arise from either neural encoding or BOLD-related factors (e.g., HRF nonlinearity, vascular saturation), which cannot be fully disentangled. While BOLD signals are known to introduce complex nonlinearities, techniques such as HRF deconvolution or explicit modeling may help mitigate their influence. We thus position JNE as a practical tool for probing nonlinear encoding in fMRI settings, and will clarify its methodological scope and limitations more explicitly in the revised manuscript.\n\n**W2:** *The work is not fully situated in the visual neuroscience literature...*\n\n**Response:** We will revise the manuscript to emphasize methodological contributions and better situate our work within the visual neuroscience literature, especially regarding cortical hierarchies.\n\n### **Response to Questions**\n\n**Q1.** *The JNE score measures...*\n\n**Response:** High JNE scores arise only when the target mapping exhibits substantial nonlinearity, as supported by our theoretical and simulation-based analyses. However, due to neural-vascular coupling, JNE cannot isolate cortical nonlinearity from BOLD signal properties. We therefore position JNE as a methodological tool for analyzing nonlinear encoding models, with fMRI serving as partial validation. In the revision, we will clarify this conceptual scope and explicitly discuss the limitations of interpreting JNE as a direct measure of neural nonlinearity.\n\n**Theoretical Justification**\n\nThe JNE metric quantifies a model's sensitivity to input perturbations, defined as the standard deviation of the Jacobian norm computed across different input samples (i.e., predicted fMRI responses). By the chain rule, the complete input–output mapping in a neural encoding framework can be expressed as:\n\n$$\nI \\xrightarrow{m(\\cdot)} X \\xrightarrow{f(\\cdot)} \\hat{y}\n$$\n\nHere, $m(\\cdot)$ denotes a pretrained ANN that extracts features $X$ from the input $I$, and $f(\\cdot)$ is a linear or nonlinear encoder that maps the features into the BOLD response space.\n\nApplying the chain rule gives:\n\n$$\n\\frac{\\partial \\hat{y}}{\\partial I} = \\frac{\\partial f(X)}{\\partial X} \\cdot \\frac{\\partial X}{\\partial I}\n$$\n\nSince the feature extractor $m(\\cdot)$ is fixed, the term $\\partial X / \\partial I$ is constant across all voxels. Therefore, voxel-wise differences in input–output sensitivity are entirely determined by $\\partial f(X) / \\partial X$. From a physiological perspective, the actual brain response pathway can be approximated as:\n\n$$\nI \\xrightarrow{g(\\cdot)} \\text{Neural Activity} \\xrightarrow{z(\\cdot)} y\n$$\n\nIf the ANN representation $X = m(I)$ partially approximates the neural state (a widely supported assumption in prior neural encoding work), then $\\partial f(X) / \\partial X$ can be viewed as a first-order approximation of $\\partial z(g(I)) / \\partial I$. That is:\n\n$$\n\\frac{\\partial f(X)}{\\partial X} \\propto \\frac{\\partial z(g(I))}{\\partial I}\n$$\n\nThus, variation in the JNE score reflects the degree of nonlinear sensitivity in the composite system $z(g(\\cdot))$. JNE only increases when this system exhibits substantial input-dependent nonlinearity.\n\n**Simulations**\n\nTo validate the above derivation, we designed a structured simulation experiment:\n\n* We constructed a tensor $A \\in \\mathbb{R}^{1000 \\times 244 \\times 768}$ to simulate the ANN-based input-to-feature mapping;\n* We defined a tensor $B \\in \\mathbb{R}^{1000 \\times 768 \\times 1}$ to represent nonlinear mappings from features to response;\n* We sampled 100 scaling coefficients $\\lambda_k = \\sin(x_k)$ uniformly from $[0, 2\\pi]$, simulating amplitude-modulated responses;\n* The final output was computed as $C_k = \\texttt{einsum}(A, \\lambda_k B)$, simulating input-dependent model responses.\n\nWe then computed JNE for both $\\lambda_k B$ and $C_k$. The standardized results are as follows:\n\n| $k$                          | 1   | 2   | 3   | 4   | 5   |\n| ------------------------------ | --- | --- | --- | --- | --- |\n| **JNE($\\lambda_k B$)** | .50 | .56 | .59 | .62 | .65 |\n| **JNE($C_k$)**         | .50 | .56 | .59 | .62 | .65 |\n\nThe two curves are closely aligned, confirming that JNE reliably captures the degree of nonlinear modulation in the output space, consistent with our theoretical expectations.\n\n**Q2:** *On this same topic...*\n\n**Response:** We agree with the reviewer that the JNE score is influenced not only by intrinsic response properties of brain regions, but also by the diversity of input features and stimulus sets.\n\nTo control for such confounding factors, we used the same CLIP-ViT representations and a unified natural image set across all brain regions, ensuring that JNE comparisons were conducted under consistent input conditions. This design allows us to focus specifically on regional differences in nonlinear response behavior, rather than input-driven variability in JNE scores.\n\nAs the reviewer pointed out, changes in stimulus distribution (e.g., using only frontal faces) can affect JNE, even if the underlying brain function remains constant. This confirms that JNE reflects the input-dependent variability or output sensitivity of a region under a particular model and stimulus condition, rather than its absolute degree of nonlinearity.\n\nThis conditional dependence is analogous to how $R^2$ also varies with stimulus diversity, yet remains a valid metric for comparing brain regions under the same experimental setup. In the revised manuscript, we will further clarify this conditional and relative nature of JNE, emphasizing that it should not be interpreted as an absolute or task-invariant measure of regional nonlinearity.\n\n**Q3:** *One more question on interpretation...*\n\n**Response:** As discussed in our response to W1, Q1 and Q2, the JNE value is jointly influenced by the structure of the model's mapping and the distribution of the input data. Moreover, it is further shaped by the nonlinear characteristics of the BOLD signal itself. Therefore, JNE should not be directly interpreted as a literal measure of the brain's intrinsic nonlinearity in response or processing. We will clarify this point explicitly in the revised manuscript.\n\n**Q4:** *I found myself wondering...*\n\n**Response:** We conducted a layer-wise analysis of JNE across all CLIP-ViT layers and found that JNE values vary systematically with model depth. We did not provide these results due to space limitations. In brief, JNE values increase sharply in the first two layers and stabilize in higher layers. Mid-to-late layers yield higher JNE scores in higher-order visual cortices, while early visual areas consistently show lower JNE values across all layers—typically below the whole-brain mean. Importantly, while the absolute JNE scores differ across layers, the relative ranking of brain regions remains highly consistent. This suggests that our core findings—such as hierarchical gradients and inter-regional differences—are robust to the choice of feature layer. The reviewer raised an insightful point that is worth further efforts. In the revision, we will focus on the methodological contribition of JNE in interpreting nonlinear encoding models.\n\n**Q5:** *In theory, with a powerful...*\n\n**Response:** In our response to **Q1**, we provided a theoretical argument suggesting that JNE can serve as an approximate proxy for nonlinearity in the target (output) space. We acknowledge that JNE is not a direct metric to measure the non-linearlity of the neural cortex, due to the complex nerual-vascular coupling mechanism.  Please refer our response to W1 and Q4 for details.\n\n**Q6.** *It would be helpful if the authors...*\n\n**Response:** We apologize for placing undue emphasis on the exploratory neuroscientific findings. In response, we will revise the manuscript to make the methodological contributions of this work more explicit. The fMRI encoding experiments serve as part of the validation instead. Additionally, we will incorporate prior studies on hierarchical processing in biological vision as the reviewer suggested. While we did observe some novel  patterns about the nonlinearity (e.g., category-specific nonlinearity preferences), the validation of those findings requires additional experiments involving multiple vision-language models, multiple fMRI datasets and more sophisticated statiscal comparisons.\n\n**Q7:** *If the goal is to characterize...*\n\n**Response:** JNE is not independent of the encoding model. As addressed in **Q1**, it serves as an approximate characterization of nonlinearity in the BOLD response space, conditioned on the given model. Furthermore, **Q4** demonstrates that the observed trends in JNE hold consistently across different representation layers, indicating that our conclusions are not sensitive to the choice of a specific layer.\n\n**Q8:** *In the Discussion...*\n\n**Response:** We thank the reviewer for their attention to this point. In the manuscript, our reference to a \"bottom-up processing mechanism\" is intended to describe the hierarchical progression by which visual information is processed from lower- to higher-level cortical areas, with increasing representational complexity. Specifically, we observed a monotonic increase in JNE values from early to higher visual regions, which supports a layered transition from linear to nonlinear processing along the cortical hierarchy. We will clarify this interpretation in the revised manuscript to avoid potential conceptual ambiguity."}, {"review_id": "wQk1mYBypf", "reviewer": "Reviewer_cJqm", "summary": "This paper introduces the Jacobian-based nonlinearity evaluation (JNE), a metric for quantifying the nonlinearity of neural encoding models.  The authors first argue that the common practice of comparing the predictive accuracy ($R^2$) of linear versus nonlinear models is an unreliable method for assessing nonlinearity, especially when using pre-trained features.  The proposed JNE metric quantifies the statistical variability of a model's local linear mappings (its Jacobian matrices) across a diverse set of input stimuli.  Applying this metric to fMRI data reveals two primary contributions: first, it provides quantitative evidence for a hierarchical gradient of increasing nonlinearity from the primary visual cortex to higher-order visual areas.  Second, it demonstrates that this nonlinearity is stimulus-selective, with regions such as the fusiform face area (FFA) and the parahippocampal place area (PPA) showing the most pronounced nonlinear responses to their preferred stimulus categories (faces and scenes, respectively), thereby linking a region's computational properties to its known functional specialization.", "questions": "I suggest that you generate a large number of simulated linear and nonlinear networks, calculate the JNE score for each network and\ntreat \"linear\" and \"nonlinear\" as the two ground-truth classes and the JNE score as the classifier's output. Could you then plot an ROC curve to see how well the JNE score distinguishes between the two classes across all possible thresholds? The AUC would quantify the metric's discriminative power in this idealized environment.", "limitations": "The discussion of limitations, including the critical BOLD confound, should not be in the appendix, but this can be fixed in the camera-ready version.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "The work is well-executed and addresses a timely and important problem\n\nA major strength is the empirical demonstration that comparing the predictive accuracy ($R^2$) of linear and nonlinear models is an unreliable method for assessing nonlinearity. By showing that both model types achieve nearly identical performance when using powerful features from foundation models like CLIP-ViT, the authors effectively invalidate a common heuristic and establish a clear need for a new approach.\n\nThe JNE metric is conceptually novel. While the Jacobian matrix is a standard mathematical tool, the proposal to use the statistical dispersion of Jacobians across a sample population as a direct, non-parametric measure of nonlinearity is an insightful formulation.  Its simplicity is noteworthy when compared to more complex learnable frameworks like LinBridge, which require training an auxiliary model to decompose the Jacobian.\n\nThe JNE metric is built on the principle that a truly linear system must have an input-invariant Jacobian. The authors then rigorously validate this premise through a series of well-designed simulation experiments. The network-level validation further confirms that JNE can reliably detect and quantify nonlinearity in more complex, multi-layer systems.\n\nThe introduction clearly articulates the problem and motivates the need for the JNE metric. The four-step computational workflow of JNE is presented with mathematical precision and is well-supported by a clear diagram (Fig. 2a,c).\n\nThe most critical limitation is that the JNE metric is applied to a model of the fMRI BOLD signal, not directly to neural activity. The BOLD signal is an indirect and sluggish hemodynamic measure that has its complex nonlinear relationship with the underlying neural events. These nonlinearities can be vascular in origin (e.g., blood flow refractory effects) or related to neural phenomena like adaptation, and are highly sensitive to experimental design, such as stimulus timing.\n\nThe paper notes that JNE lacks a defined threshold to distinguish a \"linear\" from a \"nonlinear\" region in absolute terms. This restricts the metric's utility to relative comparisons. It could, however, be quantified with simulation experiments and an ROC analysis, at least in an idealized setting.\n\nWhile the method is novel, the neuroscientific findings could be viewed as incremental. The discussion of limitations, including the critical BOLD confound, should not be in the appendix, but this can be fixed in the camera-ready version.", "was_revised": true, "final_justification": "I'm happy with the author's response and retain my score", "rebuttal": "We would like to thank the reviewer for the constructive comments and suggestions, as well as the appreciation of our study. Our point-by-point responses are detailed as follows.\n\n### **Response to Weaknesses**\n\n**W1:** *The most critical limitation is that the JNE metric is applied to a model of the fMRI BOLD signal, not directly to neural activity. The BOLD signal is an indirect and sluggish hemodynamic measure that has its complex nonlinear relationship with the underlying neural events. These nonlinearities can be vascular in origin (e.g., blood flow refractory effects) or related to neural phenomena like adaptation, and are highly sensitive to experimental design, such as stimulus timing.*\n\n**Response:** We appreciate the reviewer's comments. As noted, the BOLD signal is an indirect and temporally filtered measurement of neural activity, shaped by complex and nonlinear hemodynamic processes.\n\nAccordingly, the scope of our analysis is explicitly limited to the nonlinear structure of BOLD signals as observed via fMRI. The goal of JNE is to quantify the overall degree of nonlinearity in the mapping from input stimuli to BOLD responses, within the model–observation space. We do not claim that JNE directly reflects neuronal-level nonlinearity. Rather, it should be viewed as a tool for quantifying the explanatory behavior of the encoding model.\n\nWe also acknowledge that JNE does not distinguish whether this nonlinearity arises from neuronal or vascular mechanisms. Thus, JNE should be understood as an empirical estimate of the input–output nonlinearity of the composite function, whose interpretability is inherently constrained by the limits of fMRI as a measurement modality. In future work, we plan to incorporate HRF deconvolution strategies to attempt disentangling neural and vascular sources of nonlinearity. Furthermore, we will validate the generalizability of JNE to neural activity itself using data with higher temporal resolution (e.g., ECoG or MEG).\n\nWe will carefully revise the relevant sections to ensure that \"BOLD nonlinearity\" and \"neural nonlinearity\" are clearly distinguished, and that all terminology is used with precision.\n\n**W2:** *The paper notes that JNE lacks a defined threshold to distinguish a \"linear\" from a \"nonlinear\" region in absolute terms. This restricts the metric's utility to relative comparisons. It could, however, be quantified with simulation experiments and an ROC analysis, at least in an idealized setting.*\n\n**Response:** We sincerely thank the reviewer for the valuable suggestion. In its current form, JNE is primarily used for relative comparison of nonlinearity levels across model architectures or brain regions in empirical analyses, and that a clear absolute decision boundary between \"linear\" and \"nonlinear\" systems has not yet been established.\n\nFollow the insightful suggestion in **Q4**, we have generated 382,340 simulated linear and nonlinear networks and performed a binary classification based on the JNE metric. The classifier achieved an ROC of 1, with linear neural networks exhibited 0 JNE values while nonlinear networks  were with non-zero JNE values. However, simulations under more complex conditions such as non-uniform distribution and data noise have not been accounted for due to time limit.\n\nSpecifically, we have extended our original simulation framework (Section 3.1: *Network-level validation via ANN-based simulation*) to systematically explore the structural space of artificial networks, introducing \"layer number\" as a controlled variable $$, and jointly examining the following factors:\n\n* **Activation function type**: None, ReLU, LeakyReLU, GELU, Swish;\n* **Hidden layer width**: 16, 32, 64, 128, 256, 512, 1024;\n* **Residual connection configuration**: on / off;\n* **Activation function positions**: all $2^L$ combinations of insertion patterns per layer, where $L$ is the number of layers.\n\nBy exhaustively combining these configurations, we generated 382,340 structurally diverse networks. For each network, we performed forward passes using a fixed input distribution $\\mathcal{U}(-10, 10)$, computed the Jacobian between input and output, and extracted the JNE value.\n\nTo evaluate the discriminative power of JNE, we formulated a binary classification task:\n\n* Networks with **no nonlinearities** (activation type = None, or all insert positions set to 0) were labeled as **linear** (label = 0);\n* All other networks were labeled as **nonlinear** (label = 1).\n\n**Simulation Results:** Under these idealized conditions, JNE demonstrated excellent discriminative performance, with an AUC of 1.0.\n\n**W3:** *While the method is novel, the neuroscientific findings could be viewed as incremental. The discussion of limitations, including the critical BOLD confound, should not be in the appendix, but this can be fixed in the camera-ready version.*\n\n**Response:** We thank the reviewer's appreciation of the novely of our study. The core contribution of this study is introducing a novel metric—JNE—for quantifying the explanatory power of nonlinear neural encoding models. The incremental neuroscientific findings serve as part of the validation of the proposed method. Accordingly, we will explicitly strengthen the presentation of JNE as a methodological innovation, and move the discussion of limitations to the main text.\n\n### **Response to Questions**\n\n**Q1:** *I suggest that you generate a large number of simulated linear and nonlinear networks, calculate the JNE score for each network and treat \"linear\" and \"nonlinear\" as the two ground-truth classes and the JNE score as the classifier's output. Could you then plot an ROC curve to see how well the JNE score distinguishes between the two classes across all possible thresholds? The AUC would quantify the metric's discriminative power in this idealized environment.*\n\n**Response:** Following the reviewer's suggestion, we expanded our simulation experiments by generating over 380,000 simulated linear and nonlinear neural networks with diverse architectures. We then constructed a binary classification task based on the JNE metric. We found that JNE exhibited strong discriminative power, achieving an AUC of 1.0 under this controlled setting. Simulated linear networks exhibit 0 JNE value, and nonlinear networks are with non-zero JNE values. However, simulations under more complex conditions such as non-uniform distribution and data noise have not been accounted for due to time limit. Please refer to our response to **W2** for details.\n\n### **Response to Limitations**\n\n**L1:** *The discussion of limitations, including the critical BOLD confound, should not be in the appendix, but this can be fixed in the camera-ready version.*\n\n**Response:** We will move the discussion of limitations to the main text."}]

Author rebuttal:

Meta-review: The reviewers and I agree that this is an interesting and novel work that is thorough and clearly presented. It suggests a useful perspective on how to quantify the nonlinearity in the (mapping from model representations to) neural representations, and demonstrates through simulations and analyses of real data that this metric behaves reasonably.

I find myself agreeing with the thoughtful interpretation of reviewer VQtJ about what nonlinearity is actually being quantified here (viz., the mapping itself), so I'd encourage the authors to take that perspective into account (as well as the remaining comments e.g. moving the limitations to the main text as promised), when preparing the camera ready.
Decision: Accept (spotlight)