Multidimensional Feature Tuning in Category Selective Areas of Human Visual Cortex.

1. # %% [markdown]
2. # # Multidimensional feature tuning in category‑selective areas of human visual cortex
4. # %%
5. from scripts.run_prepro import run_prepro
6. from scripts.run_bcv import run_bcv
7. from scripts.run_bnmf import run_bnmf
8. from scripts.run_interpret import run_interpret
9. from scripts.run_encoding import run_encoding
11. from roidims.bnmf import compute_evar_consistent, compute_sim_dims_groups
12. from roidims.encoding import compute_corr_H_beta, selectivity_vs_sparseness
14. from roidims.plotting import (
15. fig_consistency,
16. fig_top_imgs,
17. fig_dprime_vs_beta,
18. fig_wordclouds,
19. fig_flatmaps_encoding_r2,
20. fig_r2_bar,
21. fig_flatmaps_encoding_betas,
22. suppfig_k_optim,
23. suppfig_consistency,
24. suppfig_sim_dims,
25. suppfig_sparseness,
26. suppfig_noiseceilings,
27. )
29. # %% [markdown]
30. # ## Background
31. # What functional organization underlies both category-selective areas and continuous feature maps in high-level visual cortex?
32. #
33. # We addressed this question by applying a data-driven decomposition to fMRI responses from face-, body-, and scene-selective areas, revealing that both views reflect a common principle: multidimensional tuning that is clustered within areas yet distributed across cortex.
35. # %% [markdown]
36. # ## fMRI data
37. # We used a subset of the Natural Scenes Dataset (NSD; Allen et al., 2022) including fMRI responses (7T) from four participants, each viewing 10,000 natural images.
38. #
39. # We defined FFA, PPA, and EBA using an independent functional localizer (fLoc; Stigliani et al., 2015), selecting voxels with a strong category preference (t > 2) and sufficient reliability (SNR > 0.2).
41. # %%
42. # Specify subjetcs and ROIs
43. subjects = ["subj01", "subj02", "subj05", "subj07"]
44. rois = ["FFA", "EBA", "PPA"]
46. # %% [markdown]
47. # ## 0. Preprocessing
48. # Single-trial responses were preprocessed with GLMdenoise and ridge regression, z-scored within sessions, averaged across repetitions, and split 70/30 into training and test sets.
50. # %%
51. # Preprocessing
52. run_prepro()
54. # %% [markdown]
55. # ## 1. Data-driven voxel decomposition
56. # First, we applied non-negative matrix factorization to reveal which underlying representational dimensions capture the fMRI activity patterns.
58. # %% [markdown]
59. # ### Background: Bayesian non-negative matrix factorization (BNMF)
60. # We decomposed voxel responses from each ROI into representational dimensions using BNMF (Schmidt et al., 2009; following Khosla et al., 2022).
61. #
62. # BNMF factorizes the image × voxel response matrix into a response matrix W (image × dimension) and a weight matrix H (dimension × voxel). The non-negativity constraint encourages sparse, additive, part-based dimensions, allowing voxels to participate in multiple dimensions simultaneously. The Bayesian extension adds exponential priors on W and H and uses Gibbs sampling to obtain robust posterior estimates.
63. #
64. # We ran 3,000 iterations (1,000 burn-in, every 5th sample retained) and baseline-shifted responses to enforce non-negativity, using the training minimum for both sets to prevent data leakage.
66. # %% [markdown]
67. # ### 1.1. Estimating the optimal dimensionality
68. # We determined the optimal number of dimensions k* for each participant and ROI using bi-cross-validation (Owen & Perry, 2009).
69. #
70. # Bi-cross-validation partitions the data into four blocks. BNMF is trained on one block, and pseudo-inverse operations on two others are used to predict the held-out block. This was repeated 5 times across ranks 1–100 (step size 3), and k* was chosen as the rank minimizing average test error.
72. # %%
73. # Bi-cross-validation
74. run_bcv(subjects, rois, k_min=1, k_max=100, k_steps=1, n_perms=5)
76. # %%
77. # Supp Fig 1: Optimal dimensionality
78. suppfig_k_optim(subjects, rois)
80. # %% [markdown]
81. # ### 1.2. Finding reliable dimensions
82. # To obtain reliable and generalizable dimensions, we used a two-step consensus approach:
83. #
84. # First, we ran BNMF 100 times with k* but random initializations per ROI, removed outlier runs, and aggregated stable dimensions via k-medoids clustering (Kotliar et al., 2019).
85. #
86. # Second, we matched dimensions across participants using pairwise correlations and a greedy selection procedure, retaining only dimensions consistent across all four participants (r > 0.3; Khosla et al., 2022).
87. #
88. # This yielded 8–20 consensus dimensions per ROI. To test the generalizability to unseen images, we projected the test set into the learned embedding using non-negative least squares regression.
90. # %%
91. # Specify inter-subject consistency threshold
92. r_thresh = 0.3
94. # %%
95. # Consensus approach
96. run_bnmf(subjects, rois, n_runs=100, r_thresh=r_thresh)
98. # %%
99. # Compute explained variance
100. compute_evar_consistent(subjects, rois)
102. # %%
103. # Fig 2a: Consistency
104. fig_consistency(rois, r_thresh)
106. # %%
107. # Supp Fig 2: Mean inter-participant consistency of all dimensions
108. suppfig_consistency(rois, r_thresh)
110. # %%
111. # Fig 2b: Top images combined across subjects
112. fig_top_imgs(subjects, rois)
114. # %%
115. # Supp Fig 3: Mean similarity matrices
116. suppfig_sim_dims(subjects, rois)
118. # %%
119. # Representational similarity of all dimensions
120. compute_sim_dims_groups(subjects, rois)
122. # %% [markdown]
123. # ## 2. Data-driven dimension labeling
124. # To interpret each dimension, we used a behavioral online experiment. We collected data from N=38 participants, who viewed collages of the 20 highest-scoring images per dimension and provided 3–5 short labels describing what the images had in common. Labels were translated from German, cleaned of generic entries, and standardized.
125. #
126. # We then used CLIP-ViT-L/14 (Radford et al., 2021) to select the most representative label per dimension from the participant-generated candidate set. For each label, we computed the cosine similarity between its text embedding and all image embeddings, weighted by each dimension's response profile, then z-scored across labels within each dimension to highlight selectively associated labels. Top-ranked labels were visualized as word clouds, yielding concise, semantically meaningful descriptors validated across the full image set.
128. # %% [markdown]
129. # ### 2.1. Finding the most representative labels
131. # %%
132. # CLIP-based interpretation
133. run_interpret(subjects, rois)
135. # %%
136. # Fig 3: Top labels for all subjects
137. fig_wordclouds(subjects, rois)
139. # %% [markdown]
140. # ## 3. Testing the relationship between multidimensional tuning and category selectivity
141. # To quantify category selectivity for the ROI's preferred category, we used a category selectivity index (d'), comparing mean responses to a preferred category against all others, weighted by response variance.
142. #
143. # To measure the extent of multidimensional tuning, we computed a sparseness index (Hoyer, 2004) based on the normalized L1/L2-norm ratio of each voxel's dimension weight profile, ranging from 0 (tuned equally to all dimensions) to 1 (tuned to a single dimension).
145. # %%
146. # Fig 3: Correlation between dimension tuning (a.u.) and category selectivity (d')
147. fig_dprime_vs_beta(subjects, rois)
149. # %%
150. # Plot sparseness for different selectivity bins
151. suppfig_sparseness(subjects, rois)
153. # %%
154. # Print correlation
155. selectivity_vs_sparseness(subjects, rois)
157. # %% [markdown]
158. # ## 4. Voxel-wise encoding models
159. # Finally, we fit voxel-wise encoding models to predict fMRI responses to held-out images from the learned dimensions, revealing the cortical topography of this organization.
160. #
161. # We used fractional ridge regression with cross-validation to avoid overfitting, projected test images into the learned BNMF space with non-negative least squares, and evaluated prediction performance on held-out data. Significance was assessed with permutation tests and FDR correction, and noise ceilings were computed to benchmark model performance.
163. # %%
164. # Fractional ridge regression
165. run_encoding(subjects, rois, n_folds=10, n_perms=3000)
167. # %%
168. # Print correlation between ROI coefficient weights
169. compute_corr_H_beta(subjects, rois)
171. # %%
172. # Fig 4a: Prediction performance maps
173. fig_flatmaps_encoding_r2
175. # %%
176. # Fig 4b: Prediction performance across ROIs
177. fig_r2_bar(subjects, rois)
179. # %%
180. # Supp Fig 6: Voxel-wise prediction performance vs. noise ceiling estimate
181. suppfig_noiseceilings(subjects, rois)
183. # %%
184. # Fig 5: Individual dimension tuning maps
185. fig_flatmaps_encoding_betas(subjects, rois)

main_analyses.ipynb at commit 05bce4d, no license · at the source