OSCR

Functional organization and natural scene responses across mouse visual cortical areas revealed with encoding manifolds.

Code ↔ Paper

18 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 18 matches
  1. [1] § Materials and methods › Data analysis › Similarity kernel and diffusion maps. ↔ allen-data-analysis/encoding-manifold-VISp.ipynb, lines 314–324 · score 0.85 · Laplace Beltrami approximation, IAN weighted graph, Diffusion maps, Laplacian, outliers, embedding
  2. [2] § Materials and methods › Data analysis › Similarity kernel and diffusion maps. ↔ encoding-manifold/encoding-manifold.ipynb, lines 263–272 · score 0.84 · Laplace Beltrami approximation, IAN weighted graph, Diffusion maps, Laplacian, embedding, matrix
  3. [3] § Materials and methods › Data analysis › Additional metrics. ↔ allen-data-analysis/read_spike_data-drifting gratings.ipynb, lines 340–404 · score 0.79 · VISam, VISrl, drifting gratings, temporal frequencies, VISal, VISpm
  4. [4] § Materials and methods › Data preprocessing ↔ allen-data-analysis/read_spike_data-drifting gratings.ipynb, lines 340–404 · score 0.75 · VISam, VISrl, drifting grating, VISal, VISpm, ISI
  5. [5] § Results ↔ allen-data-analysis/read_spike_data-static-gratings.ipynb, lines 622–696 · score 0.71 · VISam, VISrl, VISal, VISpm, static gratings, phases
  6. [6] § Materials and methods › Dataset ↔ allen-data-analysis/read_spike_data-static-gratings.ipynb, lines 418–439 · score 0.69 · visual space prior, stimulus warping, monitor, mouse
  7. [7] § Materials and methods › Data analysis › Additional metrics. ↔ allen-data-analysis/read_spike_data-static-gratings.ipynb, lines 458–534 · score 0.66 · VISam, VISrl, VISal, VISpm, metrics, spiking
  8. [8] § Materials and methods › Data analysis › Neural encoding space. ↔ CNNs/build-CNN-manifold-resnet50_block3.ipynb, lines 10–152 · score 0.63 · tensor components, neural matrices, split, lowest, reconstruction, error
  9. [9] § Materials and methods › Data analysis › Tensor decomposition. ↔ permuted-decomposition/matlab/run_permcp.m, lines 1–95 · score 0.62 · direct optimization, Tensor Toolbox, modified, permutation
  10. [10] § Materials and methods › Data analysis › Natural scene vs. static grating selectivity ratio. ↔ allen-data-analysis/read_spike_data-static-gratings.ipynb, lines 622–696 · score 0.61 · Brain Observatory, static gratings, stimulus class, phase, orientation, spikes
  11. [11] § Materials and methods › Data preprocessing ↔ allen-data-analysis/read_spike_data-static-gratings.ipynb, lines 281–391 · score 0.61 · Improved Sheather Jones, static gratings, bandwidth, algorithm, kernel, smoothed
  12. [12] § Materials and methods › Data preprocessing ↔ allen-data-analysis/read_spike_data-drifting gratings.ipynb, lines 262–318 · score 0.60 · Improved Sheather Jones, drifting gratings, bandwidth, algorithm, kernel, smoothed
  13. [13] § Materials and methods › Data analysis › Neural encoding space. ↔ CNNs/build-CNN-manifold-resnet50_block3.ipynb, lines 10–152 · score 0.59 · neural matrix, stimulus response, linear, vectors, product, encoding
  14. [14] § Materials and methods › Data analysis › Natural scene filtering. ↔ allen-data-analysis/filtering-natural-scenes/ffttools.py, lines 621–687 · score 0.57 · band pass, inner, radius, outer, Filtered, natural scenes
  15. [15] § Materials and methods › Data analysis › Neural encoding space. ↔ allen-data-analysis/encoding-manifold-VISp.ipynb, lines 56–154 · score 0.56 · neural matrix, stimulus response, linear, vectors, encoding, neuron
  16. [16] § Materials and methods › Data analysis › Neural encoding space. ↔ permuted-decomposition/choosing-n-of-components.ipynb, lines 213–225 · score 0.56 · lowest reconstruction error, components, encoding
  17. [17] § Materials and methods › Data analysis › Natural scene filtering. ↔ allen-data-analysis/filtering-natural-scenes/filtering-natural-scenes.ipynb, lines 72–134 · score 0.55 · band pass, Filtered, cpd, cropped, natural scenes
  18. [18] § Materials and methods › Data analysis › Neural encoding space. ↔ allen-data-analysis/encoding-manifold-VISp.ipynb, lines 56–154 · score 0.53 · neural matrix, neural factors, reconstruction, magnitudes, component, tensor

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Jupyter notebook · 696 lines · 24 KB · BSD-2-Clause · 5 matches

  1. # %%
  2. import os
  3. # https://allensdk.readthedocs.io/en/latest/visual_coding_neuropixels.html#
  4. #https://allensdk.readthedocs.io/en/latest/_static/examples/nb/ecephys_quickstart.html
  5. from ipywidgets import FloatProgress
  6. import numpy as np
  7. import pandas as pd
  8. import matplotlib.pyplot as plt
  9. import pickle
  10. from allensdk.brain_observatory.ecephys.ecephys_project_cache import EcephysProjectCache
  12. # %% [markdown]
  13. # #### ftns
  15. # %%
  18. #https://allensdk.readthedocs.io/en/latest/_static/examples/nb/ecephys_session.html#Stimulus-presentations
  19. from itertools import product
  20. from itertools import product
  21. def get_spike_trains(my_units, STIM_CLASS, PARAM_NAMES, prestim_len, bin_size, my_trial_len=None):
  23. trials_info = session.get_stimulus_table(STIM_CLASS)
  25. if trials_info.index.size == 0:
  26. return {}
  28. trial_len = round(np.mean(trials_info['duration'].values),2) #sec
  29. if my_trial_len is not None:
  30. assert my_trial_len <= trial_len
  31. trial_len = my_trial_len
  32. print('trial_len',trial_len)
  33. tot_len = trial_len + prestim_len
  35. time_bin_edges = np.linspace(-prestim_len, trial_len, int(tot_len/bin_size))
  36. NBINS = len(time_bin_edges)-1
  38. binarize = False
  39. if bin_size < .001:
  40. binarize = True
  42. print('NBINS',NBINS)
  44. # look at responses to a certain type of gratings
  45. stim_params = {pname:sorted(set(trials_info[pname].unique()).difference(['null'])) for pname in PARAM_NAMES}
  46. print(f'{STIM_CLASS} params:', stim_params)
  48. #get all param combinations
  49. paramvals_tuples = list(product(*[stim_params[pname] for pname in PARAM_NAMES]))
  51. all_trains = {}
  52. for ii,ui in enumerate(my_units):
  53. if ii % 5 == 0: print(ii,end=' ',flush=True)
  54. all_trains[ui] = {}
  56. for pvals_tup in paramvals_tuples:
  57. row_filter = np.prod(np.stack([trials_info[pname].values == pval for pname,pval in zip(PARAM_NAMES,pvals_tup)],axis=0),axis=0).astype('bool')
  59. sids = trials_info[row_filter].index.values
  60. Ntrials = len(sids)
  61. # print(f'{pvals_tup}, {Ntrials=}')
  63. #TODO compare speed against reading spk times directly: https://allensdk.readthedocs.io/en/latest/_static/examples/nb/ecephys_optotagging.html
  64. spike_counts_da = session.presentationwise_spike_counts(
  65. bin_edges=time_bin_edges,
  66. stimulus_presentation_ids=sids,
  67. unit_ids=[ui],
  68. binarize=binarize
  69. )
  70. spike_counts_da = np.squeeze(spike_counts_da.values)
  72. # print(spike_counts_da.shape,spike_counts_da.max(),spike_counts_da.sum())
  73. #re-convert to spike times
  75. unit_trains = []
  76. for triali in range(Ntrials):
  77. train = np.flatnonzero(spike_counts_da[triali]).astype('float32')
  78. unit_trains.append(train * bin_size * 1000) #convert bin number to time (ms)
  79. all_trains[ui][pvals_tup] = unit_trains
  80. print()
  81. return all_trains
  82. def get_spike_trains_pref_phase(my_units, STIM_CLASS, PARAM_NAMES, prestim_len, bin_size):
  84. trials_info = session.get_stimulus_table(STIM_CLASS)
  86. if trials_info.index.size == 0:
  87. return {}
  89. trial_len = round(np.mean(trials_info['duration'].values),2) #sec
  90. print('trial_len',trial_len)
  91. tot_len = trial_len + prestim_len
  93. time_bin_edges = np.linspace(-prestim_len, trial_len, int(tot_len/bin_size))
  94. NBINS = len(time_bin_edges)-1
  96. binarize = False
  97. if bin_size < .001:
  98. binarize = True
  100. print('NBINS',NBINS)
  102. # look at responses to a certain type of gratings
  103. stim_params = {pname:sorted(set(trials_info[pname].unique()).difference(['null'])) for pname in PARAM_NAMES}
  104. print(f'{STIM_CLASS} params:', stim_params)
  106. #get all param combinations
  107. paramvals_tuples = list(product(*[stim_params[pname] for pname in PARAM_NAMES]))
  109. pref_phases = cache.get_unit_analysis_metrics_for_session(session_id)['pref_phase_sg']
  111. all_trains = {}
  112. for ii,ui in enumerate(my_units):
  113. print(ii,end=' ')
  114. all_trains[ui] = {}
  116. for pvals_tup in paramvals_tuples:
  117. #append this unit's pref phase to the stimulus params tuple
  118. assert ui in pref_phases.index
  119. pvals_tup_with_phase = pvals_tup + (str(pref_phases.loc[ui]),)#must convert phase to str (!?)
  120. PARAM_NAMES_with_phase = PARAM_NAMES + ['phase']
  121. row_filter = np.prod(np.stack(
  122. [trials_info[pname].values == pval for pname,pval in \
  123. zip(PARAM_NAMES_with_phase,pvals_tup_with_phase)],axis=0),axis=0).astype('bool')
  125. sids = trials_info[row_filter].index.values
  126. Ntrials = len(sids)
  127. # print(f'{pvals_tup}, {Ntrials=}')
  129. spike_counts_da = session.presentationwise_spike_counts(
  130. bin_edges=time_bin_edges,
  131. stimulus_presentation_ids=sids,
  132. unit_ids=[ui],
  133. binarize=binarize
  134. )
  135. spike_counts_da = np.squeeze(spike_counts_da.values)
  137. # print(spike_counts_da.shape,spike_counts_da.max(),spike_counts_da.sum())
  138. #re-convert to spike times
  140. unit_trains = []
  141. for triali in range(Ntrials):
  142. train = np.flatnonzero(spike_counts_da[triali]).astype('float32')
  143. unit_trains.append(train * bin_size * 1000) #convert bin number to time (ms)
  144. all_trains[ui][pvals_tup] = unit_trains
  145. print()
  146. return all_trains
  150. def get_train_dicts(all_trains_uid, tf, dirs, trial_len_ms, prestim_len_ms):
  151. trials_traindict = {}
  152. ISI_Nspks = {}
  153. full_traindict = {}
  154. for d in dirs:
  155. ptup = (tf, d,)
  156. assert ptup in all_trains_uid
  158. ISI_Nspks[d] = []
  159. new_trains = []
  160. for train in all_trains_uid[ptup]:
  161. if len(train) == 0:
  162. new_trains.append(np.array([]))
  163. ISI_Nspks[d].append(0)
  164. continue
  165. assert max(train) < trial_len_ms + prestim_len_ms
  168. new_train = train - prestim_len_ms
  170. new_trains.append(new_train[new_train >= 0])
  171. ISI_Nspks[d].append(new_train[new_train < 0].size)
  173. trials_traindict[d] = new_trains
  174. ISI_Nspks[d] = np.asarray(ISI_Nspks[d])
  175. full_traindict[d] = all_trains_uid[ptup]
  176. return full_traindict, trials_traindict, ISI_Nspks
  178. # %%
  179. import scipy as sp
  180. import warnings
  181. def computeResponseStats(traindict, ISI_Nspks, stats_ISI_len, trial_len, verbose=False):
  183. """Runs statistical tests to compare firing rates between the ISI and a given stimulus,
  184. for any period within the stimulus trial with the same length as the ISI.
  185. It performs two comparisons against the ISI FR: one using the maximum FR found for that
  186. interval length across stimulus trials; and another using the minimum FR.
  188. The following one-sided tests are run:
  189. Mann-Whitney U
  190. https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.mannwhitneyu.html
  192. and Wilcoxon:
  193. https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.wilcoxon.html
  195. ----------------
  196. Arguments:
  197. traindict: dict, {stimulus_direction: list of spike time arrays (ms), one per trial}
  198. ISI_Nspks: dict, {stimulus_direction: list of spike counts, one per trial}
  199. stats_ISI_len: float, length of the ISI interval, in secs, used to compute ISI_Nspks
  200. trial_len: float, total length of each trial, in secs
  202. ----------------
  203. Returns:
  204. stats_results: dict, for each of 'min-interval' and 'max-interval', contains a dict
  205. containing, for each stimulus direction, the p-values found for each tests, as well
  206. as the FRs for the stimulus and the ISI
  207. """
  208. mydirs = traindict.keys()
  209. stats_results = {}
  211. interval_Nspks_for_stats = {'min':{}, 'max':{}}
  212. for d in mydirs:
  214. #combine all trains
  215. data = np.concatenate(traindict[d])
  216. counts,bins = np.histogram(data,np.arange(0,trial_len+stats_ISI_len,stats_ISI_len))
  218. maxfr = -1
  219. minfr = np.inf
  220. for i in range(0,counts.size):
  221. fr = counts[i]
  222. if fr > maxfr:
  223. maxi = i
  224. maxfr = fr
  225. if fr < minfr:
  226. mini = i
  227. minfr = fr
  230. for interval_type,i_,fr_ in [('min',mini,minfr), ('max',maxi,maxfr)]:
  232. interval_Nspks_for_stats[interval_type][d] = []
  234. for train in traindict[d]:
  235. spks_within_interval = train[(train >= stats_ISI_len*i_) & (train < stats_ISI_len*(i_+1))]
  236. interval_Nspks_for_stats[interval_type][d].append( spks_within_interval.size )
  238. interval_Nspks_for_stats[interval_type][d] = np.array(interval_Nspks_for_stats[interval_type][d])
  239. # print(d,interval_type,interval_Nspks_for_stats[interval_type][d])
  240. assert len(interval_Nspks_for_stats[interval_type][d]) == len(traindict[d])
  241. assert interval_Nspks_for_stats[interval_type][d].sum() == fr_
  244. for trial_Nspks_for_stats_,pval_type,alternative in [
  245. (interval_Nspks_for_stats['max'],'maxinterval-pval','greater'),
  246. (interval_Nspks_for_stats['min'],'mininterval-pval','less')]:
  247. stats_results[pval_type] = {}
  250. for d in mydirs:
  251. dir_grayspks = ISI_Nspks[d]
  252. dir_stimspks = trial_Nspks_for_stats_[d]
  253. assert len(dir_stimspks) == len(dir_grayspks)
  254. if len(dir_stimspks) * len(dir_grayspks) == 0:
  255. stats_results[pval_type][d] = (np.inf,0,0)
  256. stats_results[pval_type][d] = {'MANNWHITNEY':np.inf,'WILCOXON':np.inf,'isiFR':0,'stimFR':0}
  257. if verbose: print(f'{d} no spikes')
  258. continue
  260. try:
  261. with warnings.catch_warnings():
  262. warnings.simplefilter("ignore",category=RuntimeWarning)
  263. _, pval = sp.stats.mannwhitneyu(dir_stimspks, dir_grayspks, alternative=alternative)
  264. except:
  265. pval = np.inf
  267. try:
  268. with warnings.catch_warnings():
  269. warnings.simplefilter("ignore",category=RuntimeWarning)
  270. warnings.simplefilter("ignore",category=UserWarning)
  271. _, wpval = sp.stats.wilcoxon(dir_stimspks, dir_grayspks, alternative=alternative)
  272. except:
  273. wpval = np.inf
  276. stats_results[pval_type][d] = {'MANNWHITNEY':pval,'WILCOXON':wpval,'isiFR':dir_grayspks.mean()/stats_ISI_len,'stimFR':dir_stimspks.mean()/stats_ISI_len}
  277. if verbose: print(f'{d} {pval_type} ({pval:.3f},{wpval:.3f}), gray={dir_grayspks.mean()/stats_ISI_len:.2f}, stim={dir_stimspks.mean()/stats_ISI_len:.2f}')
  279. return stats_results
  281. # %%
  282. import warnings
  283. from KDEpy import FFTKDE
  284. def getResponseCurve(train_dict, total_trial_len, bw=None, samp_interval=1, MINBW=10, MAXBW=50):
  286. """Computes smooth trial-averaged response to a stim in all directions from spike trains
  287. using a kernel density estimator."""
  289. ts = np.arange(0,total_trial_len+samp_interval,samp_interval)
  290. x_ts = .5*(ts[:-1]+ts[1:]) #sample at the midpoints between sampling intervals
  292. all_ISJs = []
  295. fftkde = None
  296. if bw is None:
  297. # if no pre-specified kernel bandwidth,
  298. # auto-estimate within range [MINBW, MAXBW]
  299. bw, fftkde = fitSmoothingKernelBandwidth(train_dict, total_trial_len)
  300. if bw is not None:
  301. bw = max(min(MAXBW,bw),MINBW)
  302. fftkde.bw = bw
  303. else:
  304. bw = MAXBW
  306. for di,d in enumerate(sorted(train_dict)):
  308. full_train = []
  309. for triali,train in enumerate(train_dict[d]):
  310. if train.size == 0: continue
  311. assert max(train) < total_trial_len
  312. full_train += list(train)
  314. data = np.array(full_train)
  315. data.sort()
  316. n = data.size
  317. if fftkde is None:
  318. fftkde = FFTKDE(kernel='gaussian', bw=bw) #initialize kernel density estimator
  320. try:
  321. fftkde = fftkde.fit(data)
  322. #extend sampling one unit before and after trial
  323. ext_x_ts = np.r_[x_ts[0]-samp_interval,x_ts,x_ts[-1]+samp_interval]
  324. #then crop after applying the kernel
  325. y = fftkde.evaluate(ext_x_ts)[1:-1]
  327. except:
  328. #error occurs in the rare cases when there are zero spikes. print out to double-check
  329. # print(f'fftkde failed: {n} data points')
  330. y = np.zeros_like(x_ts)
  333. all_ISJs.append(y*n*1000/len(train_dict[d])) #convert density to spks/sec
  334. return np.array(all_ISJs), x_ts
  335. def fitSmoothingKernelBandwidth(full_traindict, total_trial_len):
  336. """Fits a spike smoothing kernel to spike train data using
  337. the improved Sheather-Jones (ISJ) algorithm:
  339. Z. I. Botev, J. F. Grotowski, and D. P. Kroese.
  340. “Kernel density estimation via diffusion.”
  341. Annals of Statistics, Volume 38, Number 5, pp. 2916-2957, 2010.
  342. https://arxiv.org/pdf/1011.2602.pdf
  344. (see https://kdepy.readthedocs.io/en/latest/index.html
  345. for more information on this implementation)
  348. ---------------
  349. Arguments:
  350. full_traindict: dict, {stimulus_direction: list of spike time arrays, one per trial}
  351. The bandwidth is computed for the stimulus direction that elicited the most spikes.
  352. total_trial_len: float or int, total length of a trial used for the trains; must
  353. use the same time unit as the spike times in `full_traindict`
  355. ---------------
  356. Returns:
  357. opt_bw: float, optimal bandwidth found
  358. fftkde: object, the fitted fftkde object, to be reused when evaluating the kernel
  359. """
  362. # 1) use stimulus direction with max n of spks to estimate optimal bandwidth
  363. maxn = -1
  364. for d, trains in full_traindict.items():
  366. data = np.concatenate(trains)
  367. assert max(data) < total_trial_len
  368. data.sort()
  369. n = data.size
  370. if n > maxn:
  371. maxn = n
  372. bestd = d
  374. data = np.concatenate(full_traindict[bestd])
  376. # 2) fit bw
  378. fftkde = FFTKDE(kernel='gaussian', bw='ISJ')
  380. try:
  381. with warnings.catch_warnings():
  382. warnings.simplefilter("ignore",category=RuntimeWarning)
  383. fftkde = fftkde.fit(data)
  384. opt_bw = fftkde.bw
  386. except:
  387. #print(f'fftkde failed: {n} data points')
  388. opt_bw = None
  389. fftkde = None
  391. return opt_bw, fftkde #return fftkde object as well, to avoid refitting
  393. # %% [markdown]
  394. # #### load cache
  396. # %%
  397. # Example cache directory path, it determines where downloaded data will be stored
  398. output_dir = 'ecephys_cache_dir/'
  400. # this path determines where downloaded data will be stored
  401. manifest_path = os.path.join(output_dir, "manifest.json")
  403. cache = EcephysProjectCache.from_warehouse(manifest=manifest_path)
  405. print(cache.get_all_session_types())
  407. # %%
  408. sessions = cache.get_session_table()
  409. brain_observatory_type_sessions = sessions[sessions["session_type"] == "brain_observatory_1.1"]
  410. print(len(brain_observatory_type_sessions))
  411. brain_observatory_type_sessions.head()
  413. # %%
  414. from glob import glob
  415. my_session_ids = [int(dirname.split('session_')[1]) for dirname in glob('ecephys_cache_dir/session_*')]
  416. print(my_session_ids)
  418. # %% [markdown]
  419. # ### get gratings imgs
  421. # %%
  422. from allensdk.brain_observatory.stimulus_info import get_spatial_grating
  423. #aspect = 1920/1200
  424. aspect = 1174/918 #same as nat scene
  425. # 21.93" wide monitor positioned at 15 cm away from the mouse's right eye and
  426. # spanned 120° x 95° of visual space prior to stimulus warping
  427. #orig_px_per_deg = 1920/120 #tot_pixels/tot_degs
  428. height = 918//2
  429. width = aspect * height
  431. grat_imgs = []
  432. for cycle_per_deg in [.02, .04, .08, 0.16, 0.32]:
  433. for ori in [0,30,60,90,120,150]:
  434. for phase in [0, 0.25, 0.5, 0.75]:
  437. tot_cycles = cycle_per_deg * 120# 1920/px_per_deg#tot_pixels/px_per_deg
  438. pix_per_cycle = width/tot_cycles
  439. grat_imgs.append( get_spatial_grating(height=height, aspect_ratio=aspect, ori=ori, pix_per_cycle=pix_per_cycle, phase=phase) )
  442. # %%
  443. f, axes = plt.subplots(11,11,figsize=(7.5*aspect,7.5))
  444. shuffled_ix = np.random.choice(range(len(grat_imgs)), len(grat_imgs), False)
  445. for i in range(len(grat_imgs)):
  446. ax = axes.ravel()[i]
  447. image = grat_imgs[shuffled_ix[i]]
  448. ax.imshow(image, cmap=plt.cm.gray)
  449. ax.set(xticks=[],yticks=[])
  450. for i in range(len(grat_imgs),11**2):
  451. ax = axes.ravel()[i]
  452. ax.axis('off')
  454. plt.subplots_adjust(wspace=0.2, hspace=0.2)
  455. plt.show()
  458. # %% [markdown]
  459. # ### collect trial N spks
  461. # %%
  462. ### COLLECTING SPIKE COUNTS -- ALL PHASES
  464. AREAs = ['VISal','VISrl','VISam','VISpm', 'VISp', 'VISl']
  465. datadir = 'data'
  466. STIM_CLASS = 'static_gratings'
  467. PARAM_NAMES = ['spatial_frequency', 'orientation', 'phase']
  469. prestim_len = 0 #sec
  470. bin_size = .0005 #sec
  473. for session_i in range(brain_observatory_type_sessions.index.values.size):
  474. session_id = brain_observatory_type_sessions.index.values[session_i]
  476. print(f'\n* {session_i}: session_id',session_id,flush=True)
  478. skip = True
  479. for AREA in AREAs:
  480. fname = f's{session_id}_{AREA}_{STIM_CLASS}_allPhases'
  481. if os.path.isfile(f'{datadir}/{fname}_trial_info.pkl'):
  482. print(AREA,'saved previously.')
  483. else:
  484. skip = False
  485. break
  486. if skip:
  487. print('skipping...')
  488. continue
  489. session = cache.get_session_data(session_id)
  491. if STIM_CLASS not in session.stimulus_names:
  492. print(f'{STIM_CLASS} not found.')
  493. continue
  494. for AREA in AREAs:
  495. fname = f's{session_id}_{AREA}_{STIM_CLASS}_allPhases'
  496. if os.path.isfile(f'{datadir}/{fname}_trial_info.pkl'):
  497. print(AREA,'saved previously.')
  498. continue
  499. my_units = session.units[session.units["ecephys_structure_acronym"].values == AREA].index.values
  500. Nunits = len(my_units)
  501. print(AREA,'Nunits',Nunits)
  502. if Nunits == 0:
  503. continue
  505. pref_phases = cache.get_unit_analysis_metrics_for_session(session_id)['pref_phase_sg']
  508. #compute for all cells
  509. all_trains = get_spike_trains(my_units, STIM_CLASS, PARAM_NAMES, prestim_len, bin_size)
  510. if not all_trains:
  511. print('0 trains')
  512. continue
  513. uis = sorted(all_trains)
  514. stims = sorted(all_trains[uis[0]])
  515. stim_ntrials = [len(all_trains[uis[0]][si]) for si in stims]
  516. X = None
  517. prefPhases = []
  518. for ii,ui in enumerate(uis):
  519. print(ii,end=' ',flush=True)
  520. prefPhase = pref_phases.loc[ui]
  521. all_frs = []
  522. for si in stims:
  523. all_frs += list(map(len,all_trains[ui][si]))
  524. if X is None:
  525. X = np.asarray(all_frs)[None,:]
  526. else:
  527. X = np.concatenate([X,np.asarray(all_frs)[None,:]])
  528. prefPhases.append(prefPhase)
  529. print()
  530. cellData = {'uis':uis, 'stims':stims, 'stim_ntrials':stim_ntrials, 'prefPhases':prefPhases}
  531. with open(f'{datadir}/{fname}_trial_info.pkl', 'wb') as f:
  532. pickle.dump(cellData, f)
  533. np.save(f'{datadir}/{fname}_trial_data.npy',X)
  534. print(fname,'saved.')
  538. # %%
  539. ### COLLECTING SPK COUNTS -- PREF PHASE
  541. AREAs = ['VISp','VISl','VISal','VISrl','VISam','VISpm']
  542. datadir = 'data'
  543. STIM_CLASS = 'static_gratings'
  544. PARAM_NAMES = ['spatial_frequency', 'orientation']
  546. prestim_len = 0 #sec
  547. bin_size = .0005 #sec
  549. trial_len_ms = 250
  550. prestim_len_ms = prestim_len * 1000
  551. dirs = [0.0, 30.0, 60.0, 90.0, 120.0, 150.0]
  553. bw = 25
  554. samp_interval = 10
  555. NDIRS = 6
  556. SFs = [0.02, 0.04, 0.08, 0.16, 0.32]
  559. for session_i in range(brain_observatory_type_sessions.index.values.size):
  560. session_id = brain_observatory_type_sessions.index.values[session_i]
  562. print(f'\n* {session_i}: session_id',session_id,flush=True)
  564. skip = True
  565. for AREA in AREAs:
  566. fname = f's{session_id}_{AREA}_{STIM_CLASS}_prefPhase'
  567. if os.path.isfile(f'{datadir}/{fname}_trial_info.pkl'):
  568. print(AREA,'saved previously.')
  569. else:
  570. skip = False
  571. break
  572. if skip:
  573. print('skipping...')
  574. continue
  575. session = cache.get_session_data(session_id)
  577. if STIM_CLASS not in session.stimulus_names:
  578. print(f'{STIM_CLASS} not found.')
  579. continue
  580. for AREA in AREAs:
  581. fname = f's{session_id}_{AREA}_{STIM_CLASS}_prefPhase'
  582. if os.path.isfile(f'{datadir}/{fname}_trial_info.pkl'):
  583. print(AREA,'saved previously.')
  584. continue
  585. my_units = session.units[session.units["ecephys_structure_acronym"].values == AREA].index.values
  586. Nunits = len(my_units)
  587. print(AREA,'Nunits',Nunits)
  588. if Nunits == 0:
  589. continue
  591. #compute for all cells
  592. all_trains = get_spike_trains_pref_phase(my_units, STIM_CLASS, PARAM_NAMES, prestim_len, bin_size)
  593. if not all_trains:
  594. print('0 trains')
  595. continue
  596. uis = sorted(all_trains)
  597. stims = sorted(product(*[stim_params[pname] for pname in PARAM_NAMES]))
  598. X = []
  599. prefPhases = []
  600. stim_ntrials = []
  601. for ii,ui in enumerate(uis):
  602. print(ii,end=' ',flush=True)
  603. prefPhase = list(all_trains[ui].keys())[0][-1]
  604. assert np.all(np.array(list(all_trains[ui].keys()))[:,2] == prefPhase)
  605. stim_ntrials.append([len(all_trains[ui][si+(prefPhase,)]) for si in stims])
  606. all_frs = []
  607. for si in stims:
  608. all_frs += list(map(len,all_trains[ui][si+(prefPhase,)]))
  610. X.append(np.asarray(all_frs))
  611. prefPhases.append(prefPhase)
  613. print()
  614. cellData = {'uis':uis, 'prefPhases':prefPhases, 'stims':stims, 'stim_ntrials':stim_ntrials}
  615. with open(f'{datadir}/{fname}_trial_info.pkl', 'wb') as f:
  616. pickle.dump(cellData, f)
  617. np.save(f'{datadir}/{fname}_trial_data.npy',X)
  618. print(fname,'saved.')
  622. # %% [markdown]
  623. # ### collect PSTHs
  625. # %%
  626. AREAs = ['VISl','VISal','VISrl','VISam','VISpm','VISp']
  627. datadir = 'data'
  628. STIM_CLASS = 'static_gratings'
  629. PARAM_NAMES = ['spatial_frequency', 'orientation']
  631. prestim_len = 0.05 #sec
  632. bin_size = .0005 #sec
  634. trial_len_ms = 300
  635. prestim_len_ms = prestim_len * 1000
  636. dirs = [0.0, 30.0, 60.0, 90.0, 120.0, 150.0]
  638. bw = 10
  639. samp_interval = 5
  640. NDIRS = 6
  641. SFs = [0.02, 0.04, 0.08, 0.16, 0.32]
  645. for session_i in range(brain_observatory_type_sessions.index.values.size):
  646. session_id = brain_observatory_type_sessions.index.values[session_i]
  648. print(f'\n* {session_i}: session_id',session_id)
  649. skip = True
  650. for AREA in AREAs:
  651. fname = f's{session_id}_{AREA}_{STIM_CLASS}_bw{bw}'
  652. if os.path.isfile(f'{datadir}/{fname}.pkl'):
  653. print(AREA,'saved previously.')
  654. else:
  655. skip = False
  656. break
  657. if skip:
  658. print('skipping...')
  659. continue
  660. session = cache.get_session_data(session_id)
  662. if STIM_CLASS not in session.stimulus_names:
  663. print(f'{STIM_CLASS} not found.')
  664. continue
  665. for AREA in AREAs:
  666. fname = f's{session_id}_{AREA}_{STIM_CLASS}_bw{bw}'
  667. if os.path.isfile(f'{datadir}/{fname}.pkl'):
  668. print(AREA,'saved previously.')
  669. continue
  670. my_units = session.units[session.units["ecephys_structure_acronym"].values == AREA].index.values
  671. Nunits = len(my_units)
  672. print(AREA,'Nunits',Nunits)
  673. if Nunits == 0:
  674. continue
  675. #compute for all cells
  676. all_trains = get_spike_trains_pref_phase(my_units, STIM_CLASS, PARAM_NAMES, prestim_len, bin_size)
  677. if not all_trains:
  678. print('0 trains')
  679. continue
  680. allData = {}
  681. for uid in my_units:
  682. allData[uid] = {}
  684. for tfi,tf in enumerate(SFs):
  685. allData[uid][tf] = {}
  687. full_traindict, trials_traindict, ISI_Nspks = get_train_dicts(all_trains[uid], tf, dirs, trial_len_ms, prestim_len_ms)
  688. all_ISJs, x_ts = getResponseCurve(full_traindict, trial_len_ms+prestim_len_ms, bw, samp_interval)
  689. allData[uid][tf]['psts'] = all_ISJs
  690. if prestim_len_ms > 0:
  691. allData[uid][tf]['stats'] = computeResponseStats(trials_traindict, ISI_Nspks, prestim_len_ms, trial_len_ms)
  692. allData[uid][tf]['n_trials'] = [full_traindict[d] for d in dirs]
  694. with open(f'{datadir}/{fname}.pkl', 'wb') as f:
  695. pickle.dump(allData, f)
  696. print(fname,'saved.')

