OSCR

Hierarchical recurrent temporal prediction as a model of the mammalian dorsal visual pathway.

Code ↔ Paper

8 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 8 matches
  1. [1] § Materials and methods › Data analysis ↔ virtual_physiology/VirtualNetworkPhysiology.py, lines 490–585 · score 0.95 · partial correlation, plaid angle, plaid component, units responded, individual components, Plaid responses
  2. [2] § Materials and methods › Data analysis ↔ analysis/surround_suppression.ipynb, lines 506–643 · score 0.76 · Surround suppression, preferred orientation, tuning curve, grating stimulus, temporal frequency, mask
  3. [3] § Results › Model units mirror the visual system’s hierarchy of 2D motion sensitivity ↔ virtual_physiology/VirtualNetworkPhysiology.py, lines 490–585 · score 0.70 · plaid angles, plaid response, tuning curves, grating stimuli, Plaid pattern, Fisher
  4. [4] § Materials and methods › Model-brain alignment ↔ neural_fitting/get_neural_data.py, lines 9–29 · score 0.58 · signal power, noise power, Neural
  5. [5] § Results › Model units mirror the visual system’s hierarchy of 2D motion sensitivity ↔ analysis/plaid_motion.ipynb, lines 68–98 · score 0.57 · plaid motion, selective units, plaid pattern, DSI, tuned
  6. [6] § Materials and methods › The hierarchical recurrent temporal prediction model ↔ models/network_hierarchical_recurrent.py, lines 149–163 · score 0.56 · enforce Dale, Law, excitatory, inhibitory, weight, recurrent
  7. [7] § Results › The hierarchical recurrent temporal prediction model matches cortical stimulus representations better than other models ↔ analysis/surround_suppression.ipynb, lines 506–643 · score 0.54 · surround suppression, tuning curve, full models, hierarchical recurrent, feedback, stimulus
  8. [8] § Results › The hierarchical recurrent temporal prediction model matches cortical stimulus representations better than other models ↔ analysis/plaid_motion.ipynb, lines 308–419 · score 0.54 · KS distance, plaid pattern, iv, iii, G3, G2

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

Python · 709 lines · 30 KB · no license · 2 matches

  1. import pickle, math, os
  3. import numpy as np
  4. import torch
  6. import scipy
  7. import scipy.stats
  8. from scipy import ndimage
  9. import scipy.optimize as opt
  10. from scipy import signal
  11. import scipy.fft as fft
  13. class VirtualPhysiology:
  14. # model trained model
  15. # hyperparameters model hyperparameters
  16. # hidden_unit_range range object of hidden units for analysis
  17. # device device for tensors (cpu or cuda)
  18. def __init__ (self, model, hyperparameters, frame_shape, hidden_units, device):
  19. self.data = []
  21. self.model = model
  22. self.model.eval()
  24. self.hyperparameters = hyperparameters
  25. self.warmup = hyperparameters["warmup"]
  26. self.frame_size = hyperparameters["frame_size"]
  27. self.t_steps = 50
  29. self.frame_shape = frame_shape
  30. self.hidden_units = hidden_units
  31. self.device = device
  33. self.data = []
  34. for group in self.hidden_units:
  35. self.data.append([])
  37. self.osi_thresh = 0.4
  38. self.dsi_thresh = 0.3
  39. self.mean_response_offset = 5
  41. min_n, max_n = min(frame_shape), max(frame_shape)
  42. self.spatial_frequencies = [i/max_n for i in range(1, min_n//2+1)] # Cycles / pixel
  43. self.orientations = np.arange(0, 360, 5) # Degrees
  44. self.temporal_frequencies = np.linspace(1/self.t_steps, 1/4, self.t_steps//4) # Cycles / frame
  46. @classmethod
  47. def load (cls, data_path, model, hyperparameters, frame_shape, hidden_units, device):
  48. virtual_physiology = cls(
  49. model=model,
  50. hyperparameters=hyperparameters,
  51. frame_shape=frame_shape,
  52. hidden_units=hidden_units,
  53. device=device
  54. )
  56. with open(data_path, 'rb') as handler:
  57. virtual_physiology.data = pickle.load(handler)
  59. return virtual_physiology
  61. # relative_path = use file name as part of data_path relative to model dir
  62. def save (self, data_path):
  63. with open(data_path, 'wb') as p:
  64. pickle.dump(self.data, p, protocol=4)
  66. with open(data_path + '.params', 'wb') as p:
  67. params = {
  68. "t_steps": self.t_steps,
  69. "spatial_frequencies": self.spatial_frequencies,
  70. "orientations": self.orientations,
  71. "temporal_frequencies": self.temporal_frequencies
  72. }
  73. pickle.dump(params, p, protocol=4)
  76. def get_unit_data (self, unit_idx):
  77. for group in self.data:
  78. for unit_data in group:
  79. if unit_data["hidden_unit_index"] == unit_idx:
  80. return unit_data
  81. return False
  84. def get_group_from_unit_idx (self, unit_idx):
  85. for group_idx in range(len(self.hidden_units)):
  86. if unit_idx < np.sum(self.hidden_units[:group_idx+1]):
  87. return group_idx
  88. return False
  90. def get_response_weighted_average (self, n_rand_stimuli=100):
  91. # Pre-allocate lists of lists to hold gaussian noise stimuli and associated unit activity
  92. stimuli = []
  93. unit_activity = []
  94. for _ in range (np.sum(self.hidden_units)):
  95. stimuli.append([])
  96. unit_activity.append([])
  98. noise_shape = (n_rand_stimuli, self.warmup+self.t_steps, self.frame_size)
  99. noise = np.random.normal(loc=0, scale=1, size=noise_shape)
  100. noise = torch.Tensor(noise).to(self.device)
  102. with torch.no_grad():
  103. _, hidden_state = self.model(noise)
  106. response = hidden_state.detach().numpy().reshape(n_rand_stimuli*(self.warmup+self.t_steps), -1)
  107. noise = noise.detach().numpy().reshape(-1, self.frame_size)
  109. rwa_arr = []
  110. for unit_idx, unit_responses in enumerate(response.T):
  111. group_idx = self.get_group_from_unit_idx(unit_idx)
  113. if unit_idx % 100 == 0:
  114. print('Processing RWA for unit', unit_idx)
  116. if np.sum(unit_responses):
  117. if group_idx == 0:
  118. rwa = np.average(noise, axis=0, weights=unit_responses)
  119. else:
  120. rwa = np.average(noise[:-group_idx], axis=0, weights=unit_responses[group_idx:])
  122. self.data[group_idx].append({
  123. "hidden_unit_index": unit_idx,
  124. "response_weighted_average": rwa
  125. })
  127. print('Finished averaging stimuli')
  129. return self
  131. # sf = cycles per pixel
  132. # tf = cycles per second
  133. # speed = tf/sf = pixels per second
  134. def get_grating_stimuli(self, spatial_frequency, orientation, temporal_frequency, grating_amplitude, frames):
  135. y_size, x_size = self.frame_shape
  137. theta = (orientation-90) * np.pi/180
  138. x, y = np.meshgrid(np.arange(0, x_size), np.arange(0, y_size))
  139. x_theta = x * np.cos(theta) + y * np.sin(theta)
  141. phase_shift = 2*np.pi*temporal_frequency
  142. phases = np.arange(frames)*phase_shift
  144. grating_frames = []
  145. for phase in phases:
  146. grating_frames.append( grating_amplitude * np.sin(2*spatial_frequency*np.pi*x_theta - phase) )
  149. gratings = np.array(grating_frames).reshape(1, frames, y_size*x_size)
  150. gratings = (gratings-np.mean(gratings))/np.std(gratings) ##
  151. gratings = torch.Tensor(gratings).to(self.device)
  153. return gratings
  155. def get_grating_responses (self):
  156. # Add array to data dictionary structures containing
  157. # complete response (for each grating phase) and mean response (averaged across phases)
  158. # for each spatial frequency/orientation/tf combination
  159. for group_data in self.data:
  160. for unit_data in group_data:
  161. unit_data["grating_responses"] = np.zeros((
  162. len(self.spatial_frequencies),
  163. len(self.orientations),
  164. len(self.temporal_frequencies),
  165. self.t_steps
  166. ))
  167. unit_data["mean_grating_responses"] = np.zeros((
  168. len(self.spatial_frequencies),
  169. len(self.orientations),
  170. len(self.temporal_frequencies)
  171. ))
  173. # Keep track of progress for display purposes
  174. param_count = 0
  175. try:
  176. total_params = self.data[0][0]["mean_grating_responses"].size
  177. except:
  178. total_params = self.data[1][0]["mean_grating_responses"].size
  180. # Loop through each parameter combination for each unit
  181. for sf_idx, sf in enumerate(self.spatial_frequencies):
  182. for ori_idx, ori in enumerate(self.orientations):
  183. for tf_idx, tf in enumerate(self.temporal_frequencies):
  185. # Generate gratings for particular param combination
  186. gratings = self.get_grating_stimuli(sf, ori, tf, grating_amplitude=1, frames=self.warmup+self.t_steps)
  188. # Feedforward pass through network
  189. with torch.no_grad():
  190. _, hidden_state = self.model(gratings)
  192. # Loop through unit responses at each time step
  193. for group_data in self.data:
  194. for unit_data in group_data:
  195. unit_idx = unit_data["hidden_unit_index"]
  197. # Discard warm up period of network's response to gratings
  198. unit_responses = hidden_state[0, self.warmup:, unit_idx].cpu().numpy()
  200. unit_data["grating_responses"][sf_idx, ori_idx, tf_idx] = unit_responses
  201. unit_data["mean_grating_responses"][sf_idx, ori_idx, tf_idx] = np.mean(unit_responses)
  203. if param_count % 100 == 99:
  204. print("Finished param combination {}/{}".format(param_count+1, total_params))
  205. param_count += 1
  207. print("Finished tuning curve")
  209. return self
  211. # Takes orientation tuning curve at max tf and sf
  212. # Returns direction selectivity (DSI)
  213. def get_DSI (self, tuning_curve):
  214. orient_pref_idx = np.where(tuning_curve == np.max(tuning_curve))[0][0]
  215. orient_pref = self.orientations[orient_pref_idx]
  216. orient_pref_resp = tuning_curve[orient_pref_idx]
  218. orient_opp = (orient_pref + 180) % 360
  219. orient_opp_idx = np.where(self.orientations == orient_opp)[0][0]
  220. orient_opp_resp = tuning_curve[orient_opp_idx]
  222. DSI = (orient_pref_resp - orient_opp_resp) / (orient_pref_resp + orient_opp_resp)
  223. #DSI = 1 - (orient_opp_resp/orient_pref_resp)
  225. return DSI
  228. # Takes orientation tuning curve at max tf and sf
  229. # Returns orientation selectivity (OSI)
  230. def get_OSI (self, tuning_curve):
  231. orient_pref_idx = np.where(tuning_curve == np.max(tuning_curve))[0][0]
  232. orient_pref = self.orientations[orient_pref_idx]
  233. orient_pref_resp = tuning_curve[orient_pref_idx]
  235. orient_orth1 = (orient_pref + 90) % 360
  236. orient_orth1_idx = np.where(self.orientations == orient_orth1)[0][0]
  237. orient_orth2 = (orient_pref - 90) % 360
  238. orient_orth2_idx = np.where(self.orientations == orient_orth2)[0][0]
  239. orient_orth_resp = (tuning_curve[orient_orth1_idx]+tuning_curve[orient_orth2_idx]) / 2
  241. OSI = (orient_pref_resp - orient_orth_resp) / (orient_pref_resp + orient_orth_resp)
  243. return OSI
  245. def get_orientation_tuning_curve(self, unit_data):
  246. mean_grating_responses = unit_data["mean_grating_responses"]
  248. # Get indices of max grating response (mean across time steps)
  249. max_mean_grating_response = np.max(mean_grating_responses)
  250. unit_data["max_mean_grating_response"] = max_mean_grating_response
  251. sf_idx, ori_idx, tf_idx = [idx[0] for idx in np.where(mean_grating_responses == max_mean_grating_response)]
  253. orientation_tuning_curve = mean_grating_responses[sf_idx, :, tf_idx]
  254. return orientation_tuning_curve
  257. # Gets OSI, DSI, CV and modulation ratio for each unit
  258. def get_grating_responses_parameters (self):
  259. for group_idx, group_data in enumerate(self.data):
  260. for unit_i, unit_data in enumerate(group_data):
  261. grating_responses = unit_data["grating_responses"]
  262. mean_grating_responses = unit_data["mean_grating_responses"]
  264. # Get indices of max grating response (mean across time steps)
  265. max_mean_grating_response = np.max(mean_grating_responses)
  266. unit_data["max_mean_grating_response"] = max_mean_grating_response
  267. sf_idx, ori_idx, tf_idx = [idx[0] for idx in np.where(mean_grating_responses == max_mean_grating_response)]
  269. # Convert these indices into the underlying parameters
  270. max_sf = unit_data["preferred_sf"] = self.spatial_frequencies[sf_idx]
  271. max_ori = unit_data["preferred_orientation"] = self.orientations[ori_idx]
  272. max_tf = unit_data["preferred_tf"] = self.temporal_frequencies[tf_idx]
  274. # Response to moving grating for parameters that give maximum response
  275. optimum_grating_response = grating_responses[sf_idx, ori_idx, tf_idx]
  276. unit_data["optimum_grating_response"] = optimum_grating_response
  278. # Get CV and DSI measures from curve
  279. # given max spatial frequency and temporal frequency parameters
  280. orientation_curve = mean_grating_responses[sf_idx, :, tf_idx]
  281. unit_data["OSI"] = self.get_OSI(orientation_curve)
  282. unit_data["DSI"] = self.get_DSI(orientation_curve)
  284. if unit_i % 50 == 49:
  285. print("Finished unit {} / {}, group {} / {}".format(
  286. unit_i+1, len(group_data), group_idx+1, len(self.data)
  287. ))
  289. return self
  291. # Filter unresponsive units and units where curve fitting failed
  292. def filter_unit_data (self, group_idx):
  293. group_data = self.data[group_idx]
  295. # Reject low mean response units (< 1% of mean, max mean response)
  296. all_max_mean = []
  297. for group in self.data:
  298. for u in group:
  299. all_max_mean.append(u["max_mean_grating_response"])
  301. mean_responses = [u["max_mean_grating_response"] for u in group_data]
  303. response_threshold = 0.01 * np.mean(all_max_mean)
  304. n_filtered = len(np.where(mean_responses < response_threshold)[0])
  305. print("{} / {} units below response threshold".format(n_filtered, len(group_data)))
  307. # Reject units where curve fitting failed (modulation_ratio set as False)
  308. n_filtered = len([u for u in group_data if not u["modulation_ratio"]])
  309. print("{} / {} units failed to fit curve for modulation ratio estimate".format(n_filtered, len(group_data)))
  311. # Now actually filter those units out
  312. return [u for u in group_data if u["modulation_ratio"] and u["max_mean_grating_response"] >= response_threshold]
  314. # Filter unresponsive units (in place method)
  315. def filter_nonresponding_units (self):
  316. for group_idx, group in enumerate(self.data):
  317. filtered = []
  318. for u in group:
  319. max_mean = np.max(u["mean_grating_responses"])
  320. if max_mean > 0:
  321. filtered.append(u)
  323. self.data[group_idx] = filtered
  325. original_n = self.hidden_units[group_idx]
  326. filtered_n = len(self.data[group_idx])
  327. print(f"{filtered_n} / {original_n} units kept after filtering non-responsive units")
  329. return self
  332. def get_moving_bar_stimuli (self, direction, x, y, bar_amplitude=1, bar_size=5, frames_len = 20):
  333. # Include warmup period for network in stimulus
  334. # Where first n frames will be gray (all 0's)
  335. total_frames = self.warmup + frames_len
  337. stimuli = np.zeros((total_frames, self.frame_shape[0], self.frame_shape[1]))
  339. for i in range(frames_len):
  340. bar_position = 0
  341. square = np.ones((bar_size, bar_size))*bar_amplitude
  343. if direction == 0 or direction == 270:
  344. bar_position = (bar_size - i) % bar_size
  345. else:
  346. bar_position = i%(bar_size)
  348. if direction == 0 or direction == 180:
  349. square[bar_position, :] = -bar_amplitude
  350. else:
  351. square[:, bar_position] = -bar_amplitude
  353. stimulus = np.zeros((self.frame_shape[0], self.frame_shape[1]))
  354. stimulus[y:y+bar_size, x:x+bar_size] = square
  356. stimuli[self.warmup+i, :, :] = stimulus
  358. # Reshape frame into flat array and convert into Tensor object
  359. stimuli = stimuli.reshape(total_frames, self.frame_size)
  360. stimuli = torch.Tensor(stimuli).unsqueeze(0).to(self.device)
  362. return stimuli
  364. def fit_sine (self, x, y):
  365. # Fit to sine
  366. def func(x, a, b, c, d):
  367. return a*np.sin(b*x + c) + d
  369. def get_mse_loss (y, y_est):
  370. return np.sum((y-y_est)**2)/ len(y)
  372. # Get r_squared from https://stackoverflow.com/a/37899817
  373. def get_rsq (y, y_est):
  374. residuals = y - y_est
  375. ss_res = np.sum(residuals**2)
  376. ss_tot = np.sum((y-np.mean(y))**2)
  377. r_squared = 1 - (ss_res / ss_tot)
  379. return r_squared
  382. best_params = []
  384. for iteration in range(5):
  385. scale = 1-0.2*iteration
  386. n_random_guesses = 10000 if iteration == 0 else 1000
  388. params = []
  389. loss_list = []
  391. for _ in range(n_random_guesses):
  392. if iteration == 0:
  393. rand_params = [
  394. np.random.uniform(low=np.mean(y)-2.5, high=np.mean(y)+2.5),
  395. np.random.uniform(low=0, high=10),
  396. np.random.uniform(low=0, high=len(y)),
  397. np.random.uniform(low=np.min(y)-2.5, high=np.max(y)+2.5)
  398. ]
  399. else:
  400. prev_best_params = best_params[-1]
  401. rand_params_ = [
  402. np.random.uniform(low=-2*scale, high=2*scale),
  403. np.random.uniform(low=-1*scale, high=1*scale),
  404. np.random.uniform(low=-2*scale, high=2*scale),
  405. np.random.uniform(low=-2*scale, high=2*scale)
  406. ]
  407. rand_params = [p+prev_best_params[idx] for idx, p in enumerate(rand_params_)]
  409. # Get the estimated curve based on fitted parameters
  410. y_est = func(x, *rand_params)
  411. loss = get_mse_loss(y, y_est) #get_rsq(y, y_est)
  413. params.append(rand_params)
  414. loss_list.append(loss)
  416. # Get the index of the lowest RSQ, use this to find the
  417. # corresponding parameters used
  418. best_params.append(params[np.argmin(loss_list)])
  420. final_params = best_params[-1]
  421. final_y_est = func(x, *final_params)
  422. final_loss = get_mse_loss(y, final_y_est)
  423. final_rsq = get_rsq(y, final_y_est)
  425. return final_params, final_y_est, final_rsq, final_loss
  428. # Takes list of response as well as a start and
  429. # end offset for where curve fitting should occur
  430. # Returns modulation ratio, estimated curve and RSQ of curve fit
  431. def get_modulation_ratio (self, activity, start_offset, end_offset):
  432. x = np.arange(start_offset, end_offset)
  433. y = activity[start_offset:end_offset]
  435. final_params, final_y_est, final_rsq, final_loss = self.fit_sine (x, y)
  436. # Try one more time if it fails
  437. if final_loss < 0.05 or final_rsq > 0.5:
  438. final_params, final_y_est, final_rsq, final_loss = self.fit_sine (x, y)
  440. # Average unit activity
  441. f0 = np.mean(activity[self.warmup:])
  442. # Absolute of the amplitude of the fitted sine
  443. f1 = (abs(final_params[0]))
  444. mod_ratio = f1/f0
  446. # Reject f values for those units with poor sine fits
  447. #if (final_loss < 0.05 or final_rsq > 0.5) and f0 != 0:
  448. cc = scipy.stats.pearsonr(final_y_est, y)[0]
  449. if cc>0.9 and f0 != 0:
  450. return {
  451. 'modulation_ratio' : mod_ratio,
  452. 'modulation_ratio_all' : mod_ratio,
  453. 'modulation_ratio_cc' : cc,
  454. 'modulation_ratio_y' : final_y_est,
  455. 'modulation_ratio_rsq' : final_rsq,
  456. 'modulation_ratio_loss' : final_loss,
  457. 'modulation_ratio_params' : final_params
  458. }
  459. else:
  460. return {
  461. 'modulation_ratio' : False,
  462. 'modulation_ratio_all' : mod_ratio,
  463. 'modulation_ratio_cc' : cc,
  464. 'modulation_ratio_y' : final_y_est,
  465. 'modulation_ratio_y_true' : final_y_est,
  466. 'modulation_ratio_rsq' : final_rsq,
  467. 'modulation_ratio_loss' : final_loss,
  468. 'modulation_ratio_params' : final_params
  469. }
  471. def get_modulation_ratio_all (self):
  472. for g_idx, g in enumerate(self.data):
  473. for unit_idx, unit_data in enumerate(g):
  474. if unit_idx % 10 == 0:
  475. print(f'Starting group {g_idx}, unit {unit_idx}')
  477. modulation_data = self.get_modulation_ratio(
  478. unit_data['optimum_grating_response'], self.warmup, self.t_steps
  479. )
  481. unit_data["modulation_ratio"] = modulation_data["modulation_ratio"]
  482. unit_data["modulation_ratio_cc"] = modulation_data["modulation_ratio_cc"]
  483. unit_data["modulation_ratio_all"] = modulation_data["modulation_ratio_all"]
  484. unit_data["modulation_ratio_y"] = modulation_data["modulation_ratio_y"]
  485. unit_data["modulation_ratio_rsq"] = modulation_data["modulation_ratio_rsq"]
  486. unit_data["modulation_ratio_loss"] = modulation_data["modulation_ratio_loss"]
  487. unit_data["modulation_ratio_params"] = modulation_data["modulation_ratio_params"]
  490. # https://www.nature.com/articles/nn1786#Sec11
  491. # plaid_angle = spread of plaid components
  492. def get_plaid_pattern_index (self, grating_amplitude=1):
  493. # Test if true response is more similar to pattern or component predictions
  494. # Partial correlation of x and y controlling for z
  495. def partial_correlation (x, y, z):
  496. r_xy, _ = scipy.stats.pearsonr(x, y)
  497. r_xz, _ = scipy.stats.pearsonr(x, z)
  498. r_yz, _ = scipy.stats.pearsonr(y, z)
  500. return (r_xy - r_xz*r_yz) / ((1-r_xz**2)*(1-r_yz**2))**0.5
  502. # Converts from distribution of r values to normal distribution
  503. # (to compare across units)
  504. def fisher_r_to_z (a, df):
  505. return 0.5*np.log((1+a)/(1-a)) / (1/df)**0.5
  507. orientation_step = 6
  508. orientations = self.orientations[::6]
  510. for group_idx, group_data in enumerate(self.data):
  511. for unit_i, unit_data in enumerate(group_data):
  512. max_sf_idx = np.where(np.array(self.spatial_frequencies) == unit_data['preferred_sf'])[0][0]
  513. max_tf_idx = np.where(np.array(self.temporal_frequencies) == unit_data['preferred_tf'])[0][0]
  515. # Response if unit responds to plaid the same as a grating of same angle
  516. pattern_prediction = unit_data['mean_grating_responses'][max_sf_idx, ::orientation_step, max_tf_idx]
  518. z_p_arr = []
  519. z_c_arr = []
  521. for plaid_angle in [60, 90, 120, 150]:
  522. # Response if unit responds to individual components of the plaid
  523. # Shift is the amount to rotate the tuning curve (+/- half the plaid angle)
  524. shift = int( (plaid_angle/2) / (orientations[1]-orientations[0]) )
  525. component_prediction_a = np.roll(pattern_prediction, shift)
  526. component_prediction_b = np.roll(pattern_prediction, -shift)
  527. component_prediction = (component_prediction_a + component_prediction_b) / 2
  529. # Now get 'true' response to the plaid (as a function of orientation)
  530. plaid_response = []
  531. with torch.no_grad():
  532. for ori_idx, ori in enumerate(orientations):
  533. # Plaid at unit's preferred tf and sf
  534. grating_a = self.get_grating_stimuli(
  535. unit_data['preferred_sf'],
  536. ori+plaid_angle/2,
  537. unit_data['preferred_tf'],
  538. grating_amplitude/2,
  539. 50
  540. )
  541. grating_b = self.get_grating_stimuli(
  542. unit_data['preferred_sf'],
  543. ori-plaid_angle/2,
  544. unit_data['preferred_tf'],
  545. grating_amplitude/2,
  546. 50
  547. )
  548. plaid = grating_a + grating_b
  550. # Feed into model
  551. _, hidden_state = self.model(plaid)
  553. # Get mean response for current unit
  554. mean_plaid_response = hidden_state[0, self.warmup:, unit_data['hidden_unit_index']] \
  555. .mean() \
  556. .detach() \
  557. .numpy()
  558. plaid_response.append(mean_plaid_response)
  560. r_p = partial_correlation(plaid_response, pattern_prediction, component_prediction)
  561. r_c = partial_correlation(plaid_response, component_prediction, pattern_prediction)
  563. z_p = fisher_r_to_z(r_p, len(pattern_prediction)-3)
  564. z_c = fisher_r_to_z(r_c, len(component_prediction)-3)
  566. z_p_arr.append(z_p)
  567. z_c_arr.append(z_c)
  569. z_p_mean = np.mean(z_p_arr)
  570. z_c_mean = np.mean(z_c_arr)
  571. pattern_idx = z_p_mean - z_c_mean
  573. # Save all these results
  574. unit_data['plaid_rp'] = r_p
  575. unit_data['plaid_rc'] = r_c
  576. unit_data['plaid_zp'] = z_p_mean
  577. unit_data['plaid_zc'] = z_c_mean
  578. unit_data['plaid_pattern_index'] = pattern_idx
  579. unit_data['plaid_response'] = plaid_response
  581. print("Finished unit {} / {}, group {} / {}".format(
  582. unit_i+1, len(group_data), group_idx+1, len(self.data)
  583. ))
  585. return self
  588. def get_moving_bar_responses (self):
  589. bar_size = 20
  590. # -5 because we are effectively convolving with moving bar "kernel" without padding
  591. convolved_size_r = self.frame_shape[0] - bar_size
  592. convolved_size_pixels_r = np.arange(convolved_size_r)
  593. convolved_size_c = self.frame_shape[1] - bar_size
  594. convolved_size_pixels_c = np.arange(convolved_size_c)
  596. # Pre-allocate array to hold mean response to stimulus at each location/direction
  597. for group_data in self.data:
  598. for unit_data in group_data:
  599. unit_data["moving_bar_responses"] = np.zeros((4, convolved_size_r, convolved_size_c))
  601. # Loop through each orientation and spatial location
  602. for orient_idx, orient in enumerate([0, 90, 180, 270]):
  603. for row_idx, row in enumerate(convolved_size_pixels_r):
  604. for col_idx, col in enumerate(convolved_size_pixels_c):
  606. stimuli = self.get_moving_bar_stimuli(orient, col, row, bar_size=bar_size)
  608. with torch.no_grad():
  609. _, hidden_state = self.model(stimuli)
  611. for group_data in self.data:
  612. for unit_data in group_data:
  613. # Get response at each time step of stimulus
  614. unit_activity = hidden_state[0, :, unit_data["hidden_unit_index"]].cpu().numpy()
  615. # Take baseline response as mean of warmup period
  616. #baseline_response = np.mean(unit_activity[:self.warmup])
  617. # Take overlal mean response as non-warmup period subtracted from baseline
  618. mean_response = np.mean(unit_activity[:]) # - baseline_response)
  620. unit_data["moving_bar_responses"][orient_idx, row_idx, col_idx] = mean_response
  622. print("Got responses for {} degrees, row {}".format(orient, row+1))
  624. # Also compute average response across orientations to produce 'heatmap'
  625. for group_data in self.data:
  626. for unit_data in group_data:
  627. unit_data["mean_moving_bar_responses"] = np.mean(unit_data["moving_bar_responses"], axis=0)
  629. return self
  632. def get_receptive_field_centres (self):
  633. def twoD_Gaussian(xy, amplitude, xo, yo, sigma_x, sigma_y, theta, offset):
  634. xy = (16, 16)
  635. x = np.linspace(0, int(xy[0]-1), int(xy[0]))
  636. y = np.linspace(0, int(xy[1]-1), int(xy[1]))
  637. x,y = np.meshgrid(x, y)
  639. xo = float(xo)
  640. yo = float(yo)
  641. a = (np.cos(theta)**2)/(2*sigma_x**2) + (np.sin(theta)**2)/(2*sigma_y**2)
  642. b = -(np.sin(2*theta))/(4*sigma_x**2) + (np.sin(2*theta))/(4*sigma_y**2)
  643. c = (np.sin(theta)**2)/(2*sigma_x**2) + (np.cos(theta)**2)/(2*sigma_y**2)
  644. g = offset + amplitude*np.exp( - (a*((x-xo)**2) + 2*b*(x-xo)*(y-yo)
  645. + c*((y-yo)**2)))
  647. return g.ravel()
  649. def weighted_centre (resp):
  650. (X,Y) = np.meshgrid(np.arange(0, resp.shape[0]), np.arange(0, resp.shape[1]))
  651. x_coord = (X*resp).sum() / resp.sum().astype("float")
  652. y_coord = (Y*resp).sum() / resp.sum().astype("float")
  653. return (y_coord, x_coord)
  656. for group_idx, group_data in enumerate(self.data):
  657. poor_fit_count = 0
  659. for unit_i, unit_data in enumerate(group_data):
  660. resp = unit_data["mean_moving_bar_responses"]
  662. max_val = np.max(resp)
  663. mean_centre = weighted_centre(resp)
  664. initial_guess = (max_val, mean_centre[1], mean_centre[0], 3, 3, 0, 0)
  666. try:
  667. popt, pcov = opt.curve_fit(
  668. twoD_Gaussian,
  669. resp.shape,
  670. resp.reshape(-1),
  671. p0=initial_guess,
  672. bounds=(
  673. [0, -np.inf, -np.inf, 0, 0, 0, 0],
  674. [np.inf, np.inf, np.inf, np.inf, np.inf, np.pi*2, np.inf]
  675. )
  676. )
  678. resp_gaussian_fit = twoD_Gaussian(resp.shape, *popt)
  680. centre = popt[2], popt[1]
  681. cc = scipy.stats.pearsonr(resp.reshape(-1), resp_gaussian_fit)[0]
  682. except Exception as e:
  683. centre = (np.nan, np.nan)
  685. # 'Convolution' with moving bar means this center needs to be scaled back up to
  686. # full size of frame
  687. scale_factor = self.frame_shape[0]/resp.shape[0]
  688. centre = (centre[0] * scale_factor, centre[1] * scale_factor)
  690. if np.isnan(centre[0]) or np.isnan(centre[1]) or cc < 0.6:
  691. unit_data["receptive_field_centre"] = (-1, -1)
  692. unit_data["receptive_field_gaussian_fit"] = np.nan
  693. unit_data["receptive_field_response"] = resp
  695. poor_fit_count += 1
  696. else:
  697. unit_data["receptive_field_centre"] = centre
  698. unit_data["receptive_field_gaussian_fit"] = resp_gaussian_fit.reshape(16, 16)
  700. if unit_i % 100 == 99:
  701. print("Finished unit {} / {}, group {} / {}".format(
  702. unit_i+1, len(group_data), group_idx+1, len(self.data)
  703. ))
  705. print("{}/{} units were rejected as Gaussian could not be fit to receptive fields in group {}".format(
  706. poor_fit_count, len(group_data), group_idx+1
  707. ))
  709. return self

VirtualNetworkPhysiology.py at commit dc9bb5f, no license · at the source

Overview

Authors: Sebastian Klavinskis-Whiting1, Andrew J. King1, Nicol S. Harper1
ORCID iDs: Andrew J. King, Nicol S. Harper
  1. Department of Physiology, Anatomy and Genetics, University of Oxford, Oxford, United Kingdom
Institutions: University of Oxford (United Kingdom)
Journal: PLoS computational biology, volume 22, issue 5, article e1013138
Dates: received 18 May 2025; accepted 7 May 2026; published online 28 May 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1371/journal.pcbi.1013138 · PMID 42207826 · PMCID PMC13252843 · OpenAlex W4410450465
Open access: gold, a free copy (OpenAlex)
Preprint: doi:10.1101/2025.05.15.654214
Status: code verified
Categories: computational modeling (no new data) (modality), systems (subfield)
Methods: Spectral & time-frequency, Connectivity, Statistics, Machine learning
MeSH: Models, Neurological*, Visual Cortex*, Visual Pathways*, Animals, Computational Biology, Computer Simulation, Neurons, Recurrent Neural Networks (* major topic)
Journal subjects: Research and Analysis Methods, Mathematical and Statistical Techniques, Statistical Methods, Forecasting, Physical Sciences, Mathematics, Statistics, Biology and life sciences, Organisms, Eukaryota, Animals, Vertebrates, Amniotes, Mammals, Primates, Monkeys, Old World monkeys, Macaque, Zoology, Anatomy, Brain, Visual Cortex, Medicine and Health Sciences, Neuroscience, Neuronal Tuning, Cognitive Science, Cognitive Psychology, Perception, Sensory Perception, Vision, Psychology, Social Sciences, Cell Biology, Cellular Types, Animal Cells, Neurons, Cellular Neuroscience, Physiology, Sensory Physiology, Visual System, Sensory Systems, Computer and Information Sciences, Neural Networks
Topic: Morphological variations and asymmetry (Geometry and Topology, Mathematics), according to OpenAlex
Funding: Wellcome Trust (WT108369/Z/2015/Z); Nuffield Department of Clinical Neurosciences, University of Oxford (Not applicable)
Citations: not cited yet (Europe PMC); 71 references in the paper

Abstract

A major goal of neuroscience is to identify general principles that can explain the diverse structures and functions of the brain. The principle of temporal prediction provides one approach, arguing that the sensory brain is optimized to represent stimulus features that efficiently predict the immediate future input. Previous work has demonstrated that feedforward hierarchical temporal prediction models can capture the tuning properties of neurons along the visual pathway, and that recurrent temporal prediction models can explain local functional connectivity within primary visual cortex. However, the visual system is also characterized by extensive inter-areal feedback recurrency, which existing models lack. We aimed to better account for the dynamic features of neurons in the visual cortex by incorporating both local recurrency and inter-areal feedback connectivity into a hierarchical temporal prediction model. The resulting model captured tuning properties along the dorsal visual pathway, including pattern motion selectivity and surround suppression, and the contribution of inter-areal connectivity to these properties. Moreover, compared with several alternative normative models, the hierarchical recurrent temporal prediction model provided the closest fit to these tuning properties and was best able to explain neuronal responses to natural stimuli. Accordingly, temporal prediction accounts well for information processing along the visual pathway.

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 8 matches between paragraphs and lines of code.

sebbkw/hierarchical_temporal_prediction

License: none: the authors keep all their rights
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: dc9bb5f145bd1969d261b1af1531ab6cc6c4d078, 2 May 2026
Languages: Python (26), Jupyter (7), Shell (1)
Size: 81 files, 34 scripts
Software Heritage: not archived
Found in: “Data Availability”
Holds: README, 7 notebooks
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (22 files), PyTorch (21 files), SciPy (12 files), Matplotlib (9 files), Pingouin (6 files), pandas (5 files), OpenCV (3 files), statsmodels (2 files), AllenSDK (1 file), h5py (1 file), scikit-image (1 file), scikit-learn (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
35 files
At the source: github.com/sebbkw/hierarchical_temporal_prediction

The paper's code and data availability statement is in the Data section.

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;
  • 34 scripts, each with its path and the digest of its content;
  • 8 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

All code supporting this study is publicly available at https://github.com/sebbkw/hierarchical_temporal_prediction. The model was trained using a novel dataset drawn from around 2.5 hours of wildlife videos, which were obtained from public sources, and are available at https://osf.io/hf2pj/ (comprising https://doi.org/10.5281/zenodo.20041063; https://doi.org/10.5281/zenodo.20041567; https://doi.org/10.5281/zenodo.20041596).

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, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 3 authors, 8 MeSH terms, 2 funders, 64 references.

Cite

This paper

Klavinskis-Whiting, S., King, A. J., & Harper, N. S. (2026). Hierarchical recurrent temporal prediction as a model of the mammalian dorsal visual pathway. PLoS computational biology, 22(5), e1013138. https://doi.org/10.1371/journal.pcbi.1013138

BibTeX

@article{klavinskiswhiting2026hierarchical,
author = {Klavinskis-Whiting, Sebastian and King, Andrew J. and Harper, Nicol S.},
title = {{Hierarchical recurrent temporal prediction as a model of the mammalian dorsal visual pathway}},
journal = {PLoS computational biology},
year = {2026},
month = may,
volume = {22},
number = {5},
pages = {e1013138},
publisher = {PLOS},
issn = {1553-734X},
doi = {10.1371/journal.pcbi.1013138},
url = {https://doi.org/10.1371/journal.pcbi.1013138},
pmid = {42207826},
pmcid = {PMC13252843}
}

RIS

TY - JOUR
AU - Klavinskis-Whiting, Sebastian
AU - King, Andrew J.
AU - Harper, Nicol S.
TI - Hierarchical recurrent temporal prediction as a model of the mammalian dorsal visual pathway
T2 - PLoS computational biology
J2 - PLoS Comput Biol
PY - 2026
DA - 2026/05/28
VL - 22
IS - 5
SP - e1013138
SN - 1553-734X
PB - PLOS
DO - 10.1371/journal.pcbi.1013138
UR - https://doi.org/10.1371/journal.pcbi.1013138
LA - en
ER -

CSL-JSON

{
"id": "10.1371/journal.pcbi.1013138",
"type": "article-journal",
"title": "Hierarchical recurrent temporal prediction as a model of the mammalian dorsal visual pathway",
"container-title": "PLoS computational biology",
"author": [
{
"family": "Klavinskis-Whiting",
"given": "Sebastian"
},
{
"family": "King",
"given": "Andrew J."
},
{
"family": "Harper",
"given": "Nicol S."
}
],
"container-title-short": "PLoS Comput Biol",
"volume": "22",
"issue": "5",
"page": "e1013138",
"DOI": "10.1371/journal.pcbi.1013138",
"PMID": "42207826",
"PMCID": "PMC13252843",
"ISSN": "1553-734X",
"publisher": "PLOS",
"URL": "https://doi.org/10.1371/journal.pcbi.1013138",
"language": "en",
"issued": {
"date-parts": [
[
2026,
5,
28
]
]
}
}

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.1038/s41467-026-76098-y [code]
A single computational objective can produce specialization of streams in visual cortex.
Journal: Nature communications
In common: scikit-image, h5py, statsmodels, 6 other tools, 4 references
[2] doi:10.1016/j.isci.2026.116825 [code]
Social hierarchy shapes behavioral and transcriptional responses to chronic stress and ketamine in male mice.
Journal: iScience
In common: Pingouin, OpenCV, scikit-image, 8 other tools, systems
[3] doi:10.1038/s41593-026-02232-0 [code]
Entorhinal cortex represents task-relevant remote locations independently of CA1.
Journal: Nature neuroscience
In common: Pingouin, OpenCV, scikit-image, 8 other tools, systems
[4] doi:10.7554/elife.108408 [code]
Frequency and laminar profile of feature-specific visual activity revealed by interleaved EEG-fMRI.
Journal: eLife
In common: Pingouin, scikit-image, h5py, 5 other tools, systems, 2 references
[5] doi:10.1371/journal.pbio.3003915 [code]
Noise-invariant representations of sound emerge along the canonical cortical hierarchy.
Journal: PLoS biology
In common: OpenCV, scikit-image, h5py, 6 other tools, 2 references
[6] doi:10.1038/s41597-026-07248-6 [code]
A large-scale fMRI dataset for vision-language semantic association.
Journal: Scientific data
In common: OpenCV, scikit-image, h5py, 6 other tools, 2 references
[7] doi:10.1038/s41593-026-02314-z [code]
Low-dimensional population dynamics in the brainstem gate REM sleep.
Journal: Nature neuroscience
In common: Pingouin, OpenCV, h5py, 6 other tools, systems, 1 reference
[8] 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, h5py, statsmodels, 6 other tools, 1 reference
[9] doi:10.1038/s42003-026-10169-0 [code]
Shared representations in brains and models reveal a two-route cortical organization during scene perception.
Journal: Communications biology
In common: h5py, statsmodels, PyTorch, 5 other tools, 3 references
[10] doi:10.1016/j.isci.2026.117375 [code]
Motor priming is associated with widespread recruitment into neural ensembles and more rapid ensemble transitions.
Journal: iScience
In common: Pingouin, OpenCV, h5py, 7 other tools, systems

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.