read_spike_data-static-gratings.ipynb at commit ddccf17, under BSD-2-Clause · at the source

Overview

Authors: Luciano Dyballa1, Greg D. Field2, Michael P. Stryker3,4, Steven W. Zucker5,6
ORCID iDs: Luciano Dyballa, Michael P. Stryker, Steven W. Zucker
  1. School of Science and Technology, IE University, Madrid, Spain
  2. Jules Stein Eye Institute, Department of Ophthalmology, David Geffen School of Medicine, University of California, Los Angeles, California, United States of America
  3. Department of Physiology, University of California, San Francisco, California, United States of America
  4. Kavli Institute for Fundamental Neuroscience, University of California, San Francisco, California, United States of America
  5. Department of Computer Science, Yale University, New Haven, Connecticut, United States of America
  6. Department of Biomedical Engineering, Yale University, New Haven, Connecticut, United States of America
Institutions: IE University (Spain); Yale University (United States)
Journal: PloS one, volume 21, issue 9, article e0356243
Dates: received 4 February 2026; accepted 1 August 2026; published online 17 September 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1371/journal.pone.0356243 · PMID 42752452 · PMCID PMC13585280 · OpenAlex W7213466761
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: mouse (organism)
Methods: Spectral & time-frequency, Statistics, Smoothing, state filtering, decompositions, Preprocessing, Connectivity, Machine learning, Single-unit activity, calcium imaging
MeSH: Primary Visual Cortex*, Visual Cortex*, Visual Perception*, Animals, Mice, Neurons, Photic Stimulation (* major topic)
Journal subjects: Biology and Life Sciences, Cell Biology, Cellular Types, Animal Cells, Neurons, Neuroscience, Cellular Neuroscience, Physical Sciences, Mathematics, Topology, Manifolds, Cognitive Science, Cognitive Psychology, Perception, Sensory Perception, Vision, Psychology, Social Sciences, Neuronal Tuning, Geometry, Non-Euclidean geometry, Discrete Mathematics, Combinatorics, Permutation, Anatomy, Brain, Visual Cortex, Medicine and Health Sciences
Topic: Visual perception and processing mechanisms (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: National Institutes of Health (NIH) (EY031059); National Science Foundation (NSF) (1822598); European Commission's Marie Skłodowska-Curie Action (101207931); Jules Stein Eye Institute from Research to Prevent Blindness (NIH) (EY000331)
Citations: cited by 1 paper (Europe PMC); 76 references in the paper

Abstract

A challenge in sensory neuroscience is understanding how populations of neurons operate in concert to represent diverse stimuli. To meet this challenge, we have created “encoding manifolds” that reveal the overall responses of brain areas to diverse stimuli and organize individual neurons in stimulus-response coordinates according to their selectivity and response dynamics. Here we use encoding manifolds to compare the population-level encoding of primary visual cortex (VISp) with that of five higher visual areas (VISam, VISal, VISpm, VISlm, and VISrl), using data from the Allen Institute Visual Coding–Neuropixels dataset from the mouse. We show that the topology of the encoding manifold for VISp and for higher visual areas is continuous, with smooth coordinates along which stimulus selectivity and response dynamics are organized with layer and cell-type specificity. Surprisingly, the manifolds revealed novel relationships between how natural scenes are encoded relative to static gratings—a relationship conserved across visual areas. Namely, neurons preferring natural scenes preferred either low or high spatial frequency gratings, but not intermediate ones. Analyzing responses by cortical layer reveals a preference for gratings concentrated in layer 6, whereas preferences for natural scenes tended to be higher in layers 2/3 and 4. The results demonstrate how machine learning approaches can be used to organize and visualize the structure of sensory coding, thereby revealing novel relationships within and across brain areas and sensory stimuli.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 18 matches between paragraphs and lines of code.

dyballa/NeuralEncodingManifolds

License: BSD-2-Clause
State: the link answers, verified on 26 September 2026
Evidence: files inventoried
Commit: ddccf17c8dc1e3d129624c05285dd5ec24a97e41, 24 October 2025
Languages: Jupyter (17), Python (5), MATLAB (4)
Size: 542 files, 26 scripts
Software Heritage: not archived
Found in: the text, “Neural encoding manifolds.”
Holds: README, license file, 17 notebooks
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Matplotlib (22 files), NumPy (22 files), SciPy (13 files), Pillow (5 files), AllenSDK (4 files), pandas (4 files), scikit-learn (4 files), Keras (3 files), TensorFlow (3 files), Plotly (1 file)
Availability: 1 check, the latest on 26 September 2026: the link answers
  • 26 September 2026: the link answers
28 files
At the source: github.com/dyballa/NeuralEncodingManifolds

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 26 scripts, each with its path and the digest of its content;
  • 18 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Data Availability

The data underlying the results presented in the study are available from the Allen Institute Neuropixels dataset: https://portal.brain-map.org/circuits-behavior/visual-coding-neuropixels.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 4 authors, 7 MeSH terms, 4 funders, 65 references.

Cite

This paper

Dyballa, L., Field, G. D., Stryker, M. P., & Zucker, S. W. (2026). Functional organization and natural scene responses across mouse visual cortical areas revealed with encoding manifolds. PloS one, 21(9), e0356243. https://doi.org/10.1371/journal.pone.0356243

BibTeX

@article{dyballa2026functional,
author = {Dyballa, Luciano and Field, Greg D. and Stryker, Michael P. and Zucker, Steven W.},
title = {{Functional organization and natural scene responses across mouse visual cortical areas revealed with encoding manifolds}},
journal = {PloS one},
year = {2026},
month = sep,
volume = {21},
number = {9},
pages = {e0356243},
publisher = {PLOS},
issn = {1932-6203},
doi = {10.1371/journal.pone.0356243},
url = {https://doi.org/10.1371/journal.pone.0356243},
pmid = {42752452},
pmcid = {PMC13585280}
}

RIS

TY - JOUR
AU - Dyballa, Luciano
AU - Field, Greg D.
AU - Stryker, Michael P.
AU - Zucker, Steven W.
TI - Functional organization and natural scene responses across mouse visual cortical areas revealed with encoding manifolds
T2 - PloS one
J2 - PLoS One
PY - 2026
DA - 2026/09/17
VL - 21
IS - 9
SP - e0356243
SN - 1932-6203
PB - PLOS
DO - 10.1371/journal.pone.0356243
UR - https://doi.org/10.1371/journal.pone.0356243
LA - en
ER -

CSL-JSON

{
"id": "10.1371/journal.pone.0356243",
"type": "article-journal",
"title": "Functional organization and natural scene responses across mouse visual cortical areas revealed with encoding manifolds",
"container-title": "PloS one",
"author": [
{
"family": "Dyballa",
"given": "Luciano"
},
{
"family": "Field",
"given": "Greg D."
},
{
"family": "Stryker",
"given": "Michael P."
},
{
"family": "Zucker",
"given": "Steven W."
}
],
"container-title-short": "PLoS One",
"volume": "21",
"issue": "9",
"page": "e0356243",
"DOI": "10.1371/journal.pone.0356243",
"PMID": "42752452",
"PMCID": "PMC13585280",
"ISSN": "1932-6203",
"publisher": "PLOS",
"URL": "https://doi.org/10.1371/journal.pone.0356243",
"language": "en",
"issued": {
"date-parts": [
[
2026,
9,
17
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1126/sciadv.aed6417 [code]
Intrinsic timing, not temporal prediction, underlies ramping dynamics in visual and parietal cortex during passive behavior.
Journal: Science advances
In common: Plotly, scikit-learn, pandas, 3 other tools, mouse, 4 references
[2] doi:10.1186/s13059-026-04177-w [code]
Genomic sequence evolution underlying human neocortical interareal diversification.
Journal: Genome biology
In common: AllenSDK, TensorFlow, Plotly, 6 other tools, mouse
[3] doi:10.1371/journal.pcbi.1013138 [code]
Hierarchical recurrent temporal prediction as a model of the mammalian dorsal visual pathway.
Journal: PLoS computational biology
In common: AllenSDK, scikit-learn, pandas, 3 other tools, 3 references
[4] doi:10.1038/s41467-026-76939-w [code]
HIPPIE: a generative model for electrophysiological analysis across species, technologies, and modalities.
Journal: Nature communications
In common: AllenSDK, Pillow, scikit-learn, 4 other tools, mouse, 2 references
[5] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: Keras, TensorFlow, Plotly, 6 other tools, mouse
[6] doi:10.1126/sciadv.aed3650 [code]
Truthful visualizations for mass spectrometry imaging enable high-spatial-resolution interactive &lt;i&gt;m/z&lt;/i&gt; mapping and exploration.
Journal: Science advances
In common: Keras, TensorFlow, Pillow, 5 other tools, mouse, 1 reference
[7] doi:10.1523/eneuro.0023-26.2026 [code]
Real-Time Segmentation and Classification of Birdsong Syllables for Learning Experiments.
Journal: eNeuro
In common: Keras, TensorFlow, Plotly, 6 other tools
[8] doi:10.1186/s12880-026-02481-2 [code]
Deep learning-based neuroanatomical profiling reveals population-specific brain changes in multiple sclerosis: a large-scale Middle Eastern study.
Journal: BMC medical imaging
In common: Keras, TensorFlow, Plotly, 6 other tools
[9] doi:10.3389/fnsys.2026.1822122 [code]
Convergence-divergence circuits for multimodal integration of innate and learned opponent valences.
Journal: Frontiers in systems neuroscience
In common: Keras, TensorFlow, Plotly, 6 other tools
[10] doi:10.1038/s41597-025-05174-7 [code]
A large-scale MEG and EEG dataset for object recognition in naturalistic scenes
Journal: —
In common: Keras, TensorFlow, Plotly, 6 other tools

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.