OSCR

Entorhinal cortex represents task-relevant remote locations independently of CA1.

Code ↔ Paper

19 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 19 matches · 2 of them tie a paragraph to a whole file, not to given lines: weak matches, whose lines are not tinted
  1. [1] § Methods › Neural data preprocessing ↔ qualityMetrics/eaj_qualityParamValues.m, the whole file · a weak match · score 0.87 · BombCell, spatial decay, noise ratio, spikes missing, halfwidth, violations
  2. [2] § Methods › LFP analysis ↔ src/ripple_detection/literature_methods.py, lines 2506–2540 · score 0.83 · selected CA1 channel, Hilbert envelope, equiripple filtered, 125 Hz, trace, SWRs
  3. [3] § Methods › Neural data preprocessing ↔ qualityMetrics/bc_qualityParamValues.m, the whole file · a weak match · score 0.81 · spatial decay, noise ratio, spikes missing, BombCell, violations, global
  4. [4] § Methods › Single unit analysis ↔ Figures_published_edition.ipynb, lines 2248–2315 · score 0.80 · Spatial aperiodic cells, Border scores, Speed scores, angle, cm, locations
  5. [5] § Methods › Single unit analysis ↔ Yggdrasil/Spikes/spikes.py, lines 394–434 · score 0.77 · co fire, CA1 spike, MEC spike, 1–10 ms, monosynaptic, shuffled
  6. [6] § Results › CA1 decouples from MEC during nonlocal coding ↔ Yggdrasil/Spikes/spikes.py, lines 394–434 · score 0.68 · co fired, CA1 spike, MEC spike, 1–10 ms, shuffle, summed
  7. [7] § Methods › Position tracking ↔ Yggdrasil/Position/arena.py, lines 438–558 · score 0.67 · arena boundaries, track graph, videos, connected, behavior, edges
  8. [8] § Methods › Electrophysiology ↔ Yggdrasil/LFP/lfp.py, lines 58–140 · score 0.65 · 0.5–500 Hz, spikeGLX, Imec, streamed, preprocessing, cat
  9. [9] § Methods › LFP analysis ↔ src/ripple_detection/detectors/_lfp.py, lines 57–124 · score 0.60 · Hilbert envelope, ripple filtered, Hz, LFP, trace, score
  10. [10] § Methods › Position tracking ↔ src/track_linearization/core.py, lines 1036–1097 · score 0.60 · track graph, jumped, network, adjacent, linearized, behavior
  11. [11] § Methods › Decoding linearized position from population spiking ↔ Yggdrasil/Sequences/sequences.py, lines 232–271 · score 0.59 · acausal posterior, position bin, stationary, fragmented, decoder, probability
  12. [12] § Methods › Behavioral training ↔ Figures_published_edition.ipynb, lines 192–275 · score 0.57 · linear track, open field, ran, surgery, day, maze
  13. [13] § Methods › LFP analysis ↔ Yggdrasil/LFP/lfp.py, lines 58–140 · score 0.57 · hardware filter, downsampled, shift, reversed, stream, temporal
  14. [14] § Methods › Single unit analysis ↔ Yggdrasil/Sequences/sequences.py, lines 629–713 · score 0.56 · consecutive bins, position bin, shuffling, Gaussian, movement, cm
  15. [15] § Results › Nonlocal content represents task-relevant information ↔ extract_sequences_published_edition.ipynb, lines 646–709 · score 0.53 · naive Bayesian, Bayesian decoder, model, decoding, MEC, immobility
  16. [16] § Results › MEC represents nonlocal positions during immobility ↔ extract_sequences_published_edition.ipynb, lines 646–709 · score 0.52 · naive Bayesian decoder, immobility bout, classified, model, MEC, decoded
  17. [17] § Results › Characterizing cells involved in nonlocal coding ↔ Figures_published_edition.ipynb, lines 725–821 · score 0.52 · nonlocal intervals, active cells, Decoded position, ratio, fields, immobility
  18. [18] § Results › Nonlocal content represents task-relevant information ↔ Figures_published_edition.ipynb, lines 725–821 · score 0.52 · acausal posterior, active cells, decoded position, immobility, Nonlocal, location
  19. [19] § Results › MEC nonlocal coding occurs largely outside of SWRs ↔ src/ripple_detection/detectors/_hse.py, lines 111–167 · score 0.51 · high synchrony events, active units, Kernel, plus, traces, movement

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 · 3,709 lines · 195 KB · GPL-3.0 · 4 matches

  1. # %% [markdown]
  2. # Calculate and display all figures for paper
  4. # %%
  5. # imports
  6. from os.path import exists, join
  7. from os import chdir
  8. import numpy as np
  9. import pandas as pd
  10. from datetime import datetime
  11. from copy import deepcopy
  12. import matplotlib.pyplot as plt
  13. import matplotlib as mpl
  14. import re
  15. from scipy import signal, stats
  16. import time
  18. from sklearn.metrics.pairwise import cosine_similarity
  19. from sklearn.linear_model import LogisticRegression, LinearRegression
  20. from sklearn.model_selection import train_test_split
  21. from sklearn.metrics import confusion_matrix
  23. from Yggdrasil.Position.position import Position
  24. import Yggdrasil.Position.position as pos
  25. from Yggdrasil.Position.arena import Box, BoxWithObject, LinearTrack, DoubleYMaze
  26. from Yggdrasil.Task.task import LinearTrackTask, DoubleYMazeTask
  27. from Yggdrasil.Spikes.spikes import Spikes
  28. import Yggdrasil.Spikes.spikes as sp
  29. from Yggdrasil.Position.spatial_functions import *
  30. from Yggdrasil.Sequences.sequences import *
  31. from Yggdrasil.Sequences.plot import *
  32. from Yggdrasil.statistics import calc_lmm
  33. from Yggdrasil.LFP.lfp import LFP
  34. import Yggdrasil.LFP.lfp as lf
  35. from Yggdrasil.Electrodes.electrodes import Electrodes
  36. from Yggdrasil.utilities import get_starts
  37. from ripple_detection import Kay_ripple_detector, multiunit_HSE_detector, get_multiunit_population_firing_rate
  40. # reload modules without restarting the kernel
  41. %load_ext autoreload
  42. %autoreload 2
  43. import warnings
  44. warnings.filterwarnings("ignore", category=DeprecationWarning)
  46. #chdir(r'C:\Users\emily\OneDrive - Stanford\GitHub\GiocomoLab')
  47. chdir(r'C:\Users\Niflheim\Documents\GitHub\Giocomo')
  49. # %%
  50. # matplotlib style sheet
  51. from cycler import cycler
  52. from matplotlib import rcParams
  53. rcParams['lines.linewidth'] = 1
  54. rcParams['axes.linewidth'] = 1
  55. rcParams['font.size'] = 9
  56. rcParams['font.family'] = 'Arial'
  57. rcParams['figure.autolayout'] = True
  58. rcParams['pdf.fonttype'] = 42
  59. #rcParams['xtick.major.pad'] = 2
  60. #rcParams['ytick.major.pad'] = 2
  61. rcParams['xtick.bottom'] = False
  62. rcParams['ytick.left'] = False
  63. rcParams['axes.spines.top'] = False
  64. rcParams['axes.spines.right'] = False
  65. rcParams['axes.grid'] = False
  66. rcParams['lines.markersize'] = 2
  68. # %% [markdown]
  69. # 👉 Set the figure output path and path to the list of sessions and list of animals
  71. # %%
  72. # figure_path = r"C:\Users\emily\Dropbox\Giocomo Lab\WT Sequences\Figures"
  73. figure_path = r"Z:\WT_Sequences\Analysis"
  74. sessions = pd.read_csv('Z:/WT_Sequences/all_sessions_Xmaze.csv')
  75. subjects = pd.read_csv('Z:/WT_Sequences/subjects.csv')
  77. n_animals = len(subjects)
  78. mpl.rcParams['axes.prop_cycle'] = plt.cycler("color", plt.cm.gray(np.linspace(0, 0.7, n_animals)))
  79. colors = ['#B21F67','#DF5B25','#FFB000','#648FFF','#785EF0']
  81. # %%
  82. # specify animals
  83. animal_dict = {
  84. 'StickPin': 0,
  85. 'TopCoat': 1,
  86. 'TopHat': 2,
  87. 'AppleBottom': 3,
  88. 'BaggySweatpants': 4,
  89. 'StrappyReeboks': 5,
  90. 'Curie': 6,
  91. 'Lovelace': 7,
  92. 'Noether': 8,
  93. 'Franklin': 9,
  94. 'Lamarr': 10,
  95. 'Payne': 11
  96. }
  97. border_dict = {
  98. 'StickPin': 1780,
  99. 'TopCoat': 1900,
  100. 'TopHat': 1800,
  101. 'AppleBottom': (2100+2280)/2,
  102. 'BaggySweatpants': (1740+1660)/2,
  103. 'StrappyReeboks': 2040,
  104. 'Curie': [np.NaN, np.NaN, 2670, np.NaN],
  105. 'Lovelace': [np.NaN, 336/2*15, 348/2*15, np.NaN],
  106. 'Noether': [np.NaN, np.NaN, 364/2*15, np.NaN],
  107. 'Franklin': [np.NaN, np.NaN, 256/2*15, 268/2*15],
  108. 'Lamarr': [np.NaN, 276/2*15, 292/2*15, 300/2*15],
  109. 'Payne': [np.NaN, 254/2*15, 268/2*15, np.NaN]
  110. }
  111. # channels with highest theta power (and in MEC near units) across all sessions, 1 per animal
  112. theta_channel = {
  113. 'StickPin': 192,
  114. 'TopCoat': 183,
  115. 'TopHat': 177,
  116. 'AppleBottom': 159,
  117. 'BaggySweatpants': 159,
  118. 'StrappyReeboks': 155,
  119. 'Curie': 107,
  120. 'Lovelace': 59,
  121. 'Noether': 122,
  122. 'Franklin': 25,
  123. 'Lamarr': 161,
  124. 'Payne': 300
  125. }
  126. LOCAL_REMOTE_THRESH = 20
  128. # %% [markdown]
  129. # # Figure 1 & S1
  131. # %% [markdown]
  132. # ### Task performance over days
  134. # %%
  135. # set these values
  136. task_type = 'Single choice' #Single choice, Reversal or Cued
  137. n_sessions = 10 #10 days of Single choice or Reversal, 3 days of Cued
  139. # extract %corr into an array
  140. # NOTE: assumes list is in order by animal
  141. # this could be re-written to load values into a df where each col is an animal
  142. # and session indices are read from file names
  143. # but I am lazy and this works fine
  144. curr_animal = sessions['Animal'].iloc[0]
  145. animal_idx = 0
  146. session_idx = 0
  147. percent_correct = np.empty((n_sessions,n_animals))
  148. percent_correct[:] = np.NaN
  149. trials_per_min = np.empty((n_sessions,n_animals))
  150. trials_per_min[:] = np.NaN
  152. for i, row in sessions.iterrows():
  153. if row['Task'] == 'X Maze':
  154. behavior_output_path = join(row['Base_Directory'], 'Preprocessed_Data/Task')
  155. ecephys_path = join(row['Base_Directory'], 'Preprocessed_Data/Spikes')
  157. task_file = join(behavior_output_path, row['File']+'_task.txt')
  158. task = DoubleYMazeTask(name=task_file)
  159. start = 0
  160. end = task.trials['Start'].iloc[-1]
  162. if task.task_type == task_type:
  163. # advance to next column if next animal
  164. if not row['Animal'] == curr_animal:
  165. curr_animal = row['Animal']
  166. animal_idx += 1
  167. session_idx = 0
  168. percent_correct[session_idx, animal_idx] = task.percent_correct
  169. trials_per_min[session_idx, animal_idx] = task.ntrials/((end-start)/60)
  170. session_idx += 1
  172. # %%
  173. fig, ax = plt.subplots()
  174. ax.plot(np.arange(1,len(percent_correct)+1), percent_correct)
  175. ax.plot(np.arange(1,len(percent_correct)+1), np.nanmean(percent_correct,1), 'k', linewidth=3)
  176. ax.plot([0, 10], [50, 50], 'k--')
  177. ax.set_xlim([0, 10])
  178. ax.set_ylim([0, 100])
  179. ax.set_xlabel('Day')
  180. ax.set_ylabel('% Correct')
  181. fig.savefig(join(figure_path,'DY_percent_correct.pdf'), format='pdf')
  183. fig, ax = plt.subplots()
  184. ax.plot(np.arange(1, len(trials_per_min)+1), trials_per_min)
  185. ax.plot(np.arange(1,len(trials_per_min)+1), np.nanmean(trials_per_min,1), 'k', linewidth=3)
  186. ax.set_xlim([0, 10])
  187. ax.set_ylim([0, 8])
  188. ax.set_xlabel('Day')
  189. ax.set_ylabel('# Trials/Min')
  190. fig.savefig(join(figure_path, 'DY_trials_per_minute.pdf'), format='pdf')
  192. # %% [markdown]
  193. # ### Simultaneous cells recorded per day
  195. # %%
  196. sessions = pd.read_csv('Z:/WT_Sequences/all_sessions_full.csv')
  198. surgery_date_list = []
  199. n_sessions = 50 # placeholder, must be ># days post-op that recordings ran
  201. # get surgery dates
  202. for i, row in subjects.iterrows():
  203. # read surgery date from csv (Excel formats as MM/DD/YYYY)
  204. surgery_date_list.append(datetime.strptime(row['Surgery_Date'], '%m/%d/%Y'))
  206. # NOTE: assumes list is in order by animal
  207. # this could be re-written to load values into a df where each col is an animal
  208. # and session indices are read from file names
  209. # but I am lazy and this works fine
  210. curr_animal = sessions['Animal'].iloc[0]
  211. animal_idx = 0
  213. #n_units = np.zeros((n_sessions,n_animals))
  214. n_units_ca1 = np.full((n_sessions,n_animals), np.NaN)
  215. n_units_ca3 = np.full((n_sessions,n_animals), np.NaN)
  216. n_units_mec = np.full((n_sessions,n_animals), np.NaN)
  218. # count # of good & MUA units on each post-op day
  219. for i, row in sessions.iterrows():
  220. if not row['Recording_Error'] and not row['Position_Error'] and \
  221. (row['Task'] == 'X Maze' or row['Task'] == 'Linear Track' or \
  222. row['Epoch_Description']=='Environment A' or row['Epoch_Description']=='Open Field Test'):
  224. rec_date = datetime.strptime(re.findall('(\d{8})', row['File'])[0], '%Y%m%d')
  226. # advance to next column if next animal
  227. if not row['Animal'] == curr_animal:
  228. curr_animal = row['Animal']
  229. animal_idx += 1
  231. session_idx = (rec_date-surgery_date_list[animal_idx]).days
  232. #n_units[session_idx, animal_idx] += len(spikes.spikes)
  234. # subset by area
  235. if ("2023_spring" in row['Base_Directory']) or ("2024_winter" in row['Base_Directory']):
  236. if (row['Animal']=='Lamarr') and (int(row['Session'][-2:])>=8):
  237. spikes = Spikes(join(row['Base_Directory'], 'Preprocessed_Data/Spikes',
  238. row['File']+'_imec0_spikes.txt'))
  239. electrodes = Electrodes(join(row['Base_Directory'], 'Preprocessed_Data/Electrodes',
  240. row['File']+'_imec0_electrodes.txt'))
  241. MEC_channels = electrodes.subset_by_location(regions=['ENTm1','ENTm2','ENTm3','ENTm4','ENTm5','ENTm6'])
  242. spikes_copy = deepcopy(spikes)
  243. spikes_copy.subset_by_channel(MEC_channels)
  244. n_units_mec[session_idx, animal_idx] = len(spikes_copy.spikes)
  245. else:
  246. spikes = Spikes(join(row['Base_Directory'], 'Preprocessed_Data/Spikes',
  247. row['File']+'_imec0_spikes.txt'))
  248. electrodes = Electrodes(join(row['Base_Directory'], 'Preprocessed_Data/Electrodes',
  249. row['File']+'_imec0_electrodes.txt'))
  250. CA1_channels = electrodes.subset_by_location(regions=['CA1','alv','ccb'])
  251. spikes_copy = deepcopy(spikes)
  252. spikes_copy.subset_by_channel(CA1_channels)
  253. n_units_ca1[session_idx, animal_idx] = len(spikes_copy.spikes)
  254. CA3_channels = electrodes.subset_by_location(regions=['CA3','DG-po'])
  255. spikes_copy = deepcopy(spikes)
  256. spikes_copy.subset_by_channel(CA3_channels)
  257. n_units_ca3[session_idx, animal_idx] = len(spikes_copy.spikes)
  259. spikes = Spikes(join(row['Base_Directory'], 'Preprocessed_Data/Spikes',
  260. row['File']+'_imec1_spikes.txt'))
  261. electrodes = Electrodes(join(row['Base_Directory'], 'Preprocessed_Data/Electrodes',
  262. row['File']+'_imec1_electrodes.txt'))
  263. MEC_channels = electrodes.subset_by_location(regions=['ENTm1','ENTm2','ENTm3','ENTm4','ENTm5','ENTm6'])
  264. spikes_copy = deepcopy(spikes)
  265. spikes_copy.subset_by_channel(MEC_channels)
  266. n_units_mec[session_idx, animal_idx] = len(spikes_copy.spikes)
  267. else:
  268. spikes = Spikes(join(row['Base_Directory'], 'Preprocessed_Data/Spikes/g1',
  269. row['File']+'_imec0_spikes.txt'))
  270. electrodes = Electrodes(join(row['Base_Directory'], 'Preprocessed_Data/Electrodes',
  271. row['File']+'_electrodes.txt'))
  272. channels = electrodes.subset_by_location(regions=['Entorhinal area medial part dorsal zone'])
  273. spikes_copy = deepcopy(spikes)
  274. spikes_copy.subset_by_channel(channels)
  275. n_units_mec[session_idx, animal_idx] = len(spikes_copy.spikes)
  277. # %%
  278. # plot
  279. days = np.arange(0,50)
  281. fig, ax = plt.subplots(figsize=(10, 5))
  282. ind = ~np.isnan(np.asarray(n_units_ca1.astype(float)))
  283. for a in range(n_animals):
  284. ax.plot(days[ind[:,a]], n_units_ca1[ind[:,a], a])
  285. ax.plot(days, np.nanmean(n_units_ca1,1), 'k', linewidth=3)
  286. ax.set(xlabel='Days Post-op', ylabel='# CA1 Units')
  287. ax.set(xlim=[0, 50], ylim=[0, 250])
  288. fig.savefig(join(figure_path,'Units_per_day_CA1.pdf'), format='pdf')
  290. fig, ax = plt.subplots(figsize=(10, 5))
  291. ind = ~np.isnan(np.asarray(n_units_ca3.astype(float)))
  292. for a in range(n_animals):
  293. ax.plot(days[ind[:,a]], n_units_ca3[ind[:,a], a])
  294. ax.plot(days, np.nanmean(n_units_ca3,1), 'k', linewidth=3)
  295. ax.set(xlabel='Days Post-op', ylabel='# CA3 Units')
  296. ax.set(xlim=[0, 50], ylim=[0, 200])
  297. fig.savefig(join(figure_path,'Units_per_day_CA3.pdf'), format='pdf')
  299. fig, ax = plt.subplots(figsize=(10, 5))
  300. ind = ~np.isnan(np.asarray(n_units_mec.astype(float)))
  301. for a in range(n_animals):
  302. ax.plot(days[ind[:,a]], n_units_mec[ind[:,a], a])
  303. ax.plot(days, np.nanmean(n_units_mec,1), 'k', linewidth=3)
  304. ax.set(xlabel='Days Post-op', ylabel='# MEC Units')
  305. ax.set(xlim=[0, 50], ylim=[0, 525])
  306. fig.savefig(join(figure_path,'Units_per_day_MEC.pdf'), format='pdf')
  308. # %%
  309. days = np.arange(0,50)
  311. n_sites_ca1 = [np.NaN, np.NaN, 130, 202, 240, 152, 118, 208, np.NaN, np.NaN, np.NaN, np.NaN]
  312. n_sites_ca3 = [np.NaN, np.NaN, 30, 108, 102, 38, 92, 6, np.NaN, np.NaN, np.NaN, np.NaN]
  313. n_sites_mec = [144, 218, 364, 364, 304, 384, 364, 198, 202, 116, 194, 204] #192,
  315. fig, ax = plt.subplots(figsize=(10, 5))
  316. ind = ~np.isnan(np.asarray(n_units_ca1.astype(float)))
  317. for a in range(n_animals):
  318. ax.plot(days[ind[:,a]], (n_units_ca1[ind[:,a], a]+n_units_ca3[ind[:,a], a]) \
  319. /(n_sites_ca1[a]+n_sites_ca3[a]))
  320. ax.plot(days, np.nanmean(n_units_ca1/n_sites_ca1,1), 'k', linewidth=3)
  321. ax.set(xlabel='Days Post-op', ylabel='# Units/Site in CA1')
  322. ax.set(xlim=[0, 50], ylim=[0, 1.75], yticks=np.arange(0,1.75,0.25))
  323. fig.savefig(join(figure_path,'Units_per_site_CA1.pdf'), format='pdf')
  325. fig, ax = plt.subplots(figsize=(10, 5))
  326. ind = ~np.isnan(np.asarray(n_units_mec.astype(float)))
  327. for a in range(n_animals):
  328. ax.plot(days[ind[:,a]], n_units_mec[ind[:,a], a]/n_sites_mec[a])
  329. ax.plot(days, np.nanmean(n_units_mec/n_sites_mec,1), 'k', linewidth=3)
  330. ax.set(xlabel='Days Post-op', ylabel='# Units/Site in MEC')
  331. ax.set(xlim=[0, 50], ylim=[0, 1.75], yticks=np.arange(0,1.75,0.25))
  332. fig.savefig(join(figure_path,'Units_per_site_MEC.pdf'), format='pdf')
  334. # %% [markdown]
  335. # ### Example heatmaps & rasters
  337. # %%
  338. bin_cm = 2
  340. for i, row in sessions.iloc[[109]].iterrows(): #118
  341. print(row['File'])
  342. for probe in range(2):
  343. spikes = Spikes(name=join(row['Base_Directory'], 'Preprocessed_Data/Spikes',
  344. row['File']+'_imec'+str(probe)+'_spikes.txt'))
  345. position = Position(name=join(row['Base_Directory'], 'Preprocessed_Data/Position',
  346. row['File']+'_position.txt'))
  348. if probe==0:
  349. electrodes = Electrodes(join(row['Base_Directory'], 'Preprocessed_Data/Electrodes',
  350. row['File']+'_imec0_electrodes.txt'))
  351. CA1_channels = electrodes.subset_by_location(regions=['CA1','alv','ccb'])
  352. spikes.subset_by_channel(CA1_channels)
  353. else:
  354. electrodes = Electrodes(join(row['Base_Directory'], 'Preprocessed_Data/Electrodes',
  355. row['File']+'_imec1_electrodes.txt'))
  357. MEC_channels = electrodes.subset_by_location(regions=['ENTm1','ENTm2','ENTm3','ENTm4','ENTm5','ENTm6'])
  358. spikes.subset_by_channel(MEC_channels)
  360. for u in range(len(spikes.spikes)):
  361. if probe==1:
  362. if u in [26,234,250]: #CA1 0,28,30; MEC 26,234,250
  363. fr_map, xbins, ybins, max_fr = calc_fr_map(spikes.spikes[u], position, bin_cm, smooth=True)
  364. if max_fr>1:
  365. print(f'Unit {u} Channel {spikes.spike_channel[u]} Max FR {max_fr}')
  366. f = position.plot_position(spikes=spikes.spikes[u])
  367. f.savefig(join(figure_path,f'Lovelace_DY01_MEC_{u}_raster.pdf'), format='pdf')
  368. f = plot_heatmap(fr_map, xbins, ybins)
  369. f.savefig(join(figure_path,f'Lovelace_DY01_MEC_{u}_heatmap.pdf'), format='pdf')
  371. # %% [markdown]
  372. # ### Speed distribution over days and locations
  374. # %%
  375. n_days = 10
  376. n_bins = 50
  377. vel_bins = np.arange(0,n_bins*2+2,2)
  379. speed_hist = np.zeros(n_bins)
  380. session_count = 0
  381. speed_hist_by_days = np.zeros((n_days,n_bins))
  382. day_counts = np.zeros(n_days)
  383. speed_hist_by_mouse = np.zeros((n_animals,n_bins))
  384. animal_counts = np.zeros(n_animals)
  386. for i, row in sessions.iterrows():
  387. session_idx = int(row['Session'][-2:])-1
  388. if not row['Recording_Error'] and not row['Position_Error'] and row['Task'] == 'X Maze' and session_idx<10:
  389. position = Position(name=join(row['Base_Directory'], 'Preprocessed_Data/Position',
  390. row['File']+'_position.txt'))
  391. this_session_hist = (pd.cut(position.position.Velocity, bins=vel_bins).value_counts().sort_index()/len(position.position.Velocity)).values
  392. speed_hist += this_session_hist
  393. session_count += 1
  394. speed_hist_by_days[session_idx, :] += this_session_hist
  395. day_counts[session_idx] += 1
  396. animal_idx = animal_dict[row['Animal']]
  397. speed_hist_by_mouse[animal_idx,:] += this_session_hist
  398. animal_counts[animal_idx] += 1
  400. # %%
  401. fig, ax = plt.subplots()
  402. for a in range(n_animals):
  403. ax.stairs(speed_hist_by_mouse[a]/animal_counts[a]*100, vel_bins)
  404. ax.stairs(speed_hist/session_count*100, vel_bins, color='k', linewidth=3)
  405. ax.set_xlim([0,80])
  406. ax.set_ylim([0,35])
  407. ax.set_xlabel('Velocity (cm/s)')
  408. ax.set_ylabel('% Time Spent')
  409. fig.savefig(join(figure_path,'Velocity_hist.pdf'), format='pdf')
  411. # %%
  412. ## Heatmap of speed over 3 examples days
  413. row = sessions.iloc[155]
  414. position = Position(name=join(row['Base_Directory'], 'Preprocessed_Data/Position',
  415. row['File']+'_position.txt'))
  416. speed_by_2D_location, x_bins, y_bins = position.calc_feature_map(bin_cm=2)
  417. f = plot_heatmap(speed_by_2D_location, x_bins, y_bins)
  418. f.savefig(join(figure_path,f'Payne_DY01_speed_heatmap.pdf'), format='pdf')
  419. print(np.nanmax(speed_by_2D_location))
  421. row = sessions.iloc[222]
  422. position = Position(name=join(row['Base_Directory'], 'Preprocessed_Data/Position',
  423. row['File']+'_position.txt'))
  424. speed_by_2D_location, x_bins, y_bins = position.calc_feature_map(bin_cm=2)
  425. f = plot_heatmap(speed_by_2D_location, x_bins, y_bins)
  426. f.savefig(join(figure_path,f'TopCoat_DY05_speed_heatmap.pdf'), format='pdf')
  427. print(np.nanmax(speed_by_2D_location))
  429. row = sessions.iloc[9]
  430. position = Position(name=join(row['Base_Directory'], 'Preprocessed_Data/Position',
  431. row['File']+'_position.txt'))
  432. speed_by_2D_location, x_bins, y_bins = position.calc_feature_map(bin_cm=2)
  433. f = plot_heatmap(speed_by_2D_location, x_bins, y_bins)
  434. f.savefig(join(figure_path,f'AppleBottom_DY10_speed_heatmap.pdf'), format='pdf')
  435. print(np.nanmax(speed_by_2D_location))
  437. # %%
  438. ## immobility by location
  439. velocity_thresh = 2 # cm/s
  440. n_sessions = 10
  441. frac_immob = np.full((n_sessions,n_animals), np.NaN)
  442. frac_immob_at_reward = np.full((n_sessions,n_animals), np.NaN)
  443. frac_immob_outside_reward = np.full((n_sessions,n_animals), np.NaN)
  444. frac_immob_at_decision = np.full((n_sessions,n_animals), np.NaN)
  446. for i, row in sessions.iterrows():
  447. session_idx = int(row['Session'][-2:])-1
  448. animal_idx = animal_dict[row['Animal']]
  449. if not row['Recording_Error'] and not row['Position_Error'] and row['Task'] == 'X Maze' and session_idx<10:
  450. position = Position(join(row['Base_Directory'], 'Preprocessed_Data/Position', row['File']+'_position.txt'))
  451. dist_to_reward = position.arena.get_dist_to_nearest_poke(position.position.X, position.position.Y)
  452. frac_immob[session_idx, animal_idx] = np.mean((position.position.Velocity.values<velocity_thresh))
  453. frac_immob_at_reward[session_idx, animal_idx] = sum((position.position.Velocity.values<velocity_thresh) & (dist_to_reward<10))/sum(dist_to_reward<10)
  454. frac_immob_outside_reward[session_idx, animal_idx] = sum((position.position.Velocity.values<velocity_thresh) & (dist_to_reward>=10))/sum(dist_to_reward>=10)
  455. dist_to_decision = position.arena.get_dist_to_decision(position.position.X, position.position.Y)
  456. frac_immob_at_decision[session_idx, animal_idx] = sum((position.position.Velocity.values<velocity_thresh) & (dist_to_decision<10))/sum(dist_to_decision<10)
  458. # %%
  459. fig, ax = plt.subplots()
  460. days = np.arange(1, 11)
  461. ind = ~np.isnan(np.asarray(frac_immob.astype(float)))
  462. for a in range(n_animals):
  463. ax.plot(days[ind[:,a]], frac_immob[ind[:,a], a]*100, linewidth=0.5)
  464. ax.plot(days, np.nanmean(frac_immob,1)*100, 'k', linewidth=2)
  465. ax.set(xlabel='Day', ylabel='% Time Immobile', ylim=[0, 100])
  466. fig.savefig(join(figure_path,f'Immobility_time_over_days.pdf'), format='pdf')
  468. fig, ax = plt.subplots()
  469. days = np.arange(1, 11)
  470. ind = ~np.isnan(np.asarray(frac_immob_at_reward.astype(float)))
  471. for a in range(n_animals):
  472. ax.plot(days[ind[:,a]], frac_immob_at_reward[ind[:,a], a]*100, linewidth=0.5)
  473. ax.plot(days, np.nanmean(frac_immob_at_reward,1)*100, 'k', linewidth=2)
  474. ax.set(xlabel='Day', ylabel='% Time Immobile at Reward', ylim=[0, 100])
  475. fig.savefig(join(figure_path,f'Immobility_time_at_reward_over_days.pdf'), format='pdf')
  477. fig, ax = plt.subplots()
  478. days = np.arange(1, 11)
  479. ind = ~np.isnan(np.asarray(frac_immob_outside_reward.astype(float)))
  480. for a in range(n_animals):
  481. ax.plot(days[ind[:,a]], frac_immob_outside_reward[ind[:,a], a]*100, linewidth=0.5)
  482. ax.plot(days, np.nanmean(frac_immob_outside_reward,1)*100, 'k', linewidth=2)
  483. ax.set(xlabel='Day', ylabel='% Time Immobile >10cm from Reward', ylim=[0, 100])
  484. fig.savefig(join(figure_path,f'Immobility_time_outside_reward_over_days.pdf'), format='pdf')
  486. # %% [markdown]
  487. # ### No overrepresentation: example heatmap of summation of normalized firing rates, summation by location over days
  489. # %%
  490. bin_cm = 2
  492. for i, row in sessions.iloc[[109]].iterrows():
  493. if not row['Recording_Error'] and not row['Position_Error'] and row['Task'] == 'X Maze':
  494. print(row['File'])
  495. ca1_spikes = None
  496. mec_spikes = None
  497. position = Position(name=join(row['Base_Directory'], 'Preprocessed_Data/Position',
  498. row['File']+'_position.txt'))
  500. if ("2023_spring" in row['Base_Directory']) or ("2024_winter" in row['Base_Directory']):
  501. if (row['Animal']=='Lamarr') and (int(row['Session'][-2:])>=8):
  502. mec_spikes = Spikes(name=join(row['Base_Directory'], 'Preprocessed_Data/Spikes',
  503. row['File']+'_imec0_spikes.txt'))
  504. electrodes = Electrodes(join(row['Base_Directory'], 'Preprocessed_Data/Electrodes',
  505. row['File']+'_imec0_electrodes.txt'))
  507. MEC_channels = electrodes.subset_by_location(regions=['ENTm1','ENTm2','ENTm3','ENTm4','ENTm5','ENTm6'])
  508. mec_spikes.subset_by_channel(MEC_channels)
  509. else:
  510. ca1_spikes = Spikes(join(row['Base_Directory'], 'Preprocessed_Data/Spikes',
  511. row['File']+'_imec0_spikes.txt'))
  512. electrodes = Electrodes(join(row['Base_Directory'], 'Preprocessed_Data/Electrodes',
  513. row['File']+'_imec0_electrodes.txt'))
  514. CA1_channels = electrodes.subset_by_location(regions=['CA1','alv','ccb'])
  515. ca1_spikes.subset_by_channel(CA1_channels)
  517. mec_spikes = Spikes(join(row['Base_Directory'], 'Preprocessed_Data/Spikes',
  518. row['File']+'_imec1_spikes.txt'))
  519. electrodes = Electrodes(join(row['Base_Directory'], 'Preprocessed_Data/Electrodes',
  520. row['File']+'_imec1_electrodes.txt'))
  521. MEC_channels = electrodes.subset_by_location(regions=['ENTm1','ENTm2','ENTm3','ENTm4','ENTm5','ENTm6'])
  522. mec_spikes.subset_by_channel(MEC_channels)
  523. else:
  524. mec_spikes = Spikes(join(row['Base_Directory'], 'Preprocessed_Data/Spikes/g1',
  525. row['File']+'_imec0_spikes.txt'))
  526. electrodes = Electrodes(join(row['Base_Directory'], 'Preprocessed_Data/Electrodes',
  527. row['File']+'_electrodes.txt'))
  528. channels = electrodes.subset_by_location(regions=['Entorhinal area medial part dorsal zone'])
  529. mec_spikes.subset_by_channel(channels)
  531. fr_map, xbins, ybins, max_fr = calc_fr_map(mec_spikes.spikes[0], position, bin_cm, smooth=True)
  532. summed_mec_map = np.zeros(fr_map.shape)
  533. summed_ca1_map = np.zeros(fr_map.shape)
  535. if ca1_spikes is not None:
  536. for u in range(len(ca1_spikes.spikes)):
  537. fr_map, xbins, ybins, max_fr = calc_fr_map(ca1_spikes.spikes[u], position, bin_cm, smooth=False)
  538. summed_ca1_map += fr_map/np.max(fr_map)
  539. f = plot_heatmap(summed_ca1_map, xbins, ybins)
  540. f.savefig(join(figure_path,'Lovelace_DY01_CA1_heatmap.pdf'), format='pdf')
  542. if mec_spikes is not None:
  543. for u in range(len(mec_spikes.spikes)):
  544. fr_map, xbins, ybins, max_fr = calc_fr_map(mec_spikes.spikes[u], position, bin_cm, smooth=False)
  545. summed_mec_map += fr_map/np.max(fr_map)
  546. f = plot_heatmap(summed_mec_map, xbins, ybins)
  547. f.savefig(join(figure_path,'Lovelace_DY01_MEC_heatmap.pdf'), format='pdf')
  549. # %%
  550. # calculate normalized ratemap in each linearized position bin
  551. # using linearized data will make it easier to combine across sessions
  552. # arena borders were re-drawn for each of the 4 cohorts
  553. # thus each cohort has slightly different linear bins,
  554. # but these can be easily combined across segments
  555. sessions = pd.read_csv(r'Z:/WT_Sequences/all_sessions_Xmaze.csv')
  556. n_sessions = 20
  557. n_animals = len(animal_dict)
  558. norm_fr_map = [[[] for _ in range(n_sessions)] for _ in range(n_animals)]
  560. for _, row in sessions.iterrows():
  561. if not row['Recording_Error'] and not row['Position_Error'] and int(row['Session'][-2:])<=20:
  562. seq = Sequences(join(row['Base_Directory'],
  563. 'Preprocessed_Data/Sequences',
  564. row['File']+'_MEC_nonlocal_immobility_sequences.txt'),
  565. load_objects=True)
  566. animal_idx = animal_dict[row['Animal']]
  567. sess_idx = int(row['Session'][-2:])-1
  569. # get linearized bin position of animal at each binned timestamp
  570. linear_to_bins = np.asarray(seq.classifier.classifier.environments[0].place_bin_centers_nodes_df_['linear_position'])
  571. linear_pos = np.asarray(seq.binned_data.position["Linear"])
  572. animal_bins = np.argmin(np.abs(linear_to_bins[:,np.newaxis]-linear_pos), axis=0)
  573. animal_bins = animal_bins.astype(dtype=float)
  574. animal_bins[np.isnan(linear_pos)] = np.NaN
  575. # same as calc_occupancy, but for 1D position
  576. bins = np.arange(np.nanmax(animal_bins)+1)
  577. occupancy, _ = np.histogram(animal_bins, bins=bins)
  578. # convert to seconds
  579. occupancy = occupancy.astype(float)
  580. occupancy *= seq.position.us_per_frame/10**3
  581. # same as get_spike_map, but for 1D position
  582. spike_map = np.zeros((len(seq.binned_data.spikes), len(bins)-1))
  583. norm_fr_map[animal_idx][sess_idx] = np.zeros(len(bins)-1)
  584. for u in np.arange(len(seq.binned_data.spikes)):
  585. spikes_per_bin = []
  586. for n_spk in np.arange(1, np.max(seq.binned_data.spikes[u,:]) + 1):
  587. spikes_per_bin += animal_bins[seq.binned_data.spikes[u,:] == n_spk].tolist() * int(n_spk)
  588. spike_map[u,:], _ = np.histogram(spikes_per_bin, bins)
  589. unit_fr_map = get_fr_map(spike_map[u,:], occupancy)
  590. norm_fr_map[animal_idx][sess_idx] += unit_fr_map/np.nanmax(unit_fr_map)
  591. norm_fr_map[animal_idx][sess_idx] /= len(seq.binned_data.spikes)
  592. print(row['File'])
  594. with open(join(figure_path,'MEC_norm_fr_map.pkl'), 'wb') as file:
  595. pickle.dump(norm_fr_map, file)
  597. # %%
  598. with open(join(figure_path,'MEC_norm_fr_map.pkl'), 'rb') as file:
  599. norm_fr_map = pickle.load(file)
  601. norm_fr_map[6][6] = norm_fr_map[6][6][:123]
  602. norm_fr_map[6][4] = norm_fr_map[6][6][:123]
  603. norm_fr_map[6][18] = norm_fr_map[6][6][:123]
  605. # sum over sessions
  606. all_sessions_fr = np.full((n_animals,n_sessions,123), np.NaN)
  607. for a in range(n_animals):
  608. for s in range(n_sessions):
  609. if len(norm_fr_map[a][s])>0:
  610. # mask bins where occupancy was 0
  611. norm_fr_map[a][s][norm_fr_map[a][s]==0] = np.NaN
  612. all_sessions_fr[a,s,0:len(norm_fr_map[a][s])] = norm_fr_map[a][s]
  613. sum_over_sessions_fr = np.nanmean(all_sessions_fr, axis=1)
  615. # mask out bins with too few sessions
  616. empty_bins = []
  617. for a in range(n_animals):
  618. empty_bins.extend(np.where(np.isnan(sum_over_sessions_fr[a]))[0])
  619. n_empty, _ = np.histogram(empty_bins, bins=np.arange(0,124))
  620. keep_locs = np.arange(0,123)[n_empty<=1]
  622. # plot over spatial bins
  623. spatial_bins = np.arange(1,len(keep_locs)+1)
  624. sum_over_sessions_fr = sum_over_sessions_fr[:,keep_locs]
  625. fig, ax = plt.subplots()
  626. for a in range(n_animals):
  627. ax.plot(spatial_bins[~np.isnan(sum_over_sessions_fr[a,:])], \
  628. sum_over_sessions_fr[a,~np.isnan(sum_over_sessions_fr[a,:])], linewidth=0.5)
  629. ax.plot(spatial_bins, np.nanmean(sum_over_sessions_fr, axis=0), 'k', linewidth=2)
  630. ax.set(xlabel='Spatial Bin', ylabel='Summed Normalized FR', ylim=[0, 0.5])
  631. fig.savefig(join(figure_path,f'MEC_norm_fr_per_spatial_bin.pdf'), format='pdf')
  633. # %%
  634. # repeat for CA1
  635. sessions = pd.read_csv(r'Z:/WT_Sequences/all_sessions_Xmaze_CA1.csv')
  636. norm_fr_map = [[[] for _ in range(n_sessions)] for _ in range(n_animals)]
  638. for _, row in sessions.iterrows():
  639. if not row['Recording_Error'] and not row['Position_Error'] and int(row['Session'][-2:])<=20:
  640. seq = Sequences(join(row['Base_Directory'],
  641. 'Preprocessed_Data/Sequences',
  642. row['File']+'_CA1_immobility_sequences.txt'),
  643. load_objects=True)
  644. animal_idx = animal_dict[row['Animal']]
  645. sess_idx = int(row['Session'][-2:])-1
  647. # get linearized bin position of animal at each binned timestamp
  648. linear_to_bins = np.asarray(seq.classifier.classifier.environments[0].place_bin_centers_nodes_df_['linear_position'])
  649. linear_pos = np.asarray(seq.binned_data.position["Linear"])
  650. animal_bins = np.argmin(np.abs(linear_to_bins[:,np.newaxis]-linear_pos), axis=0)
  651. animal_bins = animal_bins.astype(dtype=float)
  652. animal_bins[np.isnan(linear_pos)] = np.NaN
  653. # same as calc_occupancy, but for 1D position
  654. bins = np.arange(np.nanmax(animal_bins)+1)
  655. occupancy, _ = np.histogram(animal_bins, bins=bins)
  656. # convert to seconds
  657. occupancy = occupancy.astype(float)
  658. occupancy *= seq.position.us_per_frame/10**3
  659. # same as get_spike_map, but for 1D position
  660. spike_map = np.zeros((len(seq.binned_data.spikes), len(bins)-1))
  661. norm_fr_map[animal_idx][sess_idx] = np.zeros(len(bins)-1)
  662. for u in np.arange(len(seq.binned_data.spikes)):
  663. spikes_per_bin = []
  664. for n_spk in np.arange(1, np.max(seq.binned_data.spikes[u,:]) + 1):
  665. spikes_per_bin += animal_bins[seq.binned_data.spikes[u,:] == n_spk].tolist() * int(n_spk)
  666. spike_map[u,:], _ = np.histogram(spikes_per_bin, bins)
  667. unit_fr_map = get_fr_map(spike_map[u,:], occupancy)
  668. norm_fr_map[animal_idx][sess_idx] += unit_fr_map/np.nanmax(unit_fr_map)
  669. norm_fr_map[animal_idx][sess_idx] /= len(seq.binned_data.spikes)
  671. with open(join(figure_path,'CA1_norm_fr_map.pkl'), 'wb') as file:
  672. pickle.dump(norm_fr_map, file)
  674. # %%
  675. with open(join(figure_path,'CA1_norm_fr_map.pkl'), 'rb') as file:
  676. norm_fr_map = pickle.load(file)
  678. norm_fr_map[6][6] = norm_fr_map[6][6][:123]
  679. norm_fr_map[6][4] = norm_fr_map[6][6][:123]
  680. norm_fr_map[6][18] = norm_fr_map[6][6][:123]
  682. # sum over sessions
  683. all_sessions_fr = np.full((n_animals,n_sessions,123), np.NaN)
  684. for a in range(n_animals):
  685. for s in range(n_sessions):
  686. if len(norm_fr_map[a][s])>0:
  687. # mask bins where occupancy was 0
  688. norm_fr_map[a][s][norm_fr_map[a][s]==0] = np.NaN
  689. all_sessions_fr[a,s,0:len(norm_fr_map[a][s])] = norm_fr_map[a][s]
  690. sum_over_sessions_fr = np.nanmean(all_sessions_fr, axis=1)
  692. # mask out bins with too few sessions
  693. empty_bins = []
  694. for a in range(n_animals):
  695. empty_bins.extend(np.where(np.isnan(sum_over_sessions_fr[a]))[0])
  696. n_empty, _ = np.histogram(empty_bins, bins=np.arange(0,124))
  697. keep_locs = np.arange(0,123)[n_empty<=6+1]
  699. # plot over spatial bins
  700. spatial_bins = np.arange(1,len(keep_locs)+1)
  701. sum_over_sessions_fr = sum_over_sessions_fr[:,keep_locs]
  702. fig, ax = plt.subplots()
  703. for a in range(n_animals):
  704. ax.plot(spatial_bins[~np.isnan(sum_over_sessions_fr[a,:])], \
  705. sum_over_sessions_fr[a,~np.isnan(sum_over_sessions_fr[a,:])], linewidth=0.5)
  706. ax.plot(spatial_bins, np.nanmean(sum_over_sessions_fr, axis=0), 'k', linewidth=2)
  707. ax.set(xlabel='Spatial Bin', ylabel='Summed Normalized FR', ylim=[0, 0.5])
  708. fig.savefig(join(figure_path,f'CA1_norm_fr_per_spatial_bin.pdf'), format='pdf')
  710. # %% [markdown]
  711. # # Figure 2
  713. # %% [markdown]
  714. # ### Example decode
  716. # %%
  717. # open example
  718. sessions = pd.read_csv('Z:/WT_Sequences/all_sessions_Xmaze.csv')
  719. row = sessions.iloc[133] #109 #141
  720. seq = Sequences(join(row['Base_Directory'], 'Preprocessed_Data/Sequences',
  721. row['File']+'_MEC_immobility_sequences.txt'),
  722. load_objects=True)
  723. print(row['File'])
  725. # %%
  726. # plot spikes, position, decoded position, decode distance, and speed for snippet of data
  727. # subset data
  728. start_time = 1260
  729. end_time = 1290
  730. start_idx = np.argmin(np.abs(seq.binned_data.timestamps-start_time))
  731. end_idx = np.argmin(np.abs(seq.binned_data.timestamps-end_time))
  732. subset_timestamps = seq.binned_data.timestamps[start_idx:end_idx]
  734. f, ax = plt.subplots(4, 1, constrained_layout=True, sharex=False,
  735. gridspec_kw={"height_ratios": [3,3,1,1]})
  737. # order cells by peak location along track and plot raster
  738. cell_order = order_cells_by_field(seq.classifier.classifier)
  739. neuron_ind, spike_time_ind = np.nonzero(seq.binned_data.spikes[cell_order,start_idx:end_idx])
  740. ax[0].scatter(1000*(subset_timestamps[spike_time_ind]-start_time),
  741. neuron_ind, color='black', zorder=1,
  742. marker='|', linewidth=1)
  743. n_active_cells = len(np.unique(neuron_ind))
  745. # label spikes during non-local immobility
  746. time_slice = slice(start_time,end_time)
  747. dist = seq.calc_dist_from_decode(ts_to_include=seq.binned_data.position.loc[time_slice].index)
  748. ds_factor = 5
  749. binned_ds_dist = np.nanmean(np.reshape(dist[:-1],(-1,ds_factor)), axis=1)
  750. nonlocal_neuron_ind = np.empty(0, dtype=int)
  751. nonlocal_intervals = pos.filtered_timestamps_to_intervals(np.arange(0,3000), binned_ds_dist>20)
  752. for i in range(4,len(nonlocal_intervals)):
  753. neuron_ind, spike_time_ind = np.nonzero(seq.binned_data.spikes[:,
  754. int(start_idx+nonlocal_intervals[i][0]*5):int(start_idx+nonlocal_intervals[i][1]*5)])
  755. # ax[1].scatter(1000*(subset_timestamps[nonlocal_intervals[i][0]*5+spike_time_ind]-start_time),
  756. # neuron_ind, color='darkorange', zorder=1,
  757. # marker='|', linewidth=1)
  758. nonlocal_neuron_ind = np.hstack((nonlocal_neuron_ind, neuron_ind))
  759. nonlocal_neuron_ind = np.unique(nonlocal_neuron_ind)
  761. neuron_ind, spike_time_ind = np.nonzero(seq.binned_data.spikes[nonlocal_neuron_ind,start_idx:end_idx])
  762. # map onto cell order
  763. for n in range(len(neuron_ind)):
  764. neuron_ind[n] = np.argmin(np.abs(cell_order-nonlocal_neuron_ind[neuron_ind[n]]))
  765. # ax[0].scatter(1000*(subset_timestamps[spike_time_ind]-start_time),
  766. # neuron_ind, color='darkorange', zorder=1,
  767. # marker='|', linewidth=1)
  768. ax[0].set(xlim=[0,30000], ylabel='Neuron')
  770. # add classifier
  771. results = seq.classifier.classifier_results.sel(time=time_slice)
  772. time = (results.time-start_time) * 1000
  773. max_time = time.max()
  775. cmap = copy.copy(plt.cm.get_cmap('bone_r'))
  776. cmap.set_bad(color="lightgrey", alpha=1.0)
  777. (
  778. results
  779. .assign_coords(time=time)
  780. .acausal_posterior.sum("state")
  781. .plot(
  782. x="time",
  783. y="position",
  784. robust=True,
  785. add_colorbar=False,
  786. zorder=0,
  787. rasterized=True,
  788. cmap=cmap,
  789. ax=ax[1]
  790. )
  791. )
  792. seq_position = seq.binned_data.position['Linear'].loc[time_slice]
  793. max_position = int(
  794. np.ceil(seq.binned_data.position['Linear'].max()))
  795. ax[1].plot(time, seq_position, linestyle="--", linewidth=2,
  796. color="magenta", clip_on=False)
  798. rtc.plot_graph_as_1D(seq.position.arena.track_graph,
  799. edge_spacing=seq.position.arena.edge_spacing,
  800. ax=ax[1], axis="y", other_axis_start=max_time+50)
  802. # overlay: distance between decode and real positions
  803. ax[2].plot(pd.DataFrame(binned_ds_dist).interpolate())
  804. ax[2].set(xlim=[0,3000], ylabel='Distance from decoded position (cm)')
  805. # highlight nonlocal times
  806. for i in range(len(nonlocal_intervals)):
  807. ax[2].axvspan(nonlocal_intervals[i][0], nonlocal_intervals[i][1], alpha=0.4, zorder=2, color='darkorange')
  809. # overlay: velocity
  810. velocity = seq.binned_data.position['Velocity'].iloc[start_idx:end_idx]
  811. ds_factor = 5
  812. binned_ds_vel = np.nanmean(np.reshape(velocity.values,(-1,ds_factor)), axis=1)
  813. ax[3].plot(pd.DataFrame(binned_ds_vel).interpolate())
  814. ax[3].set(xlim=[0,3000], ylabel='Speed (cm/s)')
  815. # highlight immobility times
  816. immobility_intervals = pos.filtered_timestamps_to_intervals(np.arange(0,3000), binned_ds_vel<2)
  817. for i in range(len(immobility_intervals)):
  818. ax[3].axvspan(immobility_intervals[i][0], immobility_intervals[i][1], alpha=0.3, zorder=2)
  819. f.savefig(join(figure_path,'Noether_DY02_decode_snippet.pdf'), format='pdf')
  821. print(f"{len(nonlocal_neuron_ind)} of {n_active_cells} cells active during nonlocal immobility")
  823. # %% [markdown]
  824. # ### Decoder error vs velocity bins & # of cells
  826. # %%
  827. distance, velocity, session, animal, n_units = [np.array([]) for _ in range(5)]
  829. for i, row in sessions.iterrows():
  830. if not row['Recording_Error'] and not row['Position_Error'] and row['Task'] == 'X Maze':
  831. # load any sequences file, doesn't matter which interval (we're just using the
  832. # classifier & binned data associated with the file, which is the same for
  833. # immobility, movement, and SWR sequences)
  834. seq = Sequences(join(row['Base_Directory'], 'Preprocessed_Data/Sequences',
  835. row['File']+'_MEC_immobility_sequences.txt'),
  836. load_objects=True)
  838. # exclude times when the decode is fragmented
  839. cr = seq.classifier.classifier_results
  840. non_frag = cr.acausal_posterior.sum('position').argmax('state')!=1
  841. non_frag_ts = non_frag.coords['time'].values[np.where(non_frag)]
  843. # calculate distance between animal and decoded position
  844. dist = seq.calc_dist_from_decode(ts_to_include=non_frag_ts)
  846. # remove times when position couldn't be tracked
  847. vel = np.asarray(seq.binned_data.position["Velocity"].loc[non_frag_ts])
  848. dist = dist[~np.isnan(vel)]
  849. vel = vel[~np.isnan(vel)]
  851. # append
  852. distance = np.append(distance, dist)
  853. velocity = np.append(velocity, vel)
  854. session = np.append(session, [int(row['Session'][-2:])-1] * len(dist))
  855. animal = np.append(animal, [animal_dict[row['Animal']]] * len(dist))
  856. n_units = np.append(n_units, [seq.binned_data.spikes.shape[0]] * len(dist))
  858. df = pd.DataFrame({'Velocity': velocity, 'Distance': distance,
  859. 'Session': session, 'Animal': animal,
  860. 'N_Units': n_units})
  862. # %%
  863. df.to_csv(join(figure_path,'Decoded_distance_binned.csv'), index=False)
  864. #df = pd.read_csv(join(figure_path,'Decoded_distance_binned.csv'))
  866. # %%
  867. heatmap, xedges, yedges = np.histogram2d(df['Velocity'], df['Distance'], bins=[35,26])
  868. norm_heatmap = heatmap.T/np.sum(heatmap, axis=1)
  869. extent = [yedges[0], yedges[-2], xedges[0], xedges[-19]]
  870. cmap = mpl.cm.get_cmap("viridis").copy()
  871. #cmap.set_under('k')
  873. fig, ax = plt.subplots()
  874. hmap = ax.imshow(np.transpose(norm_heatmap[0:-2,0:-19]), extent=extent, origin='lower', cmap=cmap, norm=mpl.colors.LogNorm(0.01, 0.05))
  875. fig.colorbar(hmap, extend='max')
  876. ax.set(ylabel='Speed (cm/s)', xlabel='Distance from Animal (cm)')
  877. fig.savefig(join(figure_path,'Decode_distance_vs_velocity_heatmap.pdf'), format='pdf')
  879. # %%
  880. heatmap, xedges, yedges = np.histogram2d(df['Velocity'], df['Distance'], bins=[35,26])
  881. extent = [xedges[0], xedges[-19], yedges[0], yedges[-1]]
  883. fig, ax = plt.subplots()
  884. h, edges = np.histogram(df.Distance[df.Velocity<2], bins=np.arange(0,yedges[-1]+5,5))
  885. ax.stairs(h/h.sum()*100, edges, color='grey')
  886. frac_local_immob = np.sum(h[0:4]/h.sum())*100
  887. h, edges = np.histogram(df.Distance[df.Velocity>=2], bins=np.arange(0,yedges[-1]+5,5))
  888. ax.stairs(h/h.sum()*100, edges, color='black')
  889. frac_local_move = np.sum(h[0:4]/h.sum())*100
  890. ax.set(xlim=[edges[0],edges[-2]+1], ylim=[0,50], xlabel='Distance from decoded position', ylabel='% of Time')
  891. print(f"{frac_local_move:0.2f}% local during movement and {frac_local_immob:0.2f}% local during immobility")
  892. fig.savefig(join(figure_path,'Decode_distance_immob_vs_move.pdf'), format='pdf')
  894. fig, ax = plt.subplots()
  895. h, edges = np.histogram(df.Distance[df.Velocity<2], bins=np.arange(0,yedges[-1]+5,5))
  896. ax.stairs(h/h.sum()*100, edges, color='grey')
  897. h, edges = np.histogram(df.Distance[df.Velocity>=2], bins=np.arange(0,yedges[-1]+5,5))
  898. ax.stairs(h/h.sum()*100, edges, color='black')
  899. ax.set(xlim=[20,edges[-2]+1], ylim=[0,4], xlabel='Distance from decoded position', ylabel='% of Time')
  900. fig.savefig(join(figure_path,'Decode_distance_immob_vs_move_inset.pdf'), format='pdf')
  902. # %% [markdown]
  903. # ### Example raster, 1D & 2D decodes of immobility bouts with nonlocal coding
  905. # %%
  906. row = sessions.iloc[133]
  907. seq = Sequences(join(row['Base_Directory'], 'Preprocessed_Data/Sequences',
  908. row['File']+'_MEC_immobility_sequences.txt'),
  909. load_objects=True)
  910. print(seq.name)
  912. # %%
  913. seq.plot_sequence(indices=[1,33,41,305], figure_path=f"{figure_path}/Noether_DY02_MEC")
  915. # %%
  916. # repeat raster with color coding by local/nonlocal
  917. # first, load the interval
  918. start_time = seq.intervals.intervals.iloc[1].start_time
  919. end_time = seq.intervals.intervals.iloc[1].end_time
  920. start_idx = np.argmin(np.abs(seq.binned_data.timestamps-start_time))
  921. end_idx = np.argmin(np.abs(seq.binned_data.timestamps-end_time))
  922. subset_timestamps = seq.binned_data.timestamps[start_idx:end_idx]
  924. # order cells by peak location along track and plot raster
  925. cell_order, max_bin = order_cells_by_field(seq.classifier.classifier)
  927. # %%
  928. # then, determine which spatial bins are local & nonlocal
  929. cr = seq.classifier.classifier_results.sel(time=slice(start_time, end_time))
  930. map_position_ind = cr.sum("state").acausal_posterior.argmax("position").values
  931. plt.plot(map_position_ind)
  932. print(np.unique(map_position_ind))
  934. # %%
  935. # then, identify cells that have significant activity at the nonlocal (contributing_cells) and local (noncontributing_cells) positions
  936. # set decoded_bin manually for each interval
  937. # this method finds cells which have firing at the target spatial bin(s) that is >95th percentile of their firing across all spatial bins
  938. # also tried >mean across all spatial bins
  939. # also tried nonspatial approach: identify cells whose firing during this particular interval > mean or 95th percentile of shuffle across all immobility
  940. # this mainly identified high FR cells, as low FR cells are less likely to participate in any given interval
  942. # 1: local 59, nonlocal 29, 92-96
  943. # 33: local 0, nonlocal 27
  944. # 41: local 49-51 (58, 60); nonlocal 29, 78-79
  945. # 305: local 121, nonlocal 101,102
  947. field_locs, field_peaks, field_sizes, field_spacing = seq.find_fields()
  949. pf = seq.classifier.classifier.place_fields_[('', 0)].sel(neuron=cell_order)
  950. field_locs_ordered = [field_locs[i] for i in cell_order]
  951. n_cells = pf.shape[1]
  952. decoded_bin = [59]
  953. local_cells = []
  954. for u in range(n_cells):
  955. # if np.any(pf.sel(neuron=u).values[decoded_bin] - np.nanpercentile(pf.sel(neuron=u).values, 95) > 0): # Old Method
  956. for f in field_locs_ordered[u]:
  957. if np.any(np.isin(decoded_bin, f)):
  958. local_cells.append(u)
  959. local_cells = np.unique(np.asarray(local_cells))
  961. decoded_bin = [29,92,93,94,95,96]
  962. nonlocal_cells = []
  963. for u in range(n_cells):
  964. for f in field_locs_ordered[u]:
  965. if np.any(np.isin(decoded_bin, f)):
  966. nonlocal_cells.append(u)
  967. nonlocal_cells = np.unique(np.asarray(nonlocal_cells))
  969. # %%
  970. f, ax = plt.subplots(2, 1, constrained_layout=True, sharex=False,
  971. gridspec_kw={"height_ratios": [5,5]})
  973. # label spikes during non-local immobility
  974. time_slice = slice(start_time,end_time)
  975. dist = seq.calc_dist_from_decode(ts_to_include=seq.binned_data.position.loc[time_slice].index)
  976. ds_factor = 5
  977. binned_ds_dist = np.nanmean(np.reshape(dist,(-1,ds_factor)), axis=1)
  978. nonlocal_neuron_ind = np.empty(0, dtype=int)
  979. nonlocal_intervals = pos.filtered_timestamps_to_intervals(np.arange(0,len(binned_ds_dist)), binned_ds_dist>20)
  981. interval_spikes = seq.binned_data.spikes[cell_order,start_idx:end_idx]
  982. neuron_ind, spike_time_ind = np.nonzero(interval_spikes)
  983. ax[0].scatter(1000*(subset_timestamps[spike_time_ind]-start_time),
  984. neuron_ind, color='grey', zorder=1,
  985. marker='|', linewidth=1)
  987. # plot spikes from cells contributing to decode (over all time) in black
  988. neuron_ind, spike_time_ind = np.nonzero(interval_spikes[nonlocal_cells, :])
  989. neuron_ind = nonlocal_cells[neuron_ind]
  990. ax[0].scatter(1000*(subset_timestamps[spike_time_ind]-start_time),
  991. neuron_ind, color='black', zorder=1,
  992. marker='|', linewidth=1)
  993. neuron_ind, spike_time_ind = np.nonzero(interval_spikes[local_cells, :])
  994. neuron_ind = local_cells[neuron_ind]
  995. ax[0].scatter(1000*(subset_timestamps[spike_time_ind]-start_time),
  996. neuron_ind, color='black', zorder=1,
  997. marker='|', linewidth=1)
  999. # plot spikes contributing to nonlocal decode in cyan during non-local times
  1000. for i in range(len(nonlocal_intervals)):
  1001. interval_spikes = seq.binned_data.spikes[cell_order,int(start_idx+nonlocal_intervals[i][0]*5):int(start_idx+nonlocal_intervals[i][1]*5)]
  1002. neuron_ind_nl, spike_time_ind_nl = np.nonzero(interval_spikes[nonlocal_cells, :])
  1003. neuron_ind_nl = nonlocal_cells[neuron_ind_nl]
  1004. ax[0].scatter(1000*(subset_timestamps[nonlocal_intervals[i][0]*5+spike_time_ind_nl]-start_time),
  1005. neuron_ind_nl, color='cyan', zorder=1,
  1006. marker='|', linewidth=2)
  1008. # plot spikes contributing to local decode in magenta during local times
  1009. # before and after last nonlocal interval
  1010. interval_spikes_l = seq.binned_data.spikes[cell_order,start_idx:int(start_idx+nonlocal_intervals[0][0]*5)]
  1011. neuron_ind_l, spike_time_ind_l = np.nonzero(interval_spikes_l[local_cells, :])
  1012. neuron_ind_l = local_cells[neuron_ind_l]
  1013. ax[0].scatter(1000*(subset_timestamps[spike_time_ind_l]-start_time),
  1014. neuron_ind_l, color='magenta', zorder=1,
  1015. marker='|', linewidth=2)
  1016. interval_spikes_l = seq.binned_data.spikes[cell_order,int(start_idx+nonlocal_intervals[-1][1]*5):end_idx]
  1017. neuron_ind_l, spike_time_ind_l = np.nonzero(interval_spikes_l[local_cells, :])
  1018. neuron_ind_l = local_cells[neuron_ind_l]
  1019. ax[0].scatter(1000*(subset_timestamps[nonlocal_intervals[i][1]*5+spike_time_ind_l]-start_time),
  1020. neuron_ind_l, color='magenta', zorder=1,
  1021. marker='|', linewidth=2)
  1023. # between nonlocal lintervals
  1024. if len(nonlocal_intervals)>1:
  1025. for i in range(1,len(nonlocal_intervals)):
  1026. interval_spikes_l = seq.binned_data.spikes[cell_order,int(start_idx+nonlocal_intervals[i-1][1]*5):int(start_idx+nonlocal_intervals[i][0]*5)]
  1027. neuron_ind_l, spike_time_ind_l = np.nonzero(interval_spikes_l[local_cells, :])
  1028. neuron_ind_l = local_cells[neuron_ind_l]
  1029. ax[0].scatter(1000*(subset_timestamps[nonlocal_intervals[i-1][1]*5+spike_time_ind_l]-start_time),
  1030. neuron_ind_l, color='magenta', zorder=1,
  1031. marker='|', linewidth=2)
  1033. ax[0].set(xlim=[0,int((end_time-start_time)*1000)], ylabel='Neuron')
  1034. f.savefig(join(figure_path,'Noether_DY02_decode_snippet_pf_percentile_1.pdf'), format='pdf')
  1036. # %% [markdown]
  1037. # ### Nonlocal decoding per immobility bout and total % of immobility nonlocal
  1039. # %%
  1040. sessions = pd.read_csv('Z:/WT_Sequences/all_sessions_Xmaze.csv')
  1041. immob_data_by_session, immob_data_by_seq = combine_metrics_over_sessions(sessions, animal_dict, 'MEC_immobility')
  1042. local_data_by_session, local_data_by_seq = combine_metrics_over_sessions(sessions, animal_dict, 'MEC_local_immobility')
  1043. nonlocal_data_by_session, nonlocal_data_by_seq = combine_metrics_over_sessions(sessions, animal_dict, 'MEC_nonlocal_immobility')
  1044. sessions = pd.read_csv('Z:/WT_Sequences/all_sessions_Xmaze_CA1.csv')
  1045. SWR_data_by_session, SWR_data_by_seq = combine_metrics_over_sessions(sessions, animal_dict, 'CA1_SWR')
  1046. HSE_data_by_session, HSE_data_by_seq = combine_metrics_over_sessions(sessions, animal_dict, 'CA1_HSE')
  1048. # %%
  1049. plot_seq_metrics_over_events(immob_data_by_seq.Animal, immob_data_by_seq.Percent_non_local_content, figure_path,
  1050. 'perc_immob_content_nl_per_bout', '% of Each Immobility Bout Decoded Non-Locally', [0,1])
  1052. # %%
  1053. print(f"{np.nanmean(immob_data_by_session['non_local']):.02f} +/- {stats.sem(np.nanmean(immob_data_by_session['non_local'],axis=1)):.02f}" +\
  1054. f"% of immobility bouts have nonlocal coding,\n" +\
  1055. f"and {np.nanmean(immob_data_by_session['frac_non_local']):.02f} +/- {stats.sem(np.nanmean(immob_data_by_session['frac_non_local'],axis=1)):.02f}" +\
  1056. f"% of time during immobility "+\
  1057. "is spent representing nonlocal positions.")
  1059. # %%
  1060. plot_seq_metrics_over_sessions(n_animals, immob_data_by_session['non_local'], figure_path,
  1061. 'immob_with_nl', '% Immobility Bouts with NonLocal Decode', [0,100])
  1063. plot_seq_metrics_over_sessions(n_animals, immob_data_by_session['frac_non_local'], figure_path,
  1064. 'perc_immob_content_nl', '% Immobility Time with NonLocal Decode', [0,100])
  1066. # %% [markdown]
  1067. # ### Duration of non-local coding
  1069. # %%
  1070. mean_dur = np.mean(nonlocal_data_by_seq.groupby(['Animal'])['Duration'].mean().values)
  1071. sem_dur = stats.sem(nonlocal_data_by_seq.groupby(['Animal'])['Duration'].mean().values)
  1072. print(f"Periods of non-local coding during immobility last on average "+
  1073. f"{mean_dur:.03f} +/- {sem_dur:.03f} seconds.")
  1075. # %%
  1076. mean_dur = np.mean(SWR_data_by_seq.groupby(['Animal'])['Duration'].mean().values)
  1077. sem_dur = stats.sem(SWR_data_by_seq.groupby(['Animal'])['Duration'].mean().values)
  1078. print(f"SWRs last on average "+
  1079. f"{mean_dur:.03f} +/- {sem_dur:.03f} seconds.")
  1081. # %%
  1082. plot_seq_metrics_over_events(nonlocal_data_by_seq.Animal, nonlocal_data_by_seq.Duration, figure_path,
  1083. 'nl_dur', 'Length of Non-Local Decodes (s)', [0,1])
  1085. # %% [markdown]
  1086. # ### Immobility duration vs decode distance
  1088. # %%
  1089. # Non-local coding is more common during longer bouts of immobility
  1090. plot_seq_metrics_comparison_over_events(immob_data_by_seq, 'Duration', 'Non_local',
  1091. figure_path, 'immob_local_vs_nonlocal_duration', 'Duration (s)', [0,10], [0,4])
  1093. # %% [markdown]
  1094. # # Figure S3
  1096. # %% [markdown]
  1097. # ### Confusion Matrix: actual vs decoded position in MEC during movement
  1099. # %%
  1100. actual_pos = np.array([])
  1101. decoded_pos = np.array([])
  1103. start_time = time.time()
  1104. for _, row in sessions.iterrows():
  1105. sess_idx = int(row['Session'][-2:])-1
  1106. if not row['Recording_Error'] and not row['Position_Error'] and \
  1107. row['Task'] == 'X Maze' and sess_idx<=10:
  1108. animal_idx = animal_dict[row['Animal']]
  1109. binned_data = Binned_Data(join(row['Base_Directory'],
  1110. 'Preprocessed_Data/Sequences/Binned_Data',
  1111. row['File']+'_MEC_binned_data.txt'))
  1112. classifier = Classifier_Model(join(row['Base_Directory'],
  1113. 'Preprocessed_Data/Sequences/Classifier_Model',
  1114. row['File']+'_MEC_classifier_model.txt'))
  1116. # get movement times
  1117. move_ts = binned_data.position.index[binned_data.position.Velocity>=2].values
  1119. # get current position
  1120. actual_pos = np.append(actual_pos, binned_data.position.Linear[move_ts])
  1122. # get decoded linearized bins and map to linearized position
  1123. cr = classifier.classifier_results.sel(time=move_ts)
  1124. decoded_bins = cr.sum("state").acausal_posterior.argmax("position").values
  1125. bins_to_linear = classifier.classifier.environments[0].place_bin_centers_nodes_df_['linear_position']
  1126. decoded_pos = np.append(decoded_pos, bins_to_linear[decoded_bins])
  1128. print(f"{row['File']} {time.time()-start_time}")
  1130. df = pd.DataFrame({'actual_position': actual_pos,
  1131. 'MEC_decode_position': decoded_pos})
  1133. # %%
  1134. df.to_csv(join(figure_path,'position_vs_MEC_decode_during_movement.csv'), index=False)
  1136. # %%
  1137. # heatmap of actual vs decoded positions
  1138. heatmap, xedges, yedges = np.histogram2d(df.loc[:,'actual_position'], df.loc[:,'MEC_decode_position'], bins=[26,26])
  1139. extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]
  1140. cmap = mpl.cm.get_cmap("viridis").copy()
  1142. fig, ax = plt.subplots()
  1143. norm_heatmap = heatmap.T/np.sum(heatmap, axis=1)
  1144. hmap = ax.imshow(norm_heatmap, extent=extent, origin='lower', cmap=cmap, vmax=0.5)
  1145. fig.colorbar(hmap)
  1146. ax.set(xlabel='Actual Linearized Position (cm)', ylabel='MEC Decoded Linearized Position (cm)')
  1147. fig.savefig(join(figure_path,'position_vs_MEC_decode_during_movement_positionnorm.pdf'), format='pdf')
  1149. # %% [markdown]
  1150. # ### non-local coding is common across days
  1152. # %%
  1153. days = np.arange(1,11)
  1154. fig, ax = plt.subplots()
  1155. for _, animal_metric in enumerate(immob_data_by_session['non_local']):
  1156. animal_metric = animal_metric[:10]
  1157. ax.plot(days[~np.isnan(animal_metric)], animal_metric[~np.isnan(animal_metric)])
  1158. ax.plot(np.arange(1,11), np.nanmean(immob_data_by_session['non_local'].T[:10],1), 'k', linewidth=2)
  1159. ax.set(xlabel='Day', ylabel='% Immobility Bouts with NonLocal Decode', ylim=[0, 100])
  1160. fig.savefig(join(figure_path,f'immob_with_nl_over_days.pdf'), format='pdf')
  1162. fig, ax = plt.subplots()
  1163. for _, animal_metric in enumerate(immob_data_by_session['frac_non_local']):
  1164. animal_metric = animal_metric[:10]
  1165. ax.plot(days[~np.isnan(animal_metric)], animal_metric[~np.isnan(animal_metric)])
  1166. ax.plot(np.arange(1,11), np.nanmean(immob_data_by_session['frac_non_local'].T[:10],1), 'k', linewidth=2)
  1167. ax.set(xlabel='Day', ylabel='% Immobility Time with NonLocal Decode', ylim=[0, 100])
  1168. fig.savefig(join(figure_path,f'perc_immob_content_nl.pdf'), format='pdf')
  1170. # %% [markdown]
  1171. # ### pie chart: % of local & non-local immobility classified as stationary, continuous, fragmented
  1173. # %%
  1174. f, ax = plt.subplots(1,3)
  1176. ax[0].pie([np.sum(immob_data_by_seq['Percent_stationary']*immob_data_by_seq['Duration']), \
  1177. np.sum(immob_data_by_seq['Percent_continuous']*immob_data_by_seq['Duration']), \
  1178. np.sum(immob_data_by_seq['Percent_fragmented']*immob_data_by_seq['Duration'])], \
  1179. labels=['stat','cont','frag'])
  1180. ax[1].pie([np.sum(local_data_by_seq['Percent_stationary']*local_data_by_seq['Duration']), \
  1181. np.sum(local_data_by_seq['Percent_continuous']*local_data_by_seq['Duration']), \
  1182. np.sum(local_data_by_seq['Percent_fragmented']*local_data_by_seq['Duration'])], \
  1183. labels=['stat','cont','frag'])
  1184. ax[2].pie([np.sum(nonlocal_data_by_seq['Percent_stationary']*nonlocal_data_by_seq['Duration']), \
  1185. np.sum(nonlocal_data_by_seq['Percent_continuous']*nonlocal_data_by_seq['Duration']), \
  1186. np.sum(nonlocal_data_by_seq['Percent_fragmented']*nonlocal_data_by_seq['Duration'])], \
  1187. labels=['stat','cont','frag'])
  1188. print(f"{np.sum(immob_data_by_seq['Percent_fragmented']*immob_data_by_seq['Duration'])/np.sum(immob_data_by_seq['Duration']):.02f}% of all immobility, " +
  1189. f"{np.sum(local_data_by_seq['Percent_fragmented']*local_data_by_seq['Duration'])/np.sum(local_data_by_seq['Duration']):.02f}% of local immobility, " +
  1190. f"and {np.sum(nonlocal_data_by_seq['Percent_fragmented']*nonlocal_data_by_seq['Duration'])/np.sum(nonlocal_data_by_seq['Duration']):.02f}% of nonlocal immobility " +
  1191. "is fragmented.")
  1193. f.savefig(join(figure_path,'classification.pdf'), format='pdf')
  1195. # %% [markdown]
  1196. # ### local vs nonlocal quality check: variance of posteriors
  1198. # %%
  1199. # match durations
  1200. nonlocal_equivalent = []
  1201. for a in local_data_by_seq.Animal.unique():
  1202. start_idx = np.min(np.where(nonlocal_data_by_seq.Animal==a))
  1203. for __, dur in local_data_by_seq.loc[local_data_by_seq.Animal==a].Duration.items():
  1204. idx = np.abs(nonlocal_data_by_seq.loc[nonlocal_data_by_seq.Animal==a].Duration - dur).argmin()
  1205. nonlocal_equivalent.append(start_idx + idx)
  1207. is_non_local = []
  1208. is_non_local.extend([False] * len(local_data_by_seq))
  1209. is_non_local.extend([True] * len(nonlocal_data_by_seq.iloc[nonlocal_equivalent]))
  1210. local_nonlocal_equivalent_data_by_seq = pd.concat([local_data_by_seq, nonlocal_data_by_seq.iloc[nonlocal_equivalent]], ignore_index=True)
  1211. local_nonlocal_equivalent_data_by_seq.Non_local = is_non_local
  1213. local_nonlocal_equivalent_data_by_seq['Posterior_density_spread_percent_corrected'] = \
  1214. local_nonlocal_equivalent_data_by_seq.Posterior_density_spread_percent/local_nonlocal_equivalent_data_by_seq.Spatial_coverage_percent
  1216. # %%
  1217. plot_seq_metrics_comparison_over_events(local_nonlocal_equivalent_data_by_seq, 'Posterior_density_spread_percent_corrected', 'Non_local', figure_path,
  1218. 'local_vs_nonlocal_posterior_density_corrected', '% Track Covered by Posterior Density >0.95% Per % Track Covered', [0,100])
  1220. # %%
  1221. # these data structures are generated in Figure 3
  1222. with open(join(figure_path,'nonlocal_fields_data_by_seq2.pkl'), 'rb') as file:
  1223. nonlocal_fields_data_by_seq2 = pickle.load(file)
  1224. with open(join(figure_path,'local_fields_data_by_seq2.pkl'), 'rb') as file:
  1225. local_fields_data_by_seq2 = pickle.load(file)
  1227. # %%
  1228. # add posterior variance and spatial coverage
  1229. posterior_variance = []
  1230. spatial_coverage = []
  1232. for _, row in sessions.iterrows():
  1233. if not row['Recording_Error'] and not row['Position_Error'] and \
  1234. row['Task'] == 'X Maze' and int(row['Session'][-2:])<=10:
  1235. seq = Sequences(join(row['Base_Directory'], 'Preprocessed_Data/Sequences',
  1236. row['File']+'_MEC_local_immobility_sequences.txt'))
  1237. posterior_variance.extend(seq.stats.Posterior_density_spread_percent)
  1238. spatial_coverage.extend(seq.stats.Spatial_coverage_percent)
  1239. local_fields_data_by_seq2['Posterior_density_spread_percent'] = posterior_variance
  1240. local_fields_data_by_seq2['Spatial_coverage_percent'] = spatial_coverage
  1242. # repeat for nonlocal
  1243. posterior_variance = []
  1244. spatial_coverage = []
  1246. for _, row in sessions.iterrows():
  1247. if not row['Recording_Error'] and not row['Position_Error'] and \
  1248. row['Task'] == 'X Maze' and int(row['Session'][-2:])<=10:
  1249. seq = Sequences(join(row['Base_Directory'], 'Preprocessed_Data/Sequences',
  1250. row['File']+'_MEC_nonlocal_immobility_sequences.txt'))
  1251. posterior_variance.extend(seq.stats.Posterior_density_spread_percent)
  1252. spatial_coverage.extend(seq.stats.Spatial_coverage_percent)
  1253. nonlocal_fields_data_by_seq2['Posterior_density_spread_percent'] = posterior_variance
  1254. nonlocal_fields_data_by_seq2['Spatial_coverage_percent'] = spatial_coverage
  1256. # %%
  1257. # match durations
  1258. nonlocal_equivalent = []
  1259. for a in local_fields_data_by_seq2.Animal.unique():
  1260. start_idx = np.min(np.where(nonlocal_fields_data_by_seq2.Animal==a))
  1261. for __, dur in local_fields_data_by_seq2.loc[local_fields_data_by_seq2.Animal==a].Duration.items():
  1262. idx = np.abs(nonlocal_fields_data_by_seq2.loc[nonlocal_fields_data_by_seq2.Animal==a].Duration - dur).argmin()
  1263. nonlocal_equivalent.append(start_idx + idx)
  1265. is_non_local = []
  1266. is_non_local.extend([False] * len(local_fields_data_by_seq2))
  1267. is_non_local.extend([True] * len(nonlocal_fields_data_by_seq2.iloc[nonlocal_equivalent]))
  1268. local_nonlocal_fields_equivalent_data_by_seq = pd.concat([local_fields_data_by_seq2, nonlocal_fields_data_by_seq2.iloc[nonlocal_equivalent]], ignore_index=True)
  1269. local_nonlocal_fields_equivalent_data_by_seq['Non_local'] = is_non_local
  1271. local_nonlocal_fields_equivalent_data_by_seq['Posterior_density_spread_percent_corrected'] = \
  1272. local_nonlocal_fields_equivalent_data_by_seq.Posterior_density_spread_percent/local_nonlocal_fields_equivalent_data_by_seq.Spatial_coverage_percent
  1274. # %%
  1275. # calculate mean over 2D bins
  1276. local_fields_df = local_nonlocal_fields_equivalent_data_by_seq[~local_nonlocal_fields_equivalent_data_by_seq.Non_local]
  1277. local_cells = pd.cut(local_fields_df['Local_percent_cells']*100, bins=np.arange(0,110,10))
  1278. nonlocal_cells = pd.cut(local_fields_df['Nonlocal_percent_cells']*100, bins=np.arange(0,110,10))
  1279. local_fields_ppv = local_fields_df.groupby([local_cells, nonlocal_cells])['Posterior_density_spread_percent_corrected'].mean().unstack().T
  1281. # repeat for nonlocal coding events
  1282. nonlocal_fields_df = local_nonlocal_fields_equivalent_data_by_seq[local_nonlocal_fields_equivalent_data_by_seq.Non_local]
  1283. local_cells = pd.cut(nonlocal_fields_df['Local_percent_cells']*100, bins=np.arange(0,110,10))
  1284. nonlocal_cells = pd.cut(nonlocal_fields_df['Nonlocal_percent_cells']*100, bins=np.arange(0,110,10))
  1285. nonlocal_fields_ppv = nonlocal_fields_df.groupby([local_cells, nonlocal_cells])['Posterior_density_spread_percent_corrected'].mean().unstack().T
  1287. # %%
  1288. # heatmap over coding events, local vs nonlocal: x=% local cells, y=% nonlocal cells, z/color= mean corrected posterior probability variance
  1289. cmap = mpl.cm.get_cmap("viridis").copy()
  1290. cmap.set_bad(color='lightgrey')
  1292. fig, ax = plt.subplots()
  1293. sns.heatmap(local_fields_ppv, fmt=".2f", cmap=cmap, vmax=100, cbar_kws={'label': '% Track Covered by Posterior Density'})
  1294. plt.gca().invert_yaxis()
  1295. ax.set(xlabel='% Local Units', ylabel='% Non-Local Units', xticklabels=np.arange(0,100,10), yticklabels=np.arange(0,100,10))
  1296. fig.savefig(join(figure_path,'ppv_heatmap_local.pdf'), format='pdf')
  1298. fig, ax = plt.subplots()
  1299. sns.heatmap(nonlocal_fields_ppv, fmt=".2f", cmap=cmap, vmax=100, cbar_kws={'label': '% Track Covered by Posterior Density'})
  1300. plt.gca().invert_yaxis()
  1301. ax.set(xlabel='% Local Units', ylabel='% Non-Local Units', xticklabels=np.arange(0,100,10), yticklabels=np.arange(0,100,10))
  1302. fig.savefig(join(figure_path,'ppv_heatmap_nonlocal.pdf'), format='pdf')
  1304. cmap = mpl.cm.get_cmap("seismic").copy()
  1305. cmap.set_bad(color='lightgrey')
  1307. fig, ax = plt.subplots()
  1308. sns.heatmap((local_fields_ppv-nonlocal_fields_ppv), fmt=".2f", vmin=-50, vmax=50, cmap=cmap, cbar_kws={'label': '% Track Covered by Posterior Density: Local Events Minus Non-Local Events'})
  1309. plt.gca().invert_yaxis()
  1310. ax.set(xlabel='% Local Units', ylabel='% Non-Local Units', xticklabels=np.arange(0,100,10), yticklabels=np.arange(0,100,10))
  1311. fig.savefig(join(figure_path,'ppv_heatmap_delta.pdf'), format='pdf')
  1313. # %% [markdown]
  1314. # ### local vs nonlocal quality check: changes in # spikes, # active cells, or spatial info
  1316. # %%
  1317. is_non_local = []
  1318. is_non_local.extend([False] * len(local_data_by_seq))
  1319. is_non_local.extend([True] * len(nonlocal_data_by_seq))
  1320. local_nonlocal_data_by_seq = pd.concat([local_data_by_seq, nonlocal_data_by_seq], ignore_index=True)
  1321. local_nonlocal_data_by_seq.Non_local = is_non_local
  1323. # %%
  1324. # match durations
  1325. nonlocal_equivalent = []
  1326. for a in local_data_by_seq.Animal.unique():
  1327. start_idx = np.min(np.where(nonlocal_data_by_seq.Animal==a))
  1328. for __, dur in local_data_by_seq.loc[local_data_by_seq.Animal==a].Duration.items():
  1329. idx = np.abs(nonlocal_data_by_seq.loc[nonlocal_data_by_seq.Animal==a].Duration - dur).argmin()
  1330. nonlocal_equivalent.append(start_idx + idx)
  1332. is_non_local = []
  1333. is_non_local.extend([False] * len(local_data_by_seq))
  1334. is_non_local.extend([True] * len(nonlocal_data_by_seq.iloc[nonlocal_equivalent]))
  1335. local_nonlocal_equivalent_data_by_seq = pd.concat([local_data_by_seq, nonlocal_data_by_seq.iloc[nonlocal_equivalent]], ignore_index=True)
  1336. local_nonlocal_equivalent_data_by_seq.Non_local = is_non_local
  1338. # %%
  1339. plot_seq_metrics_comparison_over_events(local_nonlocal_equivalent_data_by_seq, 'Percent_units', 'Non_local', figure_path,
  1340. 'local_vs_nonlocal_active_units', '% Active Units', [0,100])
  1341. plot_seq_metrics_comparison_over_events(local_nonlocal_equivalent_data_by_seq, 'FR', 'Non_local', figure_path,
  1342. 'local_vs_nonlocal_FR_matched', 'FR', [0,50])
  1343. plot_seq_metrics_comparison_over_events(local_nonlocal_equivalent_data_by_seq, 'Spatial_information', 'Non_local', figure_path,
  1344. 'local_vs_nonlocal_spatial_info_matched', 'Spatial Information (bits/sec)', [0,1.5])
  1346. # %% [markdown]
  1347. # ### Non-local content with different decoders
  1349. # %%
  1350. # decoder trained on movement + immobility
  1351. sessions = pd.read_csv('Z:/WT_Sequences/all_sessions_Xmaze.csv')
  1352. immob_data_by_session_alt_decoder, _ = combine_metrics_over_sessions(sessions, animal_dict, 'MEC_immobility', modifier='_cv')
  1354. print(f"{np.nanmean(immob_data_by_session_alt_decoder['non_local']):.02f} +/- {stats.sem(np.nanmean(immob_data_by_session_alt_decoder['non_local'],axis=1)):.02f}" +\
  1355. f"% of immobility bouts have nonlocal coding,\n" +\
  1356. f"and {np.nanmean(immob_data_by_session_alt_decoder['frac_non_local']):.02f} +/- {stats.sem(np.nanmean(immob_data_by_session_alt_decoder['frac_non_local'],axis=1)):.02f}" +\
  1357. f"% of time during immobility "+\
  1358. "is spent representing nonlocal positions.")
  1360. plot_seq_metrics_over_sessions(n_animals, immob_data_by_session_alt_decoder['non_local'], figure_path,
  1361. 'immob_with_nl_alt_decoder_cv', '% Immobility Bouts with NonLocal Decode', [0,100])
  1363. plot_seq_metrics_over_sessions(n_animals, immob_data_by_session_alt_decoder['frac_non_local'], figure_path,
  1364. 'perc_immob_content_nl_alt_decoder_cv', '% Immobility Time with NonLocal Decode', [0,100])
  1366. # %%
  1367. # plot position and decoded position for snippet of data
  1368. # open example
  1369. row = sessions.iloc[133]
  1370. seq = Sequences(join(row['Base_Directory'], 'Preprocessed_Data/Sequences',
  1371. row['File']+'_MEC_immobility_sequences_cv.txt'),
  1372. load_objects=True)
  1374. seq.plot_sequence(indices=[1,33,41], figure_path=f"{figure_path}/Noether_DY02_MEC_alt_decoder_CV")
  1376. # %%
  1377. # decoder trained on movement variance = 6 instead of 1
  1378. sessions = pd.read_csv('Z:/WT_Sequences/all_sessions_Xmaze.csv')
  1379. immob_data_by_session_alt_decoder, _ = combine_metrics_over_sessions(sessions, animal_dict, 'MEC_immobility', modifier='_var6')
  1381. print(f"{np.nanmean(immob_data_by_session_alt_decoder['non_local']):.02f} +/- {stats.sem(np.nanmean(immob_data_by_session_alt_decoder['non_local'],axis=1)):.02f}" +\
  1382. f"% of immobility bouts have nonlocal coding,\n" +\
  1383. f"and {np.nanmean(immob_data_by_session_alt_decoder['frac_non_local']):.02f} +/- {stats.sem(np.nanmean(immob_data_by_session_alt_decoder['frac_non_local'],axis=1)):.02f}" +\
  1384. f"% of time during immobility "+\
  1385. "is spent representing nonlocal positions.")
  1387. plot_seq_metrics_over_sessions(n_animals, immob_data_by_session_alt_decoder['non_local'], figure_path,
  1388. 'immob_with_nl_alt_decoder_var6', '% Immobility Bouts with NonLocal Decode', [0,100])
  1390. plot_seq_metrics_over_sessions(n_animals, immob_data_by_session_alt_decoder['frac_non_local'], figure_path,
  1391. 'perc_immob_content_nl_alt_decoder_var6', '% Immobility Time with NonLocal Decode', [0,100])
  1393. # %%
  1394. # plot position and decoded position for snippet of data
  1395. # open example
  1396. row = sessions.iloc[133]
  1397. seq = Sequences(join(row['Base_Directory'], 'Preprocessed_Data/Sequences',
  1398. row['File']+'_MEC_immobility_sequences_var6.txt'),
  1399. load_objects=True)
  1401. seq.plot_sequence(indices=[1,33,41], figure_path=f"{figure_path}/Noether_DY02_MEC_alt_decoder_var6")
  1403. # %%
  1404. # Bayesian decoder
  1405. sessions = pd.read_csv('Z:/WT_Sequences/all_sessions_Xmaze.csv')
  1406. immob_data_by_session_alt_decoder, _ = combine_metrics_over_sessions(sessions, animal_dict, 'MEC_immobility', modifier='_bayes')
  1408. print(f"{np.nanmean(immob_data_by_session_alt_decoder['non_local']):.02f} +/- {stats.sem(np.nanmean(immob_data_by_session_alt_decoder['non_local'],axis=1)):.02f}" +\
  1409. f"% of immobility bouts have nonlocal coding,\n" +\
  1410. f"and {np.nanmean(immob_data_by_session_alt_decoder['frac_non_local']):.02f} +/- {stats.sem(np.nanmean(immob_data_by_session_alt_decoder['frac_non_local'],axis=1)):.02f}" +\
  1411. f"% of time during immobility "+\
  1412. "is spent representing nonlocal positions.")
  1414. plot_seq_metrics_over_sessions(n_animals, immob_data_by_session_alt_decoder['non_local'], figure_path,
  1415. 'immob_with_nl_alt_decoder_bayes', '% Immobility Bouts with NonLocal Decode', [0,100])
  1417. plot_seq_metrics_over_sessions(n_animals, immob_data_by_session_alt_decoder['frac_non_local'], figure_path,
  1418. 'perc_immob_content_nl_alt_decoder_bayes', '% Immobility Time with NonLocal Decode', [0,100])
  1420. # %%
  1421. # plot position and decoded position for snippet of data
  1422. # open example
  1423. row = sessions.iloc[133]
  1424. seq = Sequences(join(row['Base_Directory'], 'Preprocessed_Data/Sequences',
  1425. row['File']+'_MEC_immobility_sequences_bayes.txt'),
  1426. load_objects=True)
  1428. # subset data
  1429. start_time = 1260
  1430. end_time = 1290
  1431. start_idx = np.argmin(np.abs(seq.binned_data.timestamps-start_time))
  1432. end_idx = np.argmin(np.abs(seq.binned_data.timestamps-end_time))
  1433. subset_timestamps = seq.binned_data.timestamps[start_idx:end_idx]
  1435. f, ax = plt.subplots()
  1436. time_slice = slice(start_time,end_time)
  1438. # add classifier
  1439. results = seq.classifier.classifier_results.sel(time=time_slice)
  1440. time = (results.time-start_time) * 1000
  1441. max_time = time.max()
  1443. cmap = copy.copy(plt.cm.get_cmap('bone_r'))
  1444. cmap.set_bad(color="lightgrey", alpha=1.0)
  1445. (
  1446. results
  1447. .assign_coords(time=time)
  1448. .acausal_posterior.sum("state")
  1449. .plot(
  1450. x="time",
  1451. y="position",
  1452. robust=True,
  1453. add_colorbar=False,
  1454. zorder=0,
  1455. rasterized=True,
  1456. cmap=cmap,
  1457. ax=ax
  1458. )
  1459. )
  1460. seq_position = seq.binned_data.position['Linear'].loc[time_slice]
  1461. max_position = int(
  1462. np.ceil(seq.binned_data.position['Linear'].max()))
  1463. ax.plot(time, seq_position, linestyle="--", linewidth=2,
  1464. color="magenta", clip_on=False)
  1466. rtc.plot_graph_as_1D(seq.position.arena.track_graph,
  1467. edge_spacing=seq.position.arena.edge_spacing,
  1468. ax=ax, axis="y", other_axis_start=max_time+50)
  1470. f.savefig(join(figure_path,'Noether_DY02_decode_snippet_alt_decoder_bayes.pdf'), format='pdf')
  1472. # %%
  1473. # 2D decoder
  1474. sessions = pd.read_csv('Z:/WT_Sequences/all_sessions_Xmaze.csv')
  1475. immob_data_by_session_alt_decoder, _ = combine_metrics_over_sessions(sessions, animal_dict, 'MEC_immobility', modifier='_2D')
  1477. print(f"{np.nanmean(immob_data_by_session_alt_decoder['non_local']):.02f} +/- {stats.sem(np.nanmean(immob_data_by_session_alt_decoder['non_local'],axis=1)):.02f}" +\
  1478. f"% of immobility bouts have nonlocal coding,\n" +\
  1479. f"and {np.nanmean(immob_data_by_session_alt_decoder['frac_non_local']):.02f} +/- {stats.sem(np.nanmean(immob_data_by_session_alt_decoder['frac_non_local'],axis=1)):.02f}" +\
  1480. f"% of time during immobility "+\
  1481. "is spent representing nonlocal positions.")
  1483. plot_seq_metrics_over_sessions(n_animals, immob_data_by_session_alt_decoder['non_local'], figure_path,
  1484. 'immob_with_nl_alt_decoder_2D', '% Immobility Bouts with NonLocal Decode', [0,100])
  1486. plot_seq_metrics_over_sessions(n_animals, immob_data_by_session_alt_decoder['frac_non_local'], figure_path,
  1487. 'perc_immob_content_nl_alt_decoder_2D', '% Immobility Time with NonLocal Decode', [0,100])
  1489. # %%
  1490. # plot position and decoded position for snippet of data
  1491. # open example
  1492. row = sessions.iloc[133]
  1493. seq = Sequences(join(row['Base_Directory'], 'Preprocessed_Data/Sequences',
  1494. row['File']+'_MEC_immobility_sequences_2D.txt'),
  1495. load_objects=True)
  1497. # subset data
  1498. start_time = 1260
  1499. end_time = 1290
  1500. start_idx = np.argmin(np.abs(seq.binned_data.timestamps-start_time))
  1501. end_idx = np.argmin(np.abs(seq.binned_data.timestamps-end_time))
  1502. subset_timestamps = seq.binned_data.timestamps[start_idx:end_idx]
  1504. f, ax = plt.subplots()
  1505. time_slice = slice(start_time,end_time)
  1507. # add classifier
  1508. results = seq.classifier.classifier_results.sel(time=time_slice)
  1509. time = (results.time-start_time) * 1000
  1510. max_time = time.max()
  1512. cmap = copy.copy(plt.cm.get_cmap('bone_r'))
  1513. cmap.set_bad(color="lightgrey", alpha=1.0)
  1514. (
  1515. results
  1516. .assign_coords(time=time)
  1517. .acausal_posterior.sum("state")
  1518. .plot(
  1519. x="time",
  1520. y="position",
  1521. robust=True,
  1522. add_colorbar=False,
  1523. zorder=0,
  1524. rasterized=True,
  1525. cmap=cmap,
  1526. ax=ax
  1527. )
  1528. )
  1529. seq_position = seq.binned_data.position['Linear'].loc[time_slice]
  1530. max_position = int(
  1531. np.ceil(seq.binned_data.position['Linear'].max()))
  1532. ax.plot(time, seq_position, linestyle="--", linewidth=2,
  1533. color="magenta", clip_on=False)
  1535. rtc.plot_graph_as_1D(seq.position.arena.track_graph,
  1536. edge_spacing=seq.position.arena.edge_spacing,
  1537. ax=ax, axis="y", other_axis_start=max_time+50)
  1539. f.savefig(join(figure_path,'Noether_DY02_decode_snippet_alt_decoder_2D.pdf'), format='pdf')
  1541. # %% [markdown]
  1542. # # Figure 3 & S4
  1544. # %% [markdown]
  1545. # ### example cells with fields at local and non-local positions
  1547. # %%
  1548. # repeat raster with color coding by local/nonlocal
  1549. # first, load the interval
  1550. row = sessions.iloc[133]
  1551. seq = Sequences(join(row['Base_Directory'], 'Preprocessed_Data/Sequences',
  1552. row['File']+'_MEC_immobility_sequences.txt'),
  1553. load_objects=True)
  1554. print(seq.name)
  1556. start_time = seq.intervals.intervals.iloc[1].start_time
  1557. end_time = seq.intervals.intervals.iloc[1].end_time
  1558. start_idx = np.argmin(np.abs(seq.binned_data.timestamps-start_time))
  1559. end_idx = np.argmin(np.abs(seq.binned_data.timestamps-end_time))
  1560. subset_timestamps = seq.binned_data.timestamps[start_idx:end_idx]
  1562. # %%
  1563. # rasters: single cells
  1564. field_locs, field_peaks, field_sizes, field_spacing = seq.find_fields()
  1566. pf = seq.classifier.classifier.place_fields_[('', 0)]
  1567. n_cells = pf.shape[1]
  1568. decoded_bin = [59]
  1569. local_cells = []
  1570. for u in range(n_cells):
  1571. for f in field_locs[u]:
  1572. if np.any(np.isin(decoded_bin, f)):
  1573. local_cells.append(u)
  1574. local_cells = np.unique(np.asarray(local_cells))
  1576. decoded_bin = [29,92,93,94,95,96]
  1577. nonlocal_cells = []
  1578. for u in range(n_cells):
  1579. for f in field_locs[u]:
  1580. if np.any(np.isin(decoded_bin, f)):
  1581. nonlocal_cells.append(u)
  1582. nonlocal_cells = np.unique(np.asarray(nonlocal_cells))
  1584. # %%
  1585. for ex_cell in [20, 53, 30]:
  1586. interval_spikes = seq.binned_data.spikes[:,start_idx:end_idx]
  1587. f, ax = plt.subplots(1, figsize=(8,3))
  1588. # plot spikes from cells contributing to decode (over all time) in black
  1589. spike_time_ind = np.nonzero(interval_spikes[ex_cell, :])[0]
  1590. ax.scatter(1000*(subset_timestamps[spike_time_ind]-start_time),
  1591. np.ones(len(spike_time_ind)), color='black', zorder=1,
  1592. marker='|', linewidth=1)
  1594. # plot spikes contributing to nonlocal decode in cyan during non-local times
  1595. for i in range(len(nonlocal_intervals)):
  1596. interval_spikes = seq.binned_data.spikes[:,int(start_idx+nonlocal_intervals[i][0]*5):int(start_idx+nonlocal_intervals[i][1]*5)]
  1597. spike_time_ind_nl = np.nonzero(interval_spikes[ex_cell, :])[0]
  1598. ax.scatter(1000*(subset_timestamps[nonlocal_intervals[i][0]*5+spike_time_ind_nl]-start_time),
  1599. np.ones(len(spike_time_ind_nl)), color='cyan', zorder=1,
  1600. marker='|', linewidth=2)
  1601. ax.set(xlim=[0,int((end_time-start_time)*1000)], ylabel='Neuron')
  1602. f.savefig(join(figure_path,f'Noether_DY02_MEC_{ex_cell}_raster.pdf'), format='pdf')
  1604. # %%
  1605. bin_cm = 2
  1606. electrodes = Electrodes(join(row['Base_Directory'], 'Preprocessed_Data/Electrodes',
  1607. row['File']+'_imec1_electrodes.txt'))
  1609. MEC_channels = electrodes.subset_by_location(regions=['ENTm1','ENTm2','ENTm3','ENTm4','ENTm5','ENTm6'])
  1610. seq.spikes.subset_by_channel(MEC_channels)
  1612. for u in [20, 53, 30]:
  1613. fr_map, xbins, ybins, max_fr = calc_fr_map(seq.spikes.spikes[u], seq.position, bin_cm, smooth=True)
  1614. print(f'Unit {u} Max FR {max_fr}')
  1615. f = plot_heatmap(fr_map, xbins, ybins)
  1616. f.savefig(join(figure_path,f'Noether_DY02_MEC_{u}_heatmap.pdf'), format='pdf')
  1618. # %% [markdown]
  1619. # ### Cells that decode to non-local vs local positions during non-local vs local times
  1621. # %%
  1622. # participation of L vs NL coding cells during L vs NL intervals (L + NL multi-field cells?, L always active?)
  1624. #### for non-local times
  1625. nl_perc_cells, nl_perc_spikes, l_perc_cells, l_perc_spikes, lnl_perc_cells, lnl_perc_spikes, \
  1626. animal, session, dur = [[] for _ in range(9)]
  1628. for _, row in sessions.iterrows():
  1629. if not row['Recording_Error'] and not row['Position_Error'] and \
  1630. row['Task'] == 'X Maze' and int(row['Session'][-2:])<=10:
  1631. seq = Sequences(join(row['Base_Directory'], 'Preprocessed_Data/Sequences',
  1632. row['File']+'_MEC_nonlocal_immobility_sequences.txt'),
  1633. load_objects=True)
  1634. animal_idx = animal_dict[row['Animal']]
  1635. sess_idx = int(row['Session'][-2:])-1
  1636. animal.extend([animal_idx] * len(seq.stats))
  1637. session.extend([sess_idx] * len(seq.stats))
  1638. dur.extend(seq.stats.Duration)
  1640. linear_to_bins = np.array(seq.classifier.classifier.environments[0].place_bin_centers_nodes_df_['linear_position'])
  1641. pf = seq.classifier.classifier.place_fields_[('', 0)]
  1642. n_cells = pf.shape[1]
  1643. for s, interval in seq.intervals.intervals.iterrows():
  1644. cr = seq.classifier.classifier_results.sel(time=slice(interval.start_time, interval.end_time))
  1645. map_position_ind = cr.sum("state").acausal_posterior.argmax("position").values
  1646. decoded_bins = np.unique(map_position_ind)
  1647. linear_pos = np.asarray(seq.binned_data.position["Linear"].loc[interval.start_time:interval.end_time])
  1648. animal_bins = np.unique(np.argmin(np.abs(linear_to_bins[:,np.newaxis]-linear_pos), axis=0))
  1649. idx = seq.intervals.indices.iloc[s]
  1651. # find cells with significant fields at decoded locations
  1652. field_cells = []
  1653. for u in range(n_cells):
  1654. for f in field_locs[animal_idx][sess_idx][u]:
  1655. if np.any(np.isin(decoded_bins, f)):
  1656. field_cells.append(u)
  1657. field_cells = np.unique(np.asarray(field_cells, dtype=int))
  1659. # calculate proportion of all cells with fields,
  1660. # and of those: proportion active, fr, number of fields per cell, and size and spacing of those fields
  1661. seq_spikes = seq.binned_data.spikes[field_cells, idx['start_index']:idx['end_index']]
  1662. nl_perc_cells.append(np.sum(np.sum(seq_spikes, axis=1)>0)/seq.stats.Num_units[s])
  1663. nl_perc_spikes.append(np.sum(seq_spikes)/seq.stats.Num_spikes[s])
  1665. # repeat for cells with fields at current location
  1666. field_cells = []
  1667. for u in range(n_cells):
  1668. for f in field_locs[animal_idx][sess_idx][u]:
  1669. if np.any(np.isin(animal_bins, f)):
  1670. field_cells.append(u)
  1671. field_cells = np.unique(np.asarray(field_cells, dtype=int))
  1673. # calculate proportion of all cells with fields,
  1674. # and of those: proportion active, fr, number of fields per cell, and size and spacing of those fields
  1675. seq_spikes = seq.binned_data.spikes[field_cells, idx['start_index']:idx['end_index']]
  1676. l_perc_cells.append(np.sum(np.sum(seq_spikes, axis=1)>0)/seq.stats.Num_units[s])
  1677. l_perc_spikes.append(np.sum(seq_spikes)/seq.stats.Num_spikes[s])
  1679. # repeat for cells with fields at current AND decoded location
  1680. local_cells = field_cells
  1681. field_cells = []
  1682. for u in local_cells:
  1683. for f in field_locs[animal_idx][sess_idx][u]:
  1684. if np.any(np.isin(decoded_bins, f)):
  1685. field_cells.append(u)
  1686. field_cells = np.unique(np.asarray(field_cells, dtype=int))
  1688. # calculate proportion of all cells with fields,
  1689. # and of those: proportion active, fr, number of fields per cell, and size and spacing of those fields
  1690. seq_spikes = seq.binned_data.spikes[field_cells, idx['start_index']:idx['end_index']]
  1691. lnl_perc_cells.append(np.sum(np.sum(seq_spikes, axis=1)>0)/seq.stats.Num_units[s])
  1692. lnl_perc_spikes.append(np.sum(seq_spikes)/seq.stats.Num_spikes[s])
  1694. print(row['File'])
  1696. # append onto larger data_seq
  1697. nonlocal_fields_data_by_seq2 = pd.DataFrame({'Animal': animal, 'Session': session, 'Duration': dur,
  1698. 'Nonlocal_percent_cells': nl_perc_cells,
  1699. 'Nonlocal_percent_spikes': nl_perc_spikes,
  1700. 'Local_percent_cells': l_perc_cells,
  1701. 'Local_percent_spikes': l_perc_spikes,
  1702. 'Both_percent_cells': lnl_perc_cells,
  1703. 'Both_percent_spikes': lnl_perc_spikes})
  1705. # %%
  1706. #### for local times
  1707. nl_perc_cells, nl_perc_spikes, l_perc_cells, l_perc_spikes, lnl_perc_cells, lnl_perc_spikes, \
  1708. animal, session, dur = [[] for _ in range(9)]
  1710. for _, row in sessions.iterrows():
  1711. if not row['Recording_Error'] and not row['Position_Error'] and \
  1712. row['Task'] == 'X Maze' and int(row['Session'][-2:])<=10:
  1713. seq = Sequences(join(row['Base_Directory'], 'Preprocessed_Data/Sequences',
  1714. row['File']+'_MEC_local_immobility_sequences.txt'),
  1715. load_objects=True)
  1716. animal_idx = animal_dict[row['Animal']]
  1717. sess_idx = int(row['Session'][-2:])-1
  1718. animal.extend([animal_idx] * len(seq.stats))
  1719. session.extend([sess_idx] * len(seq.stats))
  1720. dur.extend(seq.stats.Duration)
  1722. linear_to_bins = np.array(seq.classifier.classifier.environments[0].place_bin_centers_nodes_df_['linear_position'])
  1723. pf = seq.classifier.classifier.place_fields_[('', 0)]
  1724. n_cells = pf.shape[1]
  1725. for s, interval in seq.intervals.intervals.iterrows():
  1726. cr = seq.classifier.classifier_results.sel(time=slice(interval.start_time, interval.end_time))
  1727. map_position_ind = cr.sum("state").acausal_posterior.argmax("position").values
  1728. decoded_bins = np.unique(map_position_ind)
  1729. linear_pos = np.asarray(seq.binned_data.position["Linear"].loc[interval.start_time:interval.end_time])
  1730. animal_bins = np.unique(np.argmin(np.abs(linear_to_bins[:,np.newaxis]-linear_pos), axis=0))
  1731. idx = seq.intervals.indices.iloc[s]
  1733. # find cells with significant fields at decoded locations
  1734. field_cells = []
  1735. for u in range(n_cells):
  1736. for f in field_locs[animal_idx][sess_idx][u]:
  1737. if np.any(np.isin(decoded_bins, f)):
  1738. field_cells.append(u)
  1739. field_cells = np.unique(np.asarray(field_cells, dtype=int))
  1741. # calculate proportion of all cells with fields,
  1742. # and of those: proportion active, fr, number of fields per cell, and size and spacing of those fields
  1743. seq_spikes = seq.binned_data.spikes[field_cells, idx['start_index']:idx['end_index']]
  1744. nl_perc_cells.append(np.sum(np.sum(seq_spikes, axis=1)>0)/seq.stats.Num_units[s])
  1745. nl_perc_spikes.append(np.sum(seq_spikes)/seq.stats.Num_spikes[s])
  1747. # repeat for cells with fields at current location
  1748. field_cells = []
  1749. for u in range(n_cells):
  1750. for f in field_locs[animal_idx][sess_idx][u]:
  1751. if np.any(np.isin(animal_bins, f)):
  1752. field_cells.append(u)
  1753. field_cells = np.unique(np.asarray(field_cells, dtype=int))
  1755. # calculate proportion of all cells with fields,
  1756. # and of those: proportion active, fr, number of fields per cell, and size and spacing of those fields
  1757. seq_spikes = seq.binned_data.spikes[field_cells, idx['start_index']:idx['end_index']]
  1758. l_perc_cells.append(np.sum(np.sum(seq_spikes, axis=1)>0)/seq.stats.Num_units[s])
  1759. l_perc_spikes.append(np.sum(seq_spikes)/seq.stats.Num_spikes[s])
  1761. # repeat for cells with fields at current AND decoded location
  1762. local_cells = field_cells
  1763. field_cells = []
  1764. for u in local_cells:
  1765. for f in field_locs[animal_idx][sess_idx][u]:
  1766. if np.any(np.isin(decoded_bins, f)):
  1767. field_cells.append(u)
  1768. field_cells = np.unique(np.asarray(field_cells, dtype=int))
  1770. # calculate proportion of all cells with fields,
  1771. # and of those: proportion active, fr, number of fields per cell, and size and spacing of those fields
  1772. seq_spikes = seq.binned_data.spikes[field_cells, idx['start_index']:idx['end_index']]
  1773. lnl_perc_cells.append(np.sum(np.sum(seq_spikes, axis=1)>0)/seq.stats.Num_units[s])
  1774. lnl_perc_spikes.append(np.sum(seq_spikes)/seq.stats.Num_spikes[s])
  1776. print(row['File'])
  1778. # append onto larger data_seq
  1779. local_fields_data_by_seq2 = pd.DataFrame({'Animal': animal, 'Session': session, 'Duration': dur,
  1780. 'Nonlocal_percent_cells': nl_perc_cells,
  1781. 'Nonlocal_percent_spikes': nl_perc_spikes,
  1782. 'Local_percent_cells': l_perc_cells,
  1783. 'Local_percent_spikes': l_perc_spikes,
  1784. 'Both_percent_cells': lnl_perc_cells,
  1785. 'Both_percent_spikes': lnl_perc_spikes})
  1787. # %%
  1788. with open(join(figure_path,'nonlocal_fields_data_by_seq2.pkl'), 'wb') as file:
  1789. pickle.dump(nonlocal_fields_data_by_seq2, file)
  1790. with open(join(figure_path,'local_fields_data_by_seq2.pkl'), 'wb') as file:
  1791. pickle.dump(local_fields_data_by_seq2, file)
  1793. # %%
  1794. plot_compare_paired_metrics('Local_percent_cells', 'Nonlocal_percent_cells', nonlocal_fields_data_by_seq2, figure_path,
  1795. 'nonlocal_l_vs_nl_field_cells', '% Active Cells That Have Fields at Location', [0,1])
  1796. plot_compare_paired_metrics('Nonlocal_percent_cells', 'Both_percent_cells', nonlocal_fields_data_by_seq2, figure_path,
  1797. 'nonlocal_both_vs_nl_field_cells', '% Active Cells That Have Fields at Location', [0,1])
  1798. plot_compare_paired_metrics('Local_percent_spikes', 'Nonlocal_percent_spikes', nonlocal_fields_data_by_seq2, figure_path,
  1799. 'nonlocal_l_vs_nl_field_spikes', '% Spikes from Cells That Have Fields at Location', [0,1])
  1800. plot_compare_paired_metrics('Nonlocal_percent_spikes', 'Both_percent_spikes', nonlocal_fields_data_by_seq2, figure_path,
  1801. 'nonlocal_both_vs_nl_field_spikes', '% Spikes from Cells That Have Fields at Location', [0,1])
  1803. # plot_compare_paired_metrics('Local_percent_cells', 'Nonlocal_percent_cells', local_fields_data_by_seq2, figure_path,
  1804. # 'local_l_vs_nl_field_cells', '% Active Cells That Have Fields at Location', [0,1])
  1805. # plot_compare_paired_metrics('Nonlocal_percent_cells', 'Both_percent_cells', local_fields_data_by_seq2, figure_path,
  1806. # 'local_both_vs_nl_field_cells', '% Active Cells That Have Fields at Location', [0,1])
  1807. # plot_compare_paired_metrics('Local_percent_spikes', 'Nonlocal_percent_spikes', local_fields_data_by_seq2, figure_path,
  1808. # 'local_l_vs_nl_field_spikes', '% Spikes from Cells That Have Fields at Location', [0,1])
  1809. # plot_compare_paired_metrics('Nonlocal_percent_spikes', 'Both_percent_spikes', local_fields_data_by_seq2, figure_path,
  1810. # 'local_both_vs_nl_field_spikes', '% Spikes from Cells That Have Fields at Location', [0,1])
  1812. # %%
  1813. is_non_local = []
  1814. is_non_local.extend([False] * len(local_fields_data_by_seq2))
  1815. is_non_local.extend([True] * len(nonlocal_fields_data_by_seq2))
  1816. local_nonlocal_fields_equivalent_data_by_seq2 = pd.concat([local_fields_data_by_seq2, nonlocal_fields_data_by_seq2], ignore_index=True)
  1817. local_nonlocal_fields_equivalent_data_by_seq2.loc[:,'Non_local'] = is_non_local
  1818. plot_seq_metrics_comparison_over_events(local_nonlocal_fields_equivalent_data_by_seq2, 'Local_percent_cells', 'Non_local', figure_path,
  1819. 'local_vs_nonlocal_l_perc_cells', '% Active Cells That Have Fields at Location', [0,1])
  1820. plot_seq_metrics_comparison_over_events(local_nonlocal_fields_equivalent_data_by_seq2, 'Local_percent_spikes', 'Non_local', figure_path,
  1821. 'local_vs_nonlocal_l_perc_spikes', '% Spikes from Cells That Have Fields at Location', [0,1])
  1823. # %% [markdown]
  1824. # ### are decodes less non-local without non-local spiking
  1826. # %%
  1827. # remove these spikes from test set to show now decodes are only local
  1828. for _, row in sessions.iterrows():
  1829. if not row['Recording_Error'] and not row['Position_Error'] and \
  1830. row['Task'] == 'X Maze' and int(row['Session'][-2:])<=10:
  1831. seq = Sequences(join(row['Base_Directory'], 'Preprocessed_Data/Sequences',
  1832. row['File']+'_MEC_nonlocal_immobility_sequences.txt'),
  1833. load_objects=True)
  1834. animal_idx = animal_dict[row['Animal']]
  1835. sess_idx = int(row['Session'][-2:])-1
  1837. linear_to_bins = np.array(seq.classifier.classifier.environments[0].place_bin_centers_nodes_df_['linear_position'])
  1838. pf = seq.classifier.classifier.place_fields_[('', 0)]
  1839. n_cells = pf.shape[1]
  1840. for s, interval in seq.intervals.intervals.iterrows():
  1841. cr = seq.classifier.classifier_results.sel(time=slice(interval.start_time, interval.end_time))
  1842. map_position_ind = cr.sum("state").acausal_posterior.argmax("position").values
  1843. decoded_bins = np.unique(map_position_ind)
  1844. idx = seq.intervals.indices.iloc[s]
  1846. # find cells with significant fields at decoded locations
  1847. field_cells = []
  1848. for u in range(n_cells):
  1849. for f in field_locs[animal_idx][sess_idx][u]:
  1850. if np.any(np.isin(decoded_bins, f)):
  1851. field_cells.append(u)
  1852. field_cells = np.unique(np.asarray(field_cells, dtype=int))
  1854. # remove non-local decode spikes
  1855. seq.binned_data.spikes[field_cells, idx['start_index']:idx['end_index']] = 0
  1857. # save binned data
  1858. binned_data = Binned_Data(seq.binned_data.name)
  1859. binned_data.spikes = seq.binned_data.spikes
  1860. binned_data.name = seq.binned_data.name[:-4]+'_NLspikesremoved.txt'
  1861. binned_data.save()
  1863. print(row['File'])
  1865. # %%
  1866. sessions = pd.read_csv('Z:/WT_Sequences/all_sessions_Xmaze.csv')
  1867. immob_data_by_session_CLtest, immob_data_by_seq_CLtest = combine_metrics_over_sessions(sessions, animal_dict, 'MEC_immobility_NLspikesremoved')
  1869. # %%
  1870. print(f"{np.nanmean(immob_data_by_session_CLtest['non_local']):.02f} +/- {stats.sem(np.nanmean(immob_data_by_session_CLtest['non_local'],axis=1)):.02f}" +\
  1871. f"% of immobility bouts have nonlocal coding,\n" +\
  1872. f"and {np.nanmean(immob_data_by_session_CLtest['frac_non_local']):.02f} +/- {stats.sem(np.nanmean(immob_data_by_session_CLtest['frac_non_local'],axis=1)):.02f}" +\
  1873. f"% of time during immobility "+\
  1874. "is spent representing nonlocal positions.")
  1876. # %%
  1877. plot_seq_metrics_over_sessions(n_animals, immob_data_by_session_CLtest['non_local'], figure_path,
  1878. 'immob_with_nl_CLtest', '% Immobility Bouts with NonLocal Decode', [0,100])
  1880. plot_seq_metrics_over_sessions(n_animals, immob_data_by_session_CLtest['frac_non_local'], figure_path,
  1881. 'perc_immob_content_nl_CLtest', '% Immobility Time with NonLocal Decode', [0,100])
  1883. # %% [markdown]
  1884. # ### Preferentially enriched cells
  1886. # %%
  1887. seq_type = 'MEC_nonlocal_immobility'
  1888. sessions = pd.read_csv(r'Z:/WT_Sequences/all_sessions_Xmaze.csv')
  1889. n_shuffles = 1000
  1890. n_sessions = 20
  1891. n_animals = len(animal_dict)
  1892. over_fr = np.full((n_animals, n_sessions), np.NaN)
  1893. over_fr_list = [[[] for _ in range(n_sessions)] for _ in range(n_animals)]
  1894. over_percent_intervals = np.full((n_animals, n_sessions), np.NaN)
  1895. over_percent_intervals_list = [[[] for _ in range(n_sessions)] for _ in range(n_animals)]
  1897. for _, row in sessions.iterrows():
  1898. if not row['Recording_Error'] and not row['Position_Error'] and \
  1899. row['Task'] == 'X Maze' and int(row['Session'][-2:])<=10:
  1900. # load spikes and intervals
  1901. binned_data = Binned_Data(join(row['Base_Directory'],
  1902. 'Preprocessed_Data/Sequences/Binned_Data',
  1903. row['File']+'_MEC_binned_data.txt'))
  1904. intervals = Sequence_Intervals(join(row['Base_Directory'],
  1905. 'Preprocessed_Data/Sequences/Intervals',
  1906. row['File']+'_'+seq_type+'_intervals.txt'))
  1907. seq = Sequences(join(row['Base_Directory'],
  1908. 'Preprocessed_Data/Sequences',
  1909. row['File']+'_'+seq_type+'_sequences.txt'))
  1910. interval_dur = intervals.intervals.end_time-intervals.intervals.start_time
  1911. animal_idx = animal_dict[row['Animal']]
  1912. sess_idx = int(row['Session'][-2:])-1
  1914. shuffled_fr = np.zeros((n_shuffles, binned_data.spikes.shape[0]))
  1915. shuffled_percent_intervals = np.zeros((n_shuffles, binned_data.spikes.shape[0]))
  1917. for s in range(n_shuffles):
  1918. # shuffle each cell separately
  1919. shuffled_spikes = binned_data.spikes[:, np.random.permutation(binned_data.spikes.shape[1])]
  1921. # calculate fr during intervals & participation in intervals
  1922. # (same code as calc_cell_metrics in Sequences class)
  1923. fr = []
  1924. percent_intervals = []
  1925. for c in range(binned_data.spikes.shape[0]):
  1926. fr_by_interval = 0
  1927. participating_seq = 0
  1928. for i, idx in intervals.indices.iterrows():
  1929. num_spikes = np.sum(shuffled_spikes[c, idx['start_index']:idx['end_index']])
  1930. fr_by_interval += num_spikes/interval_dur.iloc[i]
  1931. if num_spikes>0:
  1932. participating_seq += 1
  1933. fr.append(fr_by_interval/intervals.n_seq)
  1934. percent_intervals.append(participating_seq/intervals.n_seq)
  1936. # append to shuffle distribution
  1937. shuffled_fr[s,:] = fr
  1938. shuffled_percent_intervals[s,:] = percent_intervals
  1940. # calculate 95th percentile
  1941. fr_95_percentile = np.percentile(shuffled_fr, 95, axis=0)
  1942. percent_intervals_95_percentile = np.percentile(shuffled_percent_intervals, 95, axis=0)
  1944. # count how many cells are above this threshold
  1945. over_fr[animal_idx, sess_idx] = np.nanmean(seq.stats_by_cell.FR>fr_95_percentile)
  1946. over_fr_list[animal_idx][sess_idx] = np.where(seq.stats_by_cell.FR>fr_95_percentile)[0].tolist()
  1947. over_percent_intervals[animal_idx, sess_idx] = \
  1948. np.nanmean(seq.stats_by_cell.Percent_intervals>percent_intervals_95_percentile)
  1949. over_percent_intervals_list[animal_idx][sess_idx] = \
  1950. np.where(seq.stats_by_cell.Percent_intervals>percent_intervals_95_percentile)[0].tolist()
  1951. print(row['File'])
  1953. # %%
  1954. with open(join(figure_path,'over_fr.pkl'), 'wb') as file:
  1955. pickle.dump(over_fr, file)
  1956. with open(join(figure_path,'over_percent_intervals.pkl'), 'wb') as file:
  1957. pickle.dump(over_percent_intervals, file)
  1958. with open(join(figure_path,'over_fr_list.pkl'), 'wb') as file:
  1959. pickle.dump(over_fr_list, file)
  1960. with open(join(figure_path,'over_percent_intervals_list.pkl'), 'wb') as file:
  1961. pickle.dump(over_percent_intervals_list, file)
  1963. # %%
  1964. print(f"{np.nanmean(over_fr)*100:.02f} +/- {stats.sem(np.nanmean(over_fr,axis=1))*100:.02f}" +\
  1965. f"% of cells have > 95th percentile of shuffled FR during nonlocal coding.")
  1967. # %%
  1968. # plot
  1969. plot_seq_metrics_over_sessions(n_animals, over_fr*100, figure_path,
  1970. 'recruited_fr', '% Units with FR > 95th Percentile', [0,100])
  1972. # %% [markdown]
  1973. # ### DV axis
  1975. # %%
  1976. sessions = pd.read_csv(r'Z:/WT_Sequences/all_sessions_Xmaze.csv')
  1977. n_sessions = 10
  1978. n_animals = len(animal_dict)
  1979. dist_from_border = [[[] for _ in range(n_sessions)] for _ in range(n_animals)]
  1981. for _, row in sessions.iterrows():
  1982. if not row['Recording_Error'] and not row['Position_Error'] and \
  1983. row['Task'] == 'X Maze' and int(row['Session'][-2:])<=10:
  1984. animal_idx = animal_dict[row['Animal']]
  1985. sess_idx = int(row['Session'][-2:])-1
  1986. seq = Sequences(join(row['Base_Directory'],
  1987. 'Preprocessed_Data/Sequences',
  1988. row['File']+'_MEC_nonlocal_immobility_sequences.txt'))
  1990. if ('2021_pilot' in row['Base_Directory']) or ('2022_winter' in row['Base_Directory']):
  1991. electrodes = Electrodes(join(row['Base_Directory'], 'Preprocessed_Data/Electrodes',
  1992. row['File']+'_electrodes.txt'))
  1993. channels = electrodes.subset_by_location(regions=['Entorhinal area medial part dorsal zone'])
  1994. spikes = Spikes(join(row['Base_Directory'],
  1995. 'Preprocessed_Data/Spikes/g1',
  1996. row['File']+'_imec0_spikes.txt'))
  1997. spikes.subset_by_channel(channels)
  1998. dist_from_border[animal_idx][sess_idx] = [depth - border_dict[row['Animal']] for depth in electrodes.electrodes.Rel_Z.iloc[spikes.spike_channel].values.tolist()]
  1999. else:
  2000. if (row['Animal']=='Lamarr') and (sess_idx>=7):
  2001. electrodes = Electrodes(join(row['Base_Directory'], 'Preprocessed_Data/Electrodes',
  2002. row['File']+'_imec0_electrodes.txt'))
  2003. spikes = Spikes(join(row['Base_Directory'], 'Preprocessed_Data/Spikes',
  2004. row['File']+'_imec0_spikes.txt'))
  2005. else:
  2006. electrodes = Electrodes(join(row['Base_Directory'], 'Preprocessed_Data/Electrodes',
  2007. row['File']+'_imec1_electrodes.txt'))
  2008. spikes = Spikes(join(row['Base_Directory'], 'Preprocessed_Data/Spikes',
  2009. row['File']+'_imec1_spikes.txt'))
  2010. channels = electrodes.subset_by_location(regions=['ENTm1','ENTm2','ENTm3','ENTm4','ENTm5','ENTm6'])
  2011. spikes.subset_by_channel(channels)
  2012. shanks = (electrodes.electrodes.Rel_Y.iloc[spikes.spike_channel].values/250).astype(int)
  2013. depths = electrodes.electrodes.Rel_Z.iloc[spikes.spike_channel].values
  2014. for i, d in enumerate(depths):
  2015. dist_from_border[animal_idx][sess_idx].append(border_dict[row['Animal']][shanks[i]] - d)
  2017. # %%
  2018. with open(join(figure_path,'dist_from_border.pkl'), 'wb') as file:
  2019. pickle.dump(dist_from_border, file)
  2021. # %%
  2022. # bin cells by DV depth
  2023. n_sessions = 10
  2024. dv_bin1, dv_bin2, dv_bin3, dv_bin4 = [[[[] for _ in range(n_sessions)] for _ in range(n_animals)] for _ in range(4)]
  2025. for animal_idx in range(n_animals):
  2026. for sess_idx in range(n_sessions):
  2027. dv_bin1[animal_idx][sess_idx] = np.where([i<500 for i in dist_from_border[animal_idx][sess_idx]])[0]
  2028. dv_bin2[animal_idx][sess_idx] = np.where([i>=500 and i<1000 for i in dist_from_border[animal_idx][sess_idx]])[0]
  2029. dv_bin3[animal_idx][sess_idx] = np.where([i>=1000 and i<1500 for i in dist_from_border[animal_idx][sess_idx]])[0]
  2030. dv_bin4[animal_idx][sess_idx] = np.where([i>=1500 for i in dist_from_border[animal_idx][sess_idx]])[0]
  2032. # %%
  2033. # compare DV depths of enriched vs unenriched cells
  2034. sess_idx = 0
  2035. percentile = over_fr_list
  2037. enriched_cells = []
  2038. enriched_animals = []
  2039. unenriched_cells = []
  2040. unenriched_animals = []
  2041. for animal_idx in range(len(animal_dict)):
  2042. unenriched_cells.extend([dist_from_border[animal_idx][sess_idx][i] for i in range(len(dist_from_border[animal_idx][sess_idx])) \
  2043. if i not in percentile[animal_idx][sess_idx]])
  2044. unenriched_animals.extend([animal_idx]*(len(dist_from_border[animal_idx][sess_idx]) - len(percentile[animal_idx][sess_idx])))
  2045. enriched_cells.extend([dist_from_border[animal_idx][sess_idx][i] for i in percentile[animal_idx][sess_idx]])
  2046. enriched_animals.extend([animal_idx]*len(percentile[animal_idx][sess_idx]))
  2048. data_df = pd.DataFrame({'Animal': unenriched_animals + enriched_animals,
  2049. 'Session_Unique': [sess_idx]*len(unenriched_cells + enriched_cells),
  2050. 'Condition': [0]*len(unenriched_cells) + [1]*len(enriched_cells),
  2051. 'Indep_Var': unenriched_cells + enriched_cells})
  2052. data_df.loc[data_df.Indep_Var<0,'Indep_Var'] = 0
  2054. plot_compare_unpaired_metrics(unenriched_cells, enriched_cells, data_df, figure_path, \
  2055. 'nonlocal_enriched_cells_by_fr', 'Recruited > 95th percentile', [0, 3000])
  2057. # %%
  2058. # enrichment: what proportion of cells from each depth bin are on the enrichment lists?
  2059. enriched_cells = over_fr_list
  2061. enriched_dv_bin1, enriched_dv_bin2, enriched_dv_bin3, enriched_dv_bin4 = \
  2062. [np.full((len(animal_dict),n_sessions), np.NaN) for _ in range(4)]
  2063. for animal_idx in range(n_animals):
  2064. for sess_idx in range(n_sessions):
  2065. # proportion of cells in enrichment list
  2066. enriched_dv_bin1[animal_idx][sess_idx] = len(np.intersect1d(dv_bin1[animal_idx][sess_idx], enriched_cells[animal_idx][sess_idx]))/ \
  2067. len(dv_bin1[animal_idx][sess_idx]) if len(dv_bin1[animal_idx][sess_idx])>0 else np.NaN
  2068. enriched_dv_bin2[animal_idx][sess_idx] = len(np.intersect1d(dv_bin2[animal_idx][sess_idx], enriched_cells[animal_idx][sess_idx]))/ \
  2069. len(dv_bin2[animal_idx][sess_idx]) if len(dv_bin2[animal_idx][sess_idx])>0 else np.NaN
  2070. enriched_dv_bin3[animal_idx][sess_idx] = len(np.intersect1d(dv_bin3[animal_idx][sess_idx], enriched_cells[animal_idx][sess_idx]))/ \
  2071. len(dv_bin3[animal_idx][sess_idx]) if len(dv_bin3[animal_idx][sess_idx])>0 else np.NaN
  2072. enriched_dv_bin4[animal_idx][sess_idx] = len(np.intersect1d(dv_bin4[animal_idx][sess_idx], enriched_cells[animal_idx][sess_idx]))/ \
  2073. len(dv_bin4[animal_idx][sess_idx]) if len(dv_bin4[animal_idx][sess_idx])>0 else np.NaN
  2075. plot_feature_by_depth_bin(enriched_dv_bin1, enriched_dv_bin2, enriched_dv_bin3, enriched_dv_bin4, \
  2076. figure_path, 'proportion_enriched_by_depth_bin_fr', '% Units with Recruitment > 95th Percentile', [0,1])
  2078. # %%
  2079. # calculate field sizes and spacing
  2080. sessions = pd.read_csv(r'Z:/WT_Sequences/all_sessions_Xmaze.csv')
  2081. n_sessions = 10
  2082. n_animals = len(animal_dict)
  2084. field_locs = np.empty((n_animals, n_sessions), dtype=object)
  2085. field_peaks = np.empty((n_animals, n_sessions), dtype=object)
  2086. field_sizes = np.empty((n_animals, n_sessions), dtype=object)
  2087. field_spacing = np.empty((n_animals, n_sessions), dtype=object)
  2088. for _, row in sessions.iterrows():
  2089. if not row['Recording_Error'] and not row['Position_Error'] and \
  2090. row['Task'] == 'X Maze' and int(row['Session'][-2:])<=10:
  2091. # any sequence will do; we're just using its contained objects
  2092. seq = Sequences(join(row['Base_Directory'], 'Preprocessed_Data/Sequences',
  2093. row['File']+'_MEC_nonlocal_immobility_sequences.txt'),
  2094. load_objects=True)
  2095. animal_idx = animal_dict[row['Animal']]
  2096. sess_idx = int(row['Session'][-2:])-1
  2097. try:
  2098. field_locs[animal_idx, sess_idx], field_peaks[animal_idx, sess_idx], \
  2099. field_sizes[animal_idx, sess_idx], field_spacing[animal_idx, sess_idx] = seq.find_fields()
  2100. print(row['File'])
  2101. except:
  2102. print(f"ERROR: {row['File']}")
  2104. # %%
  2105. with open(join(figure_path,'field_locs.pkl'), 'wb') as file:
  2106. pickle.dump(field_locs, file)
  2107. with open(join(figure_path,'field_peaks.pkl'), 'wb') as file:
  2108. pickle.dump(field_peaks, file)
  2109. with open(join(figure_path,'field_sizes.pkl'), 'wb') as file:
  2110. pickle.dump(field_sizes, file)
  2111. with open(join(figure_path,'field_spacing.pkl'), 'wb') as file:
  2112. pickle.dump(field_spacing, file)
  2114. # %%
  2115. # is this due to enriched cells having higher field sizes?
  2116. # repeat with FR
  2117. percentile = over_fr_list
  2118. sess_idx = 0
  2120. # suppress the "Mean of an empty slice" warning that will occur thousands of times
  2121. # (due to cells with single fields having empty arrays for field_spacing)
  2122. with warnings.catch_warnings():
  2123. warnings.filterwarnings("ignore", category=RuntimeWarning)
  2125. enriched_cells = []
  2126. enriched_animals = []
  2127. unenriched_cells = []
  2128. unenriched_animals = []
  2129. for animal_idx in range(len(animal_dict)):
  2130. if field_sizes[animal_idx][sess_idx] is not None:
  2131. unenriched_cells.extend([np.nanmean(field_sizes[animal_idx][sess_idx][i]) for i in range(len(field_sizes[animal_idx][sess_idx])) \
  2132. if i not in percentile[animal_idx][sess_idx]])
  2133. unenriched_animals.extend([animal_idx]*(len(field_sizes[animal_idx][sess_idx]) - len(percentile[animal_idx][sess_idx])))
  2134. enriched_cells.extend([np.nanmean(field_sizes[animal_idx][sess_idx][i]) for i in percentile[animal_idx][sess_idx]])
  2135. enriched_animals.extend([animal_idx]*len(percentile[animal_idx][sess_idx]))
  2137. data_df = pd.DataFrame({'Animal': unenriched_animals + enriched_animals,
  2138. 'Session_Unique': [sess_idx]*len(unenriched_cells + enriched_cells),
  2139. 'Condition': [0]*len(unenriched_cells) + [1]*len(enriched_cells),
  2140. 'Indep_Var': unenriched_cells + enriched_cells})
  2141. data_df.loc[data_df.Indep_Var<0,'Indep_Var'] = 0
  2143. plot_compare_unpaired_metrics(unenriched_cells, enriched_cells, data_df, figure_path, \
  2144. 'field_sizes_nonlocal_enriched_cells_by_fr', 'Recruited > 95th percentile', [0, 100])
  2146. # %%
  2147. # bin cells by field size
  2148. n_sessions = 10
  2149. fs_bin1, fs_bin2, fs_bin3, fs_bin4 = [[[[] for _ in range(n_sessions)] for _ in range(n_animals)] for _ in range(4)]
  2151. # suppress the "Mean of an empty slice" warning that will occur thousands of times
  2152. with warnings.catch_warnings():
  2153. warnings.filterwarnings("ignore", category=RuntimeWarning)
  2154. for animal_idx in range(n_animals):
  2155. for sess_idx in range(n_sessions):
  2156. if field_sizes[animal_idx][sess_idx] is not None:
  2157. fs = [np.nanmean(field_sizes[animal_idx][sess_idx][i]) for i in range(len(field_sizes[animal_idx][sess_idx]))]
  2158. fs_bin1[animal_idx][sess_idx] = np.where([i<22 for i in fs])[0]
  2159. fs_bin2[animal_idx][sess_idx] = np.where([i>=22 and i<28 for i in fs])[0]
  2160. fs_bin3[animal_idx][sess_idx] = np.where([i>=28 and i<34 for i in fs])[0]
  2161. fs_bin4[animal_idx][sess_idx] = np.where([i>=34 for i in fs])[0]
  2163. # %%
  2164. # what proportion of cells from each field size bin are on the enrichment lists?
  2165. enriched_cells = over_fr_list
  2167. enriched_fs_bin1, enriched_fs_bin2, enriched_fs_bin3, enriched_fs_bin4 = \
  2168. [np.full((len(animal_dict),n_sessions), np.NaN) for _ in range(4)]
  2169. for animal_idx in range(n_animals):
  2170. for sess_idx in range(n_sessions):
  2171. # proportion of cells in enrichment list
  2172. enriched_fs_bin1[animal_idx][sess_idx] = len(np.intersect1d(fs_bin1[animal_idx][sess_idx], enriched_cells[animal_idx][sess_idx]))/ \
  2173. len(fs_bin1[animal_idx][sess_idx]) if len(fs_bin1[animal_idx][sess_idx])>0 else np.NaN
  2174. enriched_fs_bin2[animal_idx][sess_idx] = len(np.intersect1d(fs_bin2[animal_idx][sess_idx], enriched_cells[animal_idx][sess_idx]))/ \
  2175. len(fs_bin2[animal_idx][sess_idx]) if len(fs_bin2[animal_idx][sess_idx])>0 else np.NaN
  2176. enriched_fs_bin3[animal_idx][sess_idx] = len(np.intersect1d(fs_bin3[animal_idx][sess_idx], enriched_cells[animal_idx][sess_idx]))/ \
  2177. len(fs_bin3[animal_idx][sess_idx]) if len(fs_bin3[animal_idx][sess_idx])>0 else np.NaN
  2178. enriched_fs_bin4[animal_idx][sess_idx] = len(np.intersect1d(fs_bin4[animal_idx][sess_idx], enriched_cells[animal_idx][sess_idx]))/ \
  2179. len(fs_bin4[animal_idx][sess_idx]) if len(fs_bin4[animal_idx][sess_idx])>0 else np.NaN
  2181. plot_feature_by_depth_bin(enriched_fs_bin1, enriched_fs_bin2, enriched_fs_bin3, enriched_fs_bin4, \
  2182. figure_path, 'proportion_enriched_by_field_size_fr', '% Units with Recruitment > 95th Percentile', [0,1])
  2184. # %% [markdown]
  2185. # ### Spatial variables
  2187. # %%
  2188. # calculate enrichment during Xmaze non-local coding using Xmaze+OF-sorted sessions
  2189. start_time = time.time()
  2191. sessions = pd.read_csv(r'Z:/WT_Sequences/all_sessions_Xmaze_celltypes_g1.csv')
  2192. n_shuffles = 1000
  2193. n_sessions = 10
  2194. n_animals = len(animal_dict)
  2195. over_fr = np.full((n_animals, n_sessions), np.NaN)
  2196. over_fr_list = [[[] for _ in range(n_sessions)] for _ in range(n_animals)]
  2198. for _, row in sessions.iterrows():
  2199. if not row['Recording_Error'] and not row['Position_Error'] and \
  2200. row['Task'] == 'X Maze' and int(row['Session'][-2:])<=10:
  2201. # load spikes and intervals
  2202. binned_data = Binned_Data(join(row['Base_Directory'],
  2203. 'Preprocessed_Data/Sequences/Binned_Data',
  2204. row['File']+'_MEC_g1+3_binned_data.txt'))
  2205. intervals = Sequence_Intervals(join(row['Base_Directory'],
  2206. 'Preprocessed_Data/Sequences/Intervals',
  2207. row['File']+'_MEC_nonlocal_immobility_intervals.txt'))
  2208. seq = Sequences(join(row['Base_Directory'],
  2209. 'Preprocessed_Data/Sequences',
  2210. row['File']+'_MEC_g1+3_nonlocal_immobility_sequences.txt'))
  2211. interval_dur = intervals.intervals.end_time-intervals.intervals.start_time
  2212. animal_idx = animal_dict[row['Animal']]
  2213. sess_idx = int(row['Session'][-2:])-1
  2215. shuffled_fr = np.zeros((n_shuffles, binned_data.spikes.shape[0]))
  2217. for s in range(n_shuffles):
  2218. # shuffle each cell separately
  2219. shuffled_spikes = binned_data.spikes[:, np.random.permutation(binned_data.spikes.shape[1])]
  2221. # calculate fr during intervals & participation in intervals
  2222. # (same code as calc_cell_metrics in Sequences class)
  2223. fr = []
  2224. for c in range(binned_data.spikes.shape[0]):
  2225. fr_by_interval = 0
  2226. for i, idx in intervals.indices.iterrows():
  2227. num_spikes = np.sum(shuffled_spikes[c, idx['start_index']:idx['end_index']])
  2228. fr_by_interval += num_spikes/interval_dur.iloc[i]
  2229. fr.append(fr_by_interval/intervals.n_seq)
  2231. # append to shuffle distribution
  2232. shuffled_fr[s,:] = fr
  2234. # calculate 95th percentile
  2235. fr_95_percentile = np.percentile(shuffled_fr, 95, axis=0)
  2237. # count how many cells are above this threshold
  2238. over_fr[animal_idx, sess_idx] = np.nanmean(seq.stats_by_cell.FR>fr_95_percentile)
  2239. over_fr_list[animal_idx][sess_idx] = np.where(seq.stats_by_cell.FR>fr_95_percentile)[0].tolist()
  2240. print(f"{row['File']} {time.time() - start_time} seconds")
  2242. # %%
  2243. with open(join(figure_path,'over_fr_OF.pkl'), 'wb') as file:
  2244. pickle.dump(over_fr, file)
  2245. with open(join(figure_path,'over_fr_list_OF.pkl'), 'wb') as file:
  2246. pickle.dump(over_fr_list, file)
  2248. # %%
  2249. # refs: border (Solstad 2008), speed (Kropff 2015), HD (Taube 1990), grid (Hafting 2005), spatial aperiodic (Diehl 2017)
  2250. start_time = time.time()
  2252. bin_cm = 2.5
  2254. sessions = pd.read_csv('Z:/WT_Sequences/all_sessions_Xmaze_celltypes.csv')
  2255. n_sessions = 10
  2256. n_animals = len(animal_dict)
  2257. spatial_cells, spatial_aperiodic_cells, grid_cells, border_cells, speed_cells, hd_cells = \
  2258. [[[[] for _ in range(n_sessions)] for _ in range(n_animals)] for _ in range(6)]
  2260. for i, row in sessions.iterrows():
  2261. if not row['Recording_Error'] and not row['Position_Error']:
  2263. animal_idx = animal_dict[row['Animal']]
  2264. sess_idx = int(row['Session'][-2:])-1
  2265. position = Position(name=join(row['Base_Directory'], 'Preprocessed_Data/Position',
  2266. row['File']+'_position.txt'))
  2267. position.position = position.position[position.position.Velocity>2]
  2269. if ("2023_spring" in row['Base_Directory']) or ("2024_winter" in row['Base_Directory']):
  2270. if (row['Animal']=='Lamarr') and (int(row['Session'][-2:])>=8):
  2271. mec_spikes = Spikes(name=join(row['Base_Directory'], 'Preprocessed_Data/Spikes/g1_3',
  2272. row['File']+'_imec0_spikes.txt'))
  2273. electrodes = Electrodes(join(row['Base_Directory'], 'Preprocessed_Data/Electrodes',
  2274. row['File']+'_imec0_electrodes.txt'))
  2275. else:
  2276. mec_spikes = Spikes(join(row['Base_Directory'], 'Preprocessed_Data/Spikes/g1_3',
  2277. row['File']+'_imec1_spikes.txt'))
  2278. electrodes = Electrodes(join(row['Base_Directory'], 'Preprocessed_Data/Electrodes',
  2279. row['File']+'_imec1_electrodes.txt'))
  2280. MEC_channels = electrodes.subset_by_location(regions=['ENTm1','ENTm2','ENTm3','ENTm4','ENTm5','ENTm6'])
  2281. mec_spikes.subset_by_channel(MEC_channels)
  2282. else:
  2283. mec_spikes = Spikes(join(row['Base_Directory'], 'Preprocessed_Data/Spikes/g1_3',
  2284. row['File']+'_imec0_spikes.txt'))
  2285. electrodes = Electrodes(join(row['Base_Directory'], 'Preprocessed_Data/Electrodes',
  2286. row['File']+'_electrodes.txt'))
  2287. channels = electrodes.subset_by_location(regions=['Entorhinal area medial part dorsal zone'])
  2288. mec_spikes.subset_by_channel(channels)
  2290. spatial_sig, stability_sig, grid_sig, border_sig, speed_sig, HD_sig = \
  2291. [np.full(len(mec_spikes.spikes), np.NaN) for _ in range(6)]
  2293. for u in range(len(mec_spikes.spikes)):
  2294. fr = bin_spikes(mec_spikes.spikes[u], position.position.Timestamp)
  2295. _, spatial_sig[u] = permutation_test_spatial_info(fr, position, bin_cm)
  2296. # _, stability_sig[u] = permutation_test_spatial_stability(fr, position, bin_cm)
  2297. try: # sometimes causes timeout error
  2298. _, grid_sig[u], _, _, _ = permutation_test_grid_score(fr, position, bin_cm)
  2299. except:
  2300. pass
  2301. _, border_sig[u] = permutation_test_border_score(fr, position, bin_cm)
  2302. _, speed_sig[u] = speed_score_traditional(position.position.Velocity, fr)
  2303. _, HD_sig[u], _, _, _, _ = angle_score_traditional(position.position.HD, fr)
  2305. spatial_cells[animal_idx][sess_idx] = np.where(spatial_sig<0.05)[0]
  2306. #np.intersect1d(np.where(spatial_sig<0.05)[0], np.where(stability_sig<0.05)[0])
  2307. grid_cells[animal_idx][sess_idx] = np.where(grid_sig<0.05)[0]
  2308. border_cells[animal_idx][sess_idx] = np.where(border_sig<0.05)[0]
  2309. spatial_aperiodic_cells[animal_idx][sess_idx] = np.array(list(set(spatial_cells[animal_idx][sess_idx]) - \
  2310. set(grid_cells[animal_idx][sess_idx]) - \
  2311. set(border_cells[animal_idx][sess_idx])))
  2312. speed_cells[animal_idx][sess_idx] = np.where(speed_sig<0.05)[0]
  2313. hd_cells[animal_idx][sess_idx] = np.where(HD_sig<0.05)[0]
  2315. print(f"{row['File']} {time.time() - start_time} seconds")
  2317. # %%
  2318. # enrichment: are any of these cells on the enrichment lists?
  2319. enriched_cells = over_fr_list
  2321. enriched_spatial, enriched_grid, enriched_border, enriched_spatial_aperiodic, \
  2322. enriched_speed, enriched_hd = [np.full((len(animal_dict),len(session_range)), np.NaN) for _ in range(6)]
  2323. for animal_idx in range(n_animals):
  2324. for sess_idx in range(n_sessions):
  2325. # proportion of cells in enrichment list
  2326. enriched_spatial[animal_idx][sess_idx] = len(np.intersect1d(spatial_cells[animal_idx][sess_idx], enriched_cells[animal_idx][sess_idx]))/ \
  2327. len(spatial_cells[animal_idx][sess_idx]) if len(spatial_cells[animal_idx][sess_idx])>0 else np.NaN
  2328. enriched_grid[animal_idx][sess_idx] = len(np.intersect1d(grid_cells[animal_idx][sess_idx], enriched_cells[animal_idx][sess_idx]))/ \
  2329. len(grid_cells[animal_idx][sess_idx]) if len(grid_cells[animal_idx][sess_idx])>0 else np.NaN
  2330. enriched_border[animal_idx][sess_idx] = len(np.intersect1d(border_cells[animal_idx][sess_idx], enriched_cells[animal_idx][sess_idx]))/ \
  2331. len(border_cells[animal_idx][sess_idx]) if len(border_cells[animal_idx][sess_idx])>0 else np.NaN
  2332. enriched_spatial_aperiodic[animal_idx][sess_idx] = len(np.intersect1d(spatial_aperiodic_cells[animal_idx][sess_idx], enriched_cells[animal_idx][sess_idx]))/ \
  2333. len(spatial_aperiodic_cells[animal_idx][sess_idx]) if len(spatial_aperiodic_cells[animal_idx][sess_idx])>0 else np.NaN
  2334. enriched_speed[animal_idx][sess_idx] = len(np.intersect1d(speed_cells[animal_idx][sess_idx], enriched_cells[animal_idx][sess_idx]))/ \
  2335. len(speed_cells[animal_idx][sess_idx]) if len(speed_cells[animal_idx][sess_idx])>0 else np.NaN
  2336. enriched_hd[animal_idx][sess_idx] = len(np.intersect1d(hd_cells[animal_idx][sess_idx], enriched_cells[animal_idx][sess_idx]))/ \
  2337. len(hd_cells[animal_idx][sess_idx]) if len(hd_cells[animal_idx][sess_idx])>0 else np.NaN
  2339. plot_feature_by_cell_type(enriched_spatial, enriched_grid, enriched_border, enriched_spatial_aperiodic, enriched_speed, enriched_hd, \
  2340. figure_path, 'proportion_enriched_by_cell_type_fr', '% Units with Recruitment > 95th Percentile', [0,1])
  2341. plot_feature_by_cell_type_3types(enriched_spatial, enriched_speed, enriched_hd, \
  2342. figure_path, 'proportion_enriched_by_cell_type_3categories_fr', '% Units with Recruitment > 95th Percentile', [0,1])
  2343. plot_feature_by_cell_type_3types(enriched_grid, enriched_border, enriched_spatial_aperiodic, \
  2344. figure_path, 'proportion_enriched_by_cell_type_3spatialcategories_fr', '% Units with Recruitment > 95th Percentile', [0,1])
  2346. # %%
  2347. # pie chart of categorzied % of enriched cells that represent each cell type (including overlaps)
  2348. # pie chart of % of categorized cells open field that represent each cell type (including overlaps)
  2349. from functools import reduce
  2351. enriched_cells = over_fr_list
  2352. enriched_frac = over_fr
  2354. enriched_spatial, enriched_grid, enriched_border, enriched_spatial_aperiodic, \
  2355. enriched_speed, enriched_hd, enriched_spatial_speed, enriched_spatial_hd, enriched_speed_hd, \
  2356. enriched_spatial_speed_hd, \
  2357. spatial, grid, border, spatial_aperiodic, \
  2358. speed, hd, spatial_speed, spatial_hd, speed_hd, \
  2359. spatial_speed_hd, = \
  2360. [np.full((len(animal_dict),n_sessions), np.NaN) for _ in range(20)]
  2361. for a in range(n_animals):
  2362. for s in range(n_sessions):
  2363. n_cells = len(over_fr_list[a][s])/over_fr[a][s]
  2364. if len(enriched_cells[a][s])>0:
  2365. # proportion of enriched cells of each type
  2366. enriched_grid[a][s] = len(np.intersect1d(grid_cells[a][s], enriched_cells[a][s]))/ \
  2367. len(enriched_cells[a][s])
  2368. enriched_border[a][s] = len(np.intersect1d(border_cells[a][s], enriched_cells[a][s]))/ \
  2369. len(enriched_cells[a][s])
  2370. enriched_spatial_aperiodic[a][s] = len(np.intersect1d(spatial_aperiodic_cells[a][s], enriched_cells[a][s]))/ \
  2371. len(enriched_cells[a][s])
  2373. enriched_spatial_speed_hd[a][s] = len(reduce(np.intersect1d, (spatial_cells[a][s], speed_cells[a][s], hd_cells[a][s], enriched_cells[a][s])))/ \
  2374. len(enriched_cells[a][s])
  2375. enriched_spatial_speed[a][s] = len(reduce(np.intersect1d, (spatial_cells[a][s], speed_cells[a][s], enriched_cells[a][s])))/ \
  2376. len(enriched_cells[a][s]) - enriched_spatial_speed_hd[a][s]
  2377. enriched_spatial_hd[a][s] = len(reduce(np.intersect1d, (spatial_cells[a][s], hd_cells[a][s], enriched_cells[a][s])))/ \
  2378. len(enriched_cells[a][s]) - enriched_spatial_speed_hd[a][s]
  2379. enriched_speed_hd[a][s] = len(reduce(np.intersect1d, (speed_cells[a][s], hd_cells[a][s], enriched_cells[a][s])))/ \
  2380. len(enriched_cells[a][s]) - enriched_spatial_speed_hd[a][s]
  2381. enriched_spatial[a][s] = len(np.intersect1d(spatial_cells[a][s], enriched_cells[a][s]))/len(enriched_cells[a][s]) \
  2382. - enriched_spatial_speed_hd[a][s] - enriched_spatial_speed[a][s] - enriched_spatial_hd[a][s]
  2383. enriched_speed[a][s] = len(np.intersect1d(speed_cells[a][s], enriched_cells[a][s]))/len(enriched_cells[a][s]) \
  2384. - enriched_spatial_speed_hd[a][s] - enriched_spatial_speed[a][s] - enriched_speed_hd[a][s]
  2385. enriched_hd[a][s] = len(np.intersect1d(hd_cells[a][s], enriched_cells[a][s]))/len(enriched_cells[a][s]) \
  2386. - enriched_spatial_speed_hd[a][s] - enriched_spatial_hd[a][s] - enriched_spatial_hd[a][s]
  2388. # proportion of all cells of each type
  2389. grid[a][s] = len(grid_cells[a][s])/n_cells
  2390. border[a][s] = len(border_cells[a][s])/n_cells
  2391. spatial_aperiodic[a][s] = len(spatial_aperiodic_cells[a][s])/n_cells
  2393. spatial_speed_hd[a][s] = len(reduce(np.intersect1d, (spatial_cells[a][s], speed_cells[a][s], hd_cells[a][s])))/ \
  2394. n_cells
  2395. spatial_speed[a][s] = len(np.intersect1d(spatial_cells[a][s], speed_cells[a][s]))/ \
  2396. n_cells - spatial_speed_hd[a][s]
  2397. spatial_hd[a][s] = len(np.intersect1d(spatial_cells[a][s], hd_cells[a][s]))/ \
  2398. n_cells - spatial_speed_hd[a][s]
  2399. speed_hd[a][s] = len(np.intersect1d(speed_cells[a][s], hd_cells[a][s]))/ \
  2400. n_cells - spatial_speed_hd[a][s]
  2401. spatial[a][s] = len(spatial_cells[a][s])/n_cells \
  2402. - spatial_speed_hd[a][s] - spatial_speed[a][s] - spatial_hd[a][s]
  2403. speed[a][s] = len(speed_cells[a][s])/n_cells \
  2404. - spatial_speed_hd[a][s] - spatial_speed[a][s] - speed_hd[a][s]
  2405. hd[a][s] = len(hd_cells[a][s])/n_cells \
  2406. - spatial_speed_hd[a][s] - spatial_hd[a][s] - spatial_hd[a][s]
  2408. f, ax = plt.subplots(4)
  2410. ax[0].pie([np.nansum(enriched_spatial), np.nansum(enriched_speed), np.nansum(enriched_hd), \
  2411. np.nansum(enriched_spatial_speed), np.nansum(enriched_spatial_hd), np.nansum(enriched_speed_hd), \
  2412. np.nansum(enriched_spatial_speed_hd)], \
  2413. labels=['Position','Speed','HD','Position x Speed', 'Position x HD', 'Speed x HD', 'Position x Speed x HD'])
  2414. ax[1].pie([np.nansum(enriched_spatial_aperiodic), \
  2415. np.nansum(enriched_grid), \
  2416. np.nansum(enriched_border)], \
  2417. labels=['Spatial Aperiodic', 'Grid', 'Border'])
  2418. ax[2].pie([np.nansum(spatial), np.nansum(speed), np.nansum(hd), \
  2419. np.nansum(spatial_speed), np.nansum(spatial_hd), np.nansum(speed_hd), \
  2420. np.nansum(spatial_speed_hd)], \
  2421. labels=['Position','Speed','HD','Position x Speed', 'Position x HD', 'Speed x HD', 'Position x Speed x HD'])
  2422. ax[3].pie([np.nansum(spatial_aperiodic), \
  2423. np.nansum(grid), \
  2424. np.nansum(border)], \
  2425. labels=['Spatial Aperiodic', 'Grid', 'Border'])
  2426. f.savefig(join(figure_path,'cell_types_pie_chart.pdf'), format='pdf')
  2428. # %%
  2429. # representative plots of spatial variable coding
  2430. # pick a session
  2431. row = sessions.iloc[17]
  2432. print(row['File'])
  2433. position = Position(name=join(row['Base_Directory'], 'Preprocessed_Data/Position',
  2434. row['File']+'_position.txt'))
  2435. if ("2023_spring" in row['Base_Directory']) or ("2024_winter" in row['Base_Directory']):
  2436. if (row['Animal']=='Lamarr') and (int(row['Session'][-2:])>=8):
  2437. mec_spikes = Spikes(name=join(row['Base_Directory'], 'Preprocessed_Data/Spikes/g1_3',
  2438. row['File']+'_imec0_spikes.txt'))
  2439. electrodes = Electrodes(join(row['Base_Directory'], 'Preprocessed_Data/Electrodes',
  2440. row['File']+'_imec0_electrodes.txt'))
  2441. else:
  2442. mec_spikes = Spikes(join(row['Base_Directory'], 'Preprocessed_Data/Spikes/g1_3',
  2443. row['File']+'_imec1_spikes.txt'))
  2444. electrodes = Electrodes(join(row['Base_Directory'], 'Preprocessed_Data/Electrodes',
  2445. row['File']+'_imec1_electrodes.txt'))
  2446. MEC_channels = electrodes.subset_by_location(regions=['ENTm1','ENTm2','ENTm3','ENTm4','ENTm5','ENTm6'])
  2447. mec_spikes.subset_by_channel(MEC_channels)
  2448. else:
  2449. mec_spikes = Spikes(join(row['Base_Directory'], 'Preprocessed_Data/Spikes/g1_3',
  2450. row['File']+'_imec0_spikes.txt'))
  2451. electrodes = Electrodes(join(row['Base_Directory'], 'Preprocessed_Data/Electrodes',
  2452. row['File']+'_electrodes.txt'))
  2453. channels = electrodes.subset_by_location(regions=['Entorhinal area medial part dorsal zone'])
  2454. mec_spikes.subset_by_channel(channels)
  2456. sess = "BaggySweatpants_DY08"
  2457. # spatial firing
  2458. unit = 67
  2459. #[ 0, 136, 149] #[ 63, 67, 76, 95, 111, 163] #[70, 72, 43, 44, 46, 48, 54, 89, 27, 94]
  2460. #0, 2, 8, 117, 125, 136, 141, 149, 156, 180
  2461. fr_map, xbins, ybins, max_fr = calc_fr_map(mec_spikes.spikes[unit], position, bin_cm, smooth=True)
  2462. f = plot_heatmap(fr_map, xbins, ybins)
  2463. f.savefig(join(figure_path,sess+'_grid_example.pdf'), format='pdf', transparent=True,
  2464. dpi=300, bbox_inches='tight')
  2465. print(max_fr)
  2466. unit = 46
  2467. #28, 33, 63, 67, 76, 93, 95, 99, 100, 108, 111, 123, 130, 139, 147, 164, 169, 173, 179
  2468. fr_map, xbins, ybins, max_fr = calc_fr_map(mec_spikes.spikes[unit], position, bin_cm, smooth=True)
  2469. f = plot_heatmap(fr_map, xbins, ybins)
  2470. f.savefig(join(figure_path,sess+'_border_example.pdf'), format='pdf', transparent=True,
  2471. dpi=300, bbox_inches='tight')
  2472. print(max_fr)
  2473. unit = 94
  2474. #[ 11, 140, 142, 14, 22, 153, 26, 27, 25, 29, 161, 38, 43, 44, 45, 46, 48, 52, 54, 55, 58, 64, 70, 72, 78, 89, 91, 94, 107, 115]
  2475. fr_map, xbins, ybins, max_fr = calc_fr_map(mec_spikes.spikes[unit], position, bin_cm, smooth=True)
  2476. f = plot_heatmap(fr_map, xbins, ybins)
  2477. f.savefig(join(figure_path,sess+'_spatial_aperiodic_example.pdf'), format='pdf', transparent=True,
  2478. dpi=300, bbox_inches='tight')
  2479. print(max_fr)
  2481. # stability
  2482. fr = bin_spikes(mec_spikes.spikes[unit], position.position.Timestamp)
  2483. n_samps = len(fr)
  2484. split_idx = int(n_samps/2)
  2485. first_half_pos = deepcopy(position)
  2486. first_half_pos.position = first_half_pos.position.iloc[:split_idx]
  2487. second_half_pos = deepcopy(position)
  2488. second_half_pos.position = second_half_pos.position.iloc[split_idx:]
  2489. fr_map, xbins, ybins, max_fr = calc_fr_map(fr[:split_idx], first_half_pos, bin_cm, smooth=True, spikes_are_binned=True)
  2490. f = plot_heatmap(fr_map, xbins, ybins)
  2491. f.savefig(join(figure_path,sess+'_spatial_aperiodic_example_first_half.pdf'), format='pdf', transparent=True,
  2492. dpi=300, bbox_inches='tight')
  2493. fr_map, xbins, ybins, max_fr = calc_fr_map(fr[split_idx:], second_half_pos, bin_cm, smooth=True, spikes_are_binned=True)
  2494. f = plot_heatmap(fr_map, xbins, ybins)
  2495. f.savefig(join(figure_path,sess+'_spatial_aperiodic_example_second_half.pdf'), format='pdf', transparent=True,
  2496. dpi=300, bbox_inches='tight')
  2498. # speed
  2499. unit = 29 #[ 12, 29, 30, 48, 50, 53, 69, 72, 87, 98, 103, 126, 133, 137, 144, 146, 148, 161, 166, 171, 182, 192]
  2500. fr = bin_spikes(mec_spikes.spikes[unit], position.position.Timestamp)*60
  2501. bins = np.arange(0,100,5)
  2502. fr_by_speed = pd.DataFrame({'FR': fr, 'Speed': pd.cut(position.position.Velocity, bins=bins)})
  2503. fr_by_speed_mean = fr_by_speed.groupby(['Speed']).mean().reset_index()
  2504. fr_by_speed_sem = fr_by_speed.groupby(['Speed']).sem().reset_index()
  2505. f, ax = plt.subplots()
  2506. ax.plot(bins[:-1], fr_by_speed_mean.FR)
  2507. ax.fill_between(bins[:-1], fr_by_speed_mean.FR - fr_by_speed_sem.FR, fr_by_speed_mean.FR + fr_by_speed_sem.FR, \
  2508. color='grey', alpha=0.5, rasterized=True)
  2509. ax.set(xlim=[0,50])
  2510. f.savefig(join(figure_path,sess+'_speed_example.pdf'), format='pdf', transparent=True,
  2511. dpi=300, bbox_inches='tight')
  2513. # HD
  2514. unit = 32 #[ 10, 11, 14, 16, 19, 20, 22, 24, 26, 28, 32, 47, 49,
  2515. #53, 55, 62, 65, 66, 79, 90, 103, 107, 111, 121, 130, 139,
  2516. #150, 156, 164, 167, 175, 193, 194, 197]
  2517. fr = bin_spikes(mec_spikes.spikes[unit], position.position.Timestamp)*60
  2518. bins = np.linspace(-1*math.pi,math.pi,18)
  2519. fr_by_hd = pd.DataFrame({'FR': fr, 'HD': pd.cut(position.position.HD, bins=bins)})
  2520. fr_by_hd_mean = fr_by_hd.groupby(['HD']).mean().reset_index()
  2521. fr_by_hd_sem = fr_by_hd.groupby(['HD']).sem().reset_index()
  2522. means = np.append(fr_by_hd_mean.FR,fr_by_hd_mean.FR[0])
  2523. sems = np.append(fr_by_hd_sem.FR,fr_by_hd_sem.FR[0])
  2524. f, ax = plt.subplots(subplot_kw={'projection': 'polar'})
  2525. ax.plot(bins, means)
  2526. ax.fill_between(bins, means - sems, means + sems, \
  2527. color='grey', alpha=0.5, rasterized=True)
  2528. f.savefig(join(figure_path,sess+'_HD_example.pdf'), format='pdf', transparent=True,
  2529. dpi=300, bbox_inches='tight')
  2531. # %% [markdown]
  2532. # ### Reward coding
  2534. # %%
  2535. # find cells whose firing rates are significantly increased at every reward location
  2536. sessions = pd.read_csv('//oak-smb-giocomo.stanford.edu/groups/giocomo/emijones/WT_Sequences/all_sessions_Xmaze.csv')
  2537. n_shuffles = 1000
  2538. n_sessions = 20
  2539. n_animals = len(animal_dict)
  2540. all_rewards = np.full((n_animals, n_sessions), np.NaN)
  2541. reward_increase = np.full((n_animals, n_sessions), np.NaN)
  2542. all_rewards_list = [[[] for _ in range(n_sessions)] for _ in range(n_animals)]
  2543. reward_increase_list = [[[] for _ in range(n_sessions)] for _ in range(n_animals)]
  2545. for _, row in sessions.iterrows():
  2546. sess_idx = int(row['Session'][-2:])-1
  2547. if not row['Recording_Error'] and not row['Position_Error'] and \
  2548. row['Task'] == 'X Maze' and sess_idx<=10:
  2549. # load spikes
  2550. binned_data = Binned_Data(join(row['Base_Directory'],
  2551. 'Preprocessed_Data/Sequences/Binned_Data',
  2552. row['File']+'_MEC_binned_data.txt'))
  2553. position = Position(join(row['Base_Directory'], 'Preprocessed_Data/Position',
  2554. row['File']+'_position.txt'))
  2555. task = DoubleYMazeTask(join(row['Base_Directory'], 'Preprocessed_Data/Task',
  2556. row['File']+'_task.txt'))
  2557. hd_shift = position.position.HD+math.pi
  2558. animal_idx = animal_dict[row['Animal']]
  2560. # calculate reward FR for each cell at each rewarded position
  2561. # set start and end indices and reward locations for later shuffles
  2562. # loop over reward visits
  2563. true_fr_by_loc = np.zeros((binned_data.spikes.shape[0],4))
  2564. true_fr_rew_delta = np.zeros(binned_data.spikes.shape[0])
  2565. reward_location, reward_start_idx, reward_end_idx = [np.zeros(len(task.trials), dtype=int) for _ in range(3)]
  2566. num_spikes_by_loc = np.zeros((binned_data.spikes.shape[0],4))
  2567. reward_dur_by_loc = np.zeros(4)
  2568. num_spikes_rew_or_not = np.zeros((binned_data.spikes.shape[0],2))
  2569. reward_dur_rew_or_not = np.zeros(2)
  2570. for t, trial in task.trials.iterrows():
  2571. # find timepoints of first entering poke + turning around to leave poke
  2572. reward_start_idx[t] = np.argmin(np.abs(position.position.Timestamp-trial.Start))
  2573. # if recording ended before mouse left reward, take last index
  2574. if (np.where(np.abs(hd_shift[reward_start_idx[t]:]-hd_shift[reward_start_idx[t]])>math.pi/2)[0]).size==0:
  2575. reward_end_idx[t] = len(hd_shift)-1
  2576. else:
  2577. reward_end_idx[t] = np.where(np.abs(hd_shift[reward_start_idx[t]:]-hd_shift[reward_start_idx[t]])>math.pi/2)[0][0] + reward_start_idx[t]
  2579. # which reward location
  2580. # TL, BL, TR, BR
  2581. if (trial['Inbound/Outbound']=='Outbound') and (trial['Trajectory']=='Top'):
  2582. reward_location[t] = 0
  2583. elif (trial['Inbound/Outbound']=='Outbound') and (trial['Trajectory']=='Bottom'):
  2584. reward_location[t] = 1
  2585. elif (trial['Inbound/Outbound']=='Inbound') and (trial['Trajectory']=='Top'):
  2586. reward_location[t] = 2
  2587. else:
  2588. reward_location[t] = 3
  2590. num_spikes_by_loc[:, reward_location[t]] += np.sum(binned_data.spikes[:, reward_start_idx[t]:reward_end_idx[t]], axis=1)
  2591. reward_dur_by_loc[reward_location[t]] += position.position.Timestamp[reward_end_idx[t]] - position.position.Timestamp[reward_start_idx[t]]
  2593. # choice reward, which could be rewarded or not
  2594. if (reward_location[t]==2) or (reward_location[t]==3):
  2595. if trial['Correct']:
  2596. num_spikes_rew_or_not[:,0] += num_spikes_by_loc[:, reward_location[t]]
  2597. reward_dur_rew_or_not[0] += reward_dur_by_loc[reward_location[t]]
  2598. else:
  2599. num_spikes_rew_or_not[:,1] += num_spikes_by_loc[:, reward_location[t]]
  2600. reward_dur_rew_or_not[1] += reward_dur_by_loc[reward_location[t]]
  2602. for r in range(4):
  2603. true_fr_by_loc[:,r] = num_spikes_by_loc[:,r]/reward_dur_by_loc[r]
  2604. true_fr_rew_delta = num_spikes_rew_or_not[:,0]/reward_dur_rew_or_not[0] - num_spikes_rew_or_not[:,1]/reward_dur_rew_or_not[1]
  2606. # loop over shuffles
  2607. shuffled_fr_list_by_loc = np.zeros((n_shuffles, binned_data.spikes.shape[0], 4))
  2608. shuffled_fr_list_rew_delta = np.zeros((n_shuffles, binned_data.spikes.shape[0]))
  2609. for s in range(n_shuffles):
  2610. # shuffle each cell separately
  2611. shuffled_spikes = binned_data.spikes[:, np.random.permutation(binned_data.spikes.shape[1])]
  2613. # calculate shuffled fr during rewards
  2614. num_spikes_by_loc = np.zeros((binned_data.spikes.shape[0],4))
  2615. reward_dur_by_loc = np.zeros(4)
  2616. num_spikes_rew_or_not = np.zeros((binned_data.spikes.shape[0],2))
  2617. reward_dur_rew_or_not = np.zeros(2)
  2618. for t, trial in task.trials.iterrows():
  2619. num_spikes_by_loc[:, reward_location[t]] += np.sum(shuffled_spikes[:, reward_start_idx[t]:reward_end_idx[t]], axis=1)
  2620. reward_dur_by_loc[reward_location[t]] += position.position.Timestamp[reward_end_idx[t]] - position.position.Timestamp[reward_start_idx[t]]
  2622. # choice reward, which could be rewarded or not
  2623. if (reward_location[t]==2) or (reward_location[t]==3):
  2624. if trial['Correct']:
  2625. num_spikes_rew_or_not[:,0] += num_spikes_by_loc[:, reward_location[t]]
  2626. reward_dur_rew_or_not[0] += reward_dur_by_loc[reward_location[t]]
  2627. else:
  2628. num_spikes_rew_or_not[:,1] += num_spikes_by_loc[:, reward_location[t]]
  2629. reward_dur_rew_or_not[1] += reward_dur_by_loc[reward_location[t]]
  2631. # append to shuffle distribution
  2632. for r in range(4):
  2633. shuffled_fr_list_by_loc[s,:,r] = num_spikes_by_loc[:,r]/reward_dur_by_loc[r]
  2634. shuffled_fr_list_rew_delta[s,:] = num_spikes_rew_or_not[:,0]/reward_dur_rew_or_not[0] - num_spikes_rew_or_not[:,1]/reward_dur_rew_or_not[1]
  2636. # calculate 95th percentile
  2637. fr_95_percentile_by_loc = np.zeros((binned_data.spikes.shape[0],4))
  2638. fr_95_percentile_rew_delta = np.zeros((binned_data.spikes.shape[0]))
  2639. for r in range(4):
  2640. fr_95_percentile_by_loc[:,r] = np.percentile(shuffled_fr_list_by_loc[:,:,r], 95, axis=0)
  2641. fr_95_percentile_rew_delta = np.percentile(shuffled_fr_list_rew_delta, 95, axis=0)
  2643. # count how many cells are above this threshold for all 4 rewards or delta FR (rewarded vs not)
  2644. all_rewards[animal_idx, sess_idx] = np.nanmean((true_fr_by_loc[:,0]>fr_95_percentile_by_loc[:,0]) &
  2645. (true_fr_by_loc[:,1]>fr_95_percentile_by_loc[:,1]) &
  2646. (true_fr_by_loc[:,2]>fr_95_percentile_by_loc[:,2]) &
  2647. (true_fr_by_loc[:,3]>fr_95_percentile_by_loc[:,3]))
  2648. all_rewards_list[animal_idx][sess_idx] = np.where((true_fr_by_loc[:,0]>fr_95_percentile_by_loc[:,0]) &
  2649. (true_fr_by_loc[:,1]>fr_95_percentile_by_loc[:,1]) &
  2650. (true_fr_by_loc[:,2]>fr_95_percentile_by_loc[:,2]) &
  2651. (true_fr_by_loc[:,3]>fr_95_percentile_by_loc[:,3]))[0].tolist()
  2652. reward_increase[animal_idx][sess_idx] = np.nanmean(true_fr_rew_delta>fr_95_percentile_rew_delta)
  2653. reward_increase_list[animal_idx][sess_idx] = np.where(true_fr_rew_delta>fr_95_percentile_rew_delta)[0].tolist()
  2655. print(row['File'])
  2657. # %%
  2658. with open(join(figure_path,'all_rewards.pkl'), 'wb') as file:
  2659. pickle.dump(all_rewards, file)
  2660. with open(join(figure_path,'all_rewards_list.pkl'), 'wb') as file:
  2661. pickle.dump(all_rewards_list, file)
  2662. with open(join(figure_path,'reward_increase.pkl'), 'wb') as file:
  2663. pickle.dump(reward_increase, file)
  2664. with open(join(figure_path,'reward_increase_list.pkl'), 'wb') as file:
  2665. pickle.dump(reward_increase_list, file)
  2667. # %%
  2668. print(f"{np.nanmean(all_rewards)*100:.02f} +/- {stats.sem(np.nanmean(all_rewards,axis=1))*100:.02f}" +\
  2669. f"% of cells have > 95th percentile of shuffled FR during all rewards,\n" +\
  2670. f"and {np.nanmean(reward_increase)*100:.02f} +/- {stats.sem(np.nanmean(reward_increase,axis=1))*100:.02f}" +\
  2671. f"% of cells have > 95th percentile of shuffled increase in FR during reward consumption vs no reward.")
  2673. # %%
  2674. with open(join(figure_path,'over_fr_list.pkl'), 'rb') as file:
  2675. over_fr_list = pickle.load(file)
  2677. # %%
  2678. # are preferentially recruited cells more likely to be reward coding?
  2679. n_sessions = 10
  2680. reward_enriched, reward_inc_enriched, reward_unenriched, reward_inc_unenriched = \
  2681. [np.full((n_animals, n_sessions), np.NaN) for _ in range(4)]
  2682. for animal_idx in range(n_animals):
  2683. for sess_idx in range(n_sessions):
  2684. # what proportion of recruited cells are reward coding
  2685. # number of reward coding cells in enriched list divided by number of total cells in enriched list
  2686. reward_enriched[animal_idx][sess_idx] = len(np.intersect1d(over_fr_list[animal_idx][sess_idx], all_rewards_list[animal_idx][sess_idx]))/ \
  2687. len(over_fr_list[animal_idx][sess_idx]) if len(over_fr_list[animal_idx][sess_idx])>0 else np.NaN
  2688. reward_inc_enriched[animal_idx][sess_idx] = len(np.intersect1d(over_fr_list[animal_idx][sess_idx], reward_increase_list[animal_idx][sess_idx]))/ \
  2689. len(over_fr_list[animal_idx][sess_idx]) if len(over_fr_list[animal_idx][sess_idx])>0 else np.NaN
  2691. # what proportion of non-recruited cells are reward coding
  2692. # number of reward coding cells not in enriched list divided by number of total cells not in enriched list
  2693. reward_unenriched[animal_idx][sess_idx] = np.sum(~np.isin(all_rewards_list[animal_idx][sess_idx], over_fr_list[animal_idx][sess_idx]))/ \
  2694. (len(all_rewards_list[animal_idx][sess_idx])/all_rewards[animal_idx][sess_idx]) if len(over_fr_list[animal_idx][sess_idx])>0 else np.NaN
  2695. reward_inc_unenriched[animal_idx][sess_idx] = np.sum(~np.isin(reward_increase_list[animal_idx][sess_idx], over_fr_list[animal_idx][sess_idx]))/ \
  2696. (len(all_rewards_list[animal_idx][sess_idx])/all_rewards[animal_idx][sess_idx]) if len(over_fr_list[animal_idx][sess_idx])>0 else np.NaN
  2699. # %%
  2700. # compare
  2701. data_by_session = pd.DataFrame({'Animal': np.repeat(np.arange(n_animals),n_sessions),
  2702. 'Session': np.tile(np.arange(n_sessions),n_animals),
  2703. 'Reward_enriched': reward_enriched.flatten()*100,
  2704. 'Reward_unenriched': reward_unenriched.flatten()*100,
  2705. 'Reward_inc_enriched': reward_inc_enriched.flatten()*100,
  2706. 'Reward_inc_unenriched': reward_inc_unenriched.flatten()*100})
  2707. plot_compare_paired_metrics('Reward_enriched', 'Reward_unenriched', data_by_session, \
  2708. figure_path, 'percent_rew_coding_recruited_vs_not', '% of Cells that Represent Reward', [0,100])
  2709. plot_compare_paired_metrics('Reward_inc_enriched', 'Reward_inc_unenriched', data_by_session, \
  2710. figure_path, 'percent_rew_inc_recruited_vs_not', '% of Cells that Represent Reward Consumption', [0,100])
  2712. # %% [markdown]
  2713. # # Figure S5
  2715. # %% [markdown]
  2716. # ### Decode with SWR trace
  2718. # %%
  2719. sessions = pd.read_csv('Z:/WT_Sequences/all_sessions_Xmaze.csv')
  2720. row = sessions.iloc[156] #64
  2721. print(row['File'])
  2722. seq = Sequences(join(row['Base_Directory'], 'Preprocessed_Data/Sequences',
  2723. row['File']+'_MEC_immobility_sequences.txt'),
  2724. load_objects=True)
  2725. CA1_seq = Sequences(join(row['Base_Directory'], 'Preprocessed_Data/Sequences',
  2726. row['File']+'_CA1_immobility_sequences.txt'),
  2727. load_objects=True)
  2728. SWR_intervals = Sequence_Intervals(join(row['Base_Directory'],
  2729. 'Preprocessed_Data/Sequences/Intervals',
  2730. row['File']+'_CA1_SWR_intervals.txt'))
  2731. lfp = LFP(join(row['Base_Directory'], 'Preprocessed_Data/LFP',
  2732. row['File']+'_CA1_SWR_channel_lfp.txt'))
  2733. ripple_trace = lf.bandpass_filter(lfp.lfp, lfp.rate, lf.RIPPLE[0], lf.RIPPLE[1], lf.RIPPLE[2])
  2735. #seq.stats[seq.stats['Non_local'] & (seq.stats['Duration']<5)]
  2736. #SWR_intervals.intervals
  2738. # %%
  2739. # plot decode
  2740. idx = 108
  2741. seq.plot_sequence(indices=[idx], figure_path=f"{figure_path}/Payne_DY02_MEC")
  2742. CA1_seq.plot_sequence(indices=[idx], figure_path=f"{figure_path}/Payne_DY02_CA1")
  2744. # plot SWR trace
  2745. f, ax = plt.subplots(figsize=(6.3, 3))
  2746. ax.plot(lfp.timestamps, ripple_trace)
  2747. ax.set(xlabel='Time (s)', ylabel='Amplitude')
  2749. # label SWRs
  2750. for ripple in SWR_intervals.intervals.itertuples():
  2751. ax.axvspan(ripple.start_time, ripple.end_time, alpha=0.3, zorder=2)
  2753. # display only overlap with sequence
  2754. start_time = seq.intervals.intervals['start_time'].iloc[idx]
  2755. end_time = seq.intervals.intervals['end_time'].iloc[idx]
  2756. ax.set(xlim=[start_time, end_time])
  2758. f.savefig(join(figure_path,'Payne_DY02_SWR_example_ripple_trace_immobinterval108.pdf'), format='pdf')
  2760. # %% [markdown]
  2761. # ### % of nonlocal coding overlapping with SWRs and vice versa
  2762. # SWRs are not enriched during nonlocal coding. Nonlocal coding is equally likely to occur during an SWR as it is during an immobility bout. SWRs only overlap with a small fraction of nonlocal content.
  2764. # %%
  2765. # averaged over sessions
  2766. sessions = pd.read_csv('Z:/WT_Sequences/all_sessions_Xmaze_CA1.csv')
  2767. n_sessions = 20
  2768. n_animals = len(animal_dict)
  2769. ripple_overlap = np.zeros((n_animals, n_sessions))
  2770. ripple_time = np.full((n_animals, n_sessions), np.NaN)
  2772. for _, row in sessions.iterrows():
  2773. if not row['Recording_Error'] and not row['Position_Error'] and \
  2774. row['Task'] == 'X Maze' and int(row['Session'][-2:])<=20:
  2776. animal_idx = animal_dict[row['Animal']]
  2777. sess_idx = int(row['Session'][-2:])-1
  2779. # load intervals
  2780. nonlocal_immobility_intervals = Sequence_Intervals(join(row['Base_Directory'],
  2781. 'Preprocessed_Data/Sequences/Intervals',
  2782. row['File']+'_MEC_nonlocal_immobility_intervals.txt'))
  2783. immobility_intervals = Sequence_Intervals(join(row['Base_Directory'],
  2784. 'Preprocessed_Data/Sequences/Intervals',
  2785. row['File']+'_immobility_intervals.txt'))
  2786. SWR_intervals = Sequence_Intervals(join(row['Base_Directory'],
  2787. 'Preprocessed_Data/Sequences/Intervals',
  2788. row['File']+'_CA1_SWR_intervals.txt'))
  2790. # iterate over nonlocal intervals, then SWRs
  2791. for nl_immob in nonlocal_immobility_intervals.intervals.itertuples():
  2792. # find overlaps
  2793. for ripple in SWR_intervals.intervals.itertuples():
  2794. if ripple.start_time<nl_immob.end_time and ripple.end_time>nl_immob.start_time:
  2795. if ripple.start_time<=nl_immob.start_time:
  2796. ripple_overlap[animal_idx,sess_idx] += ripple.end_time-nl_immob.start_time
  2797. elif ripple.end_time<=nl_immob.end_time:
  2798. ripple_overlap[animal_idx,sess_idx] += ripple.end_time-ripple.start_time
  2799. else:
  2800. ripple_overlap[animal_idx,sess_idx] += nl_immob.end_time-ripple.start_time
  2801. elif ripple.end_time>end_time:
  2802. break
  2803. ripple_overlap[animal_idx,sess_idx] /= np.sum(nonlocal_immobility_intervals.intervals.end_time -\
  2804. nonlocal_immobility_intervals.intervals.start_time)
  2805. ripple_time[animal_idx,sess_idx] = (np.sum(SWR_intervals.intervals.end_time - SWR_intervals.intervals.start_time))/ \
  2806. (np.sum(immobility_intervals.intervals.end_time -\
  2807. immobility_intervals.intervals.start_time))
  2809. # remove sessions that don't have SWRs (e.g. no CA1 data)
  2810. ripple_overlap = np.where(ripple_overlap==0, np.nan, ripple_overlap)
  2811. ripple_overlap *= 100
  2812. ripple_time *=100
  2814. # %%
  2815. print(f"{np.nanmean(np.nanmean(SWR_data_by_session['non_local'])):.02f} +/- {stats.sem(np.nanmean(SWR_data_by_session['non_local'],axis=1), nan_policy='omit'):.02f}" +\
  2816. f"% of SWRs have nonlocal coding in MEC,\n" +\
  2817. f"and {np.nanmean(np.nanmean(ripple_overlap)):.02f} +/- {stats.sem(np.nanmean(ripple_overlap,axis=1), nan_policy='omit'):.02f}" +\
  2818. f"% of nonlocal decoded time overlaps with SWRs, \n" +\
  2819. f"and {np.nanmean(np.nanmean(ripple_time)):.02f} +/- {stats.sem(np.nanmean(ripple_time,axis=1), nan_policy='omit'):.02f}" +\
  2820. f"% of time during immobility is spent having SWRs.")
  2822. # %%
  2823. # plot
  2824. plot_seq_metrics_over_sessions(n_animals, ripple_overlap, figure_path,
  2825. 'nl_overlap_with_swr', '% NonLocal Decode Overlapping with SWRs', [0,100])
  2826. plot_seq_metrics_over_sessions(n_animals, SWR_data_by_session['non_local'], figure_path,
  2827. 'swr_overlap_with_nl', '% SWRs with NonLocal Decode', [0,100])
  2829. # %% [markdown]
  2830. # ### compare to shuffle: % nonlocal coding overlapping with SWRs and vice versa
  2832. # %%
  2833. n_shuffles = 1000
  2834. rng = np.random.default_rng()
  2836. n_sessions = 20
  2837. n_animals = len(animal_dict)
  2838. nonlocal_contain_SWRs_5, nonlocal_contain_SWRs_95, nonlocal_contain_SWRs_median, \
  2839. SWRs_contain_nonlocal_5, SWRs_contain_nonlocal_95, SWRs_contain_nonlocal_median =\
  2840. [np.full((n_animals, n_sessions), np.NaN) for _ in range(6)]
  2842. sessions = pd.read_csv('Z:/WT_Sequences/all_sessions_Xmaze_CA1.csv')
  2843. for _, row in sessions.iterrows():
  2844. if not row['Recording_Error'] and not row['Position_Error'] and \
  2845. row['Task'] == 'X Maze' and int(row['Session'][-2:])<=20:
  2847. animal_idx = animal_dict[row['Animal']]
  2848. sess_idx = int(row['Session'][-2:])-1
  2850. seq = Sequences(join(row['Base_Directory'],
  2851. 'Preprocessed_Data/Sequences',
  2852. row['File']+'_MEC_immobility_sequences.txt'),
  2853. load_objects=True)
  2854. dist = seq.calc_dist_from_decode()
  2855. ts = seq.binned_data.timestamps
  2856. SWR_intervals = Sequence_Intervals(join(row['Base_Directory'],
  2857. 'Preprocessed_Data/Sequences/Intervals',
  2858. row['File']+'_CA1_SWR_intervals.txt'))
  2860. # 1000 shuffles
  2861. perc_SWRs_contain_nonlocal = np.zeros(1000)
  2862. perc_nonlocal_contain_SWRs = np.zeros(1000)
  2863. for b in range(n_shuffles):
  2864. ts_nonlocal = np.array([])
  2865. for ripple in SWR_intervals.indices.itertuples():
  2866. # pick a random immobility interval
  2867. seq_idx = int(rng.random()*len(seq.stats))
  2868. selected_interval = seq.intervals.indices.iloc[seq_idx]
  2869. # pick a random time within that interval of the same duration as the SWR
  2870. SWR_idx_len = ripple.end_index-ripple.start_index
  2871. index_range = selected_interval.end_index-selected_interval.start_index-SWR_idx_len
  2872. shuffled_start_idx = int(rng.random()*index_range)
  2873. ts_subset = ts[shuffled_start_idx:shuffled_start_idx+SWR_idx_len]
  2874. is_nonlocal = dist[shuffled_start_idx:shuffled_start_idx+SWR_idx_len]>20
  2875. count_nonlocal = np.nansum(is_nonlocal)
  2876. ts_nonlocal = np.concatenate((ts_nonlocal, ts_subset[is_nonlocal]))
  2877. if count_nonlocal>0:
  2878. perc_SWRs_contain_nonlocal[b] += 1
  2879. perc_nonlocal_contain_SWRs[b] = len(ts_nonlocal) #len(np.unique(ts_nonlocal))
  2880. perc_SWRs_contain_nonlocal = perc_SWRs_contain_nonlocal/len(SWR_intervals.intervals) * 100
  2881. perc_nonlocal_contain_SWRs = perc_nonlocal_contain_SWRs / \
  2882. (np.sum(nonlocal_immobility_intervals.indices.end_index -\
  2883. nonlocal_immobility_intervals.indices.start_index)) * 100
  2885. SWRs_contain_nonlocal_5[animal_idx, sess_idx] = np.percentile(perc_SWRs_contain_nonlocal, 5)
  2886. SWRs_contain_nonlocal_95[animal_idx, sess_idx] = np.percentile(perc_SWRs_contain_nonlocal, 95)
  2887. SWRs_contain_nonlocal_median[animal_idx, sess_idx] = np.percentile(perc_SWRs_contain_nonlocal, 50)
  2888. nonlocal_contain_SWRs_5[animal_idx, sess_idx] = np.percentile(perc_nonlocal_contain_SWRs, 5)
  2889. nonlocal_contain_SWRs_95[animal_idx, sess_idx] = np.percentile(perc_nonlocal_contain_SWRs, 95)
  2890. nonlocal_contain_SWRs_median[animal_idx, sess_idx] = np.percentile(perc_nonlocal_contain_SWRs, 50)
  2892. # %%
  2893. # % NonLocal Decode Overlapping with SWRs
  2894. n_nonnan = np.count_nonzero(~np.isnan(ripple_overlap))
  2895. n_more_than_95 = np.nansum(nonlocal_contain_SWRs_95 < ripple_overlap)
  2896. print(f"{n_more_than_95} of {n_nonnan} sessions ({n_more_than_95/n_nonnan*100:0.2f}%) had SWRs during non-local content than if SWRs were randomly allocated across immobility.")
  2898. # % SWRs with NonLocal Decode
  2899. n_more_than_95 = np.nansum(SWRs_contain_nonlocal_95 < SWR_data_by_session['non_local'])
  2900. n_nonnan = np.count_nonzero(~np.isnan(SWR_data_by_session['non_local']))
  2901. print(f"{n_more_than_95} of {n_nonnan} sessions ({n_more_than_95/n_nonnan*100:0.2f}%) had more non-local content during SWRs than if SWRs were randomly allocated across immobility.")
  2903. # %%
  2904. f, ax = plt.subplots(1, figsize=(4,4))
  2905. sns.kdeplot(nonlocal_contain_SWRs_95.flatten(), bw_adjust=0.75, color='red')
  2906. ax.set(xlabel='% SWRs with NonLocal Decode', ylabel='Proportion', xlim=[0,100])
  2907. f.savefig(join(figure_path,'swr_overlap_with_nl_over_events_shuffle.pdf'), format='pdf', transparent=True,
  2908. dpi=300, bbox_inches='tight')
  2910. f, ax = plt.subplots(1, figsize=(4,4))
  2911. sns.kdeplot(SWRs_contain_nonlocal_95.flatten(), bw_adjust=0.75, color='red')
  2912. ax.set(xlabel='% NonLocal Decode Overlapping with SWRs', ylabel='Proportion', xlim=[0,100])
  2913. f.savefig(join(figure_path,'nl_overlap_with_swr_over_events_shuffle.pdf'), format='pdf', transparent=True,
  2914. dpi=300, bbox_inches='tight')
  2916. # %% [markdown]
  2917. # ### SWR rate during immobility intervals with any nonlocal content vs only local content
  2919. # %%
  2920. sessions = pd.read_csv('Z:/WT_Sequences/all_sessions_Xmaze_CA1.csv')
  2922. animal = []
  2923. session = []
  2924. ripple_rate = []
  2925. seq_stats = None
  2927. for _, row in sessions.iterrows():
  2928. if not row['Recording_Error'] and not row['Position_Error'] and \
  2929. row['Task'] == 'X Maze' and int(row['Session'][-2:])<=20:
  2931. animal_idx = animal_dict[row['Animal']]
  2932. sess_idx = int(row['Session'][-2:])-1
  2934. # load intervals
  2935. seq = Sequences(join(row['Base_Directory'],
  2936. 'Preprocessed_Data/Sequences',
  2937. row['File']+'_MEC_immobility_sequences.txt'))
  2938. intervals = Sequence_Intervals(join(row['Base_Directory'],
  2939. 'Preprocessed_Data/Sequences/Intervals',
  2940. row['File']+'_immobility_intervals.txt'))
  2941. SWR_intervals = Sequence_Intervals(join(row['Base_Directory'],
  2942. 'Preprocessed_Data/Sequences/Intervals',
  2943. row['File']+'_CA1_SWR_intervals.txt'))
  2944. ripple_rate_session = np.zeros(len(seq.stats))
  2946. # iterate over nonlocal intervals, then SWRs
  2947. if len(SWR_intervals.intervals>0):
  2948. for nl in range(len(seq.stats)):
  2949. # find overlaps
  2950. for ripple in SWR_intervals.intervals.itertuples():
  2951. if ripple.start_time<intervals.intervals.iloc[nl].end_time and \
  2952. ripple.end_time>intervals.intervals.iloc[nl].start_time:
  2953. ripple_rate_session[nl] +=1
  2954. elif ripple.end_time>end_time:
  2955. break
  2956. ripple_rate_session[nl] /= seq.stats.Duration.iloc[nl]
  2958. # build arrays over all sequences
  2959. animal_idx = animal_dict[row['Animal']]
  2960. animal.extend([animal_idx] * len(seq.stats))
  2961. session_idx = int(row['Session'][-2:])-1
  2962. session.extend([session_idx] * len(seq.stats))
  2963. ripple_rate.extend(ripple_rate_session)
  2964. if seq_stats is not None:
  2965. seq_stats = pd.concat([seq_stats, seq.stats], axis=0, ignore_index=True)
  2966. else:
  2967. seq_stats = seq.stats
  2969. # build full df
  2970. ripple_rate_by_seq = pd.DataFrame({'Animal': animal, 'Session': session, 'SWR_rate': ripple_rate})
  2971. ripple_rate_by_seq = pd.concat([ripple_rate_by_seq, seq_stats], axis=1)
  2973. # %%
  2974. plot_seq_metrics_comparison_over_events(ripple_rate_by_seq, 'SWR_rate', 'Non_local',
  2975. figure_path, 'local_vs_nonlocal_SWR_rate', 'SWR rate (Hz)', [0,4])
  2977. # %% [markdown]
  2978. # ### Distance between decoded positions in MEC and CA1 during SWRs
  2980. # %%
  2981. sessions = pd.read_csv(r'Z:/WT_Sequences/all_sessions_Xmaze_CA1.csv')
  2982. n_sessions = 20
  2983. n_animals = len(animal_dict)
  2984. n_shuffles = 1000
  2986. start_time = time.time()
  2987. loc_dist, loc_dist_vs_shuffle = \
  2988. [np.full((n_animals, n_sessions), np.NaN) for _ in range(2)]
  2989. # subset data
  2990. for _, row in sessions.iterrows():
  2991. if not row['Recording_Error'] and not row['Position_Error'] and int(row['Session'][-2:])<=20:
  2992. # load seq objects (any seq, we're just using its objects)
  2993. MEC_seq = Sequences(join(row['Base_Directory'], 'Preprocessed_Data/Sequences',
  2994. row['File']+'_MEC_SWR_sequences.txt'),
  2995. load_objects=True)
  2996. CA1_seq = Sequences(join(row['Base_Directory'], 'Preprocessed_Data/Sequences',
  2997. row['File']+'_CA1_SWR_sequences.txt'),
  2998. load_objects=True)
  3000. animal_idx = animal_dict[row['Animal']]
  3001. sess_idx = int(row['Session'][-2:])-1
  3003. # collect all CA1 decoded positions for shuffling
  3004. all_SWR_CA1_decoded_pos = np.array([])
  3005. for i, interval in CA1_seq.intervals.intervals.iterrows():
  3006. cr = CA1_seq.classifier.classifier_results.sel(time=slice(interval.start_time, interval.end_time))
  3007. all_SWR_CA1_decoded_pos = np.append(all_SWR_CA1_decoded_pos, cr.sum("state").acausal_posterior.argmax("position").values)
  3009. # calculate distance between decoded positions during each interval
  3010. loc_dist_by_seq, loc_dist_by_seq_sig = \
  3011. [np.full(len(CA1_seq.intervals.intervals), np.NaN) for _ in range(2)]
  3012. n_samps = len(all_SWR_CA1_decoded_pos)
  3013. for i, interval in CA1_seq.intervals.intervals.iterrows():
  3014. cr = MEC_seq.classifier.classifier_results.sel(time=slice(interval.start_time, interval.end_time))
  3015. MEC_decoded_pos = cr.sum("state").acausal_posterior.argmax("position").values
  3016. cr = CA1_seq.classifier.classifier_results.sel(time=slice(interval.start_time, interval.end_time))
  3017. CA1_decoded_pos = cr.sum("state").acausal_posterior.argmax("position").values
  3018. loc_dist_by_seq[i] = np.nanmean(np.abs(MEC_decoded_pos-CA1_decoded_pos))
  3020. # circularly shuffle decoded positions
  3021. loc_dist_shuffle = np.full(n_shuffles, np.NaN)
  3022. for s in range(n_shuffles):
  3023. CA1_decoded_pos_shuffle = np.roll(all_SWR_CA1_decoded_pos, np.random.randint(n_samps))
  3024. loc_dist_shuffle[s] = np.nanmean(np.abs(MEC_decoded_pos-CA1_decoded_pos_shuffle[0:len(MEC_decoded_pos)]))
  3025. loc_dist_by_seq_sig[i] = np.nanmean(loc_dist_by_seq[i]<loc_dist_shuffle)
  3026. loc_dist[animal_idx][sess_idx] = np.nanmean(loc_dist_by_seq)
  3027. loc_dist_vs_shuffle[animal_idx][sess_idx] = np.nanmean(loc_dist_by_seq_sig<0.05)
  3029. print(f"{row['File']} {time.time()-start_time}")
  3031. # %%
  3032. with open(join(figure_path,'SWR_loc_dist.pkl'), 'wb') as file:
  3033. pickle.dump(loc_dist, file)
  3034. with open(join(figure_path,'SWR_loc_dist_vs_shuffle.pkl'), 'wb') as file:
  3035. pickle.dump(loc_dist_vs_shuffle, file)
  3037. # %%
  3038. # plot: KDE of distance between decodes
  3039. printf(f"Decodes in MEC and CA1 were {np.nanmean(loc_dist):.02f} +/- {stats.sem(np.nanmean(loc_dist,axis=1), nan_policy='omit'):.02f}cm apart from each other during SWRs, " +\
  3040. f"and {np.nanmean(loc_dist_vs_shuffle)*100:.02f} +/- {stats.sem(np.nanmean(loc_dist_vs_shuffle,axis=1), nan_policy='omit'):.02f}% " +\
  3041. "of SWRs had distances below shuffle.")
  3043. plot_seq_metrics_over_sessions(n_animals, loc_dist, figure_path,
  3044. 'MEC_vs_CA1_decode_dist_during_SWRs', 'Distance between MEC and CA1 Decoded Positions (cm)', [0,130])
  3045. plot_seq_metrics_over_sessions(n_animals, loc_dist_vs_shuffle*100, figure_path,
  3046. 'MEC_vs_CA1_decode_dist_during_SWRs', 'Percent of Distances < 95th Percentile', [0,100])
  3048. # %% [markdown]
  3049. # # Figure S6
  3051. # %% [markdown]
  3052. # ### Repeat with HSEs
  3054. # %%
  3055. HSE_data_by_session, HSE_data_by_seq = combine_metrics_over_sessions(sessions, animal_dict, 'CA1_HSE')
  3057. # %%
  3058. # averaged over sessions
  3059. sessions = pd.read_csv('Z:/WT_Sequences/all_sessions_Xmaze_CA1.csv')
  3060. n_sessions = 20
  3061. n_animals = len(animal_dict)
  3062. hse_overlap = np.zeros((n_animals, n_sessions))
  3064. for _, row in sessions.iterrows():
  3065. if not row['Recording_Error'] and not row['Position_Error'] and \
  3066. row['Task'] == 'X Maze' and int(row['Session'][-2:])<=20:
  3068. animal_idx = animal_dict[row['Animal']]
  3069. sess_idx = int(row['Session'][-2:])-1
  3071. # load intervals
  3072. nonlocal_immobility_intervals = Sequence_Intervals(join(row['Base_Directory'],
  3073. 'Preprocessed_Data/Sequences/Intervals',
  3074. row['File']+'_MEC_nonlocal_immobility_intervals.txt'))
  3075. HSE_intervals = Sequence_Intervals(join(row['Base_Directory'],
  3076. 'Preprocessed_Data/Sequences/Intervals',
  3077. row['File']+'_CA1_HSE_intervals.txt'))
  3079. # iterate over nonlocal intervals, then SWRs
  3080. for nl_immob in nonlocal_immobility_intervals.intervals.itertuples():
  3081. # find overlaps
  3082. for ripple in HSE_intervals.intervals.itertuples():
  3083. if ripple.start_time<nl_immob.end_time and ripple.end_time>nl_immob.start_time:
  3084. if ripple.start_time<=nl_immob.start_time:
  3085. hse_overlap[animal_idx,sess_idx] += ripple.end_time-nl_immob.start_time
  3086. elif ripple.end_time<=nl_immob.end_time:
  3087. hse_overlap[animal_idx,sess_idx] += ripple.end_time-ripple.start_time
  3088. else:
  3089. hse_overlap[animal_idx,sess_idx] += nl_immob.end_time-ripple.start_time
  3090. elif ripple.end_time>end_time:
  3091. break
  3092. hse_overlap[animal_idx,sess_idx] /= np.sum(nonlocal_immobility_intervals.intervals.end_time -\
  3093. nonlocal_immobility_intervals.intervals.start_time)
  3094. # remove sessions that don't have SWRs (e.g. no CA1 data)
  3095. hse_overlap = np.where(hse_overlap==0, np.nan, hse_overlap)
  3096. hse_overlap *= 100
  3098. # %%
  3099. print(f"{np.nanmean(np.nanmean(HSE_data_by_session['non_local'])):.02f} +/- {stats.sem(np.nanmean(HSE_data_by_session['non_local'],axis=1), nan_policy='omit'):.02f}" +\
  3100. f"% of HSEs have nonlocal coding in MEC,\n" +\
  3101. f"and {np.nanmean(np.nanmean(hse_overlap)):.02f} +/- {stats.sem(np.nanmean(hse_overlap,axis=1), nan_policy='omit'):.02f}" +\
  3102. f"% of nonlocal decoded time overlaps with HSEs.")
  3104. # %%
  3105. # plot
  3106. plot_seq_metrics_over_sessions(n_animals, hse_overlap, figure_path,
  3107. 'nl_overlap_with_hse', '% NonLocal Decode Overlapping with HSEs', [0,100])
  3108. plot_seq_metrics_over_sessions(n_animals, HSE_data_by_session['non_local'], figure_path,
  3109. 'hse_overlap_with_nl', '% HSEs with NonLocal Decode', [0,100])
  3111. # %% [markdown]
  3112. # ### compare to shuffle: % nonlocal coding overlapping with HSEs and vice versa
  3114. # %%
  3115. n_shuffles = 1000
  3116. rng = np.random.default_rng()
  3118. n_sessions = 20
  3119. n_animals = len(animal_dict)
  3120. nonlocal_contain_HSEs_5, nonlocal_contain_HSEs_95, nonlocal_contain_HSEs_median, \
  3121. HSEs_contain_nonlocal_5, HSEs_contain_nonlocal_95, HSEs_contain_nonlocal_median =\
  3122. [np.full((n_animals, n_sessions), np.NaN) for _ in range(6)]
  3124. sessions = pd.read_csv('Z:/WT_Sequences/all_sessions_Xmaze_CA1.csv')
  3125. for _, row in sessions.iterrows():
  3126. if not row['Recording_Error'] and not row['Position_Error'] and \
  3127. row['Task'] == 'X Maze' and int(row['Session'][-2:])<=20:
  3129. animal_idx = animal_dict[row['Animal']]
  3130. sess_idx = int(row['Session'][-2:])-1
  3132. seq = Sequences(join(row['Base_Directory'],
  3133. 'Preprocessed_Data/Sequences',
  3134. row['File']+'_MEC_immobility_sequences.txt'),
  3135. load_objects=True)
  3136. dist = seq.calc_dist_from_decode()
  3137. ts = seq.binned_data.timestamps
  3138. SWR_intervals = Sequence_Intervals(join(row['Base_Directory'],
  3139. 'Preprocessed_Data/Sequences/Intervals',
  3140. row['File']+'_CA1_HSE_intervals.txt'))
  3142. # 1000 shuffles
  3143. perc_SWRs_contain_nonlocal = np.zeros(1000)
  3144. perc_nonlocal_contain_SWRs = np.zeros(1000)
  3145. for b in range(n_shuffles):
  3146. ts_nonlocal = np.array([])
  3147. for ripple in SWR_intervals.indices.itertuples():
  3148. # pick a random immobility interval
  3149. seq_idx = int(rng.random()*len(seq.stats))
  3150. selected_interval = seq.intervals.indices.iloc[seq_idx]
  3151. # pick a random time within that interval of the same duration as the SWR
  3152. SWR_idx_len = ripple.end_index-ripple.start_index
  3153. index_range = selected_interval.end_index-selected_interval.start_index-SWR_idx_len
  3154. shuffled_start_idx = int(rng.random()*index_range)
  3155. ts_subset = ts[shuffled_start_idx:shuffled_start_idx+SWR_idx_len]
  3156. is_nonlocal = dist[shuffled_start_idx:shuffled_start_idx+SWR_idx_len]>20
  3157. count_nonlocal = np.nansum(is_nonlocal)
  3158. ts_nonlocal = np.concatenate((ts_nonlocal, ts_subset[is_nonlocal]))
  3159. if count_nonlocal>0:
  3160. perc_SWRs_contain_nonlocal[b] += 1
  3161. perc_nonlocal_contain_SWRs[b] = len(ts_nonlocal) #len(np.unique(ts_nonlocal))
  3162. perc_SWRs_contain_nonlocal = perc_SWRs_contain_nonlocal/len(SWR_intervals.intervals) * 100
  3163. perc_nonlocal_contain_SWRs = perc_nonlocal_contain_SWRs / \
  3164. (np.sum(nonlocal_immobility_intervals.indices.end_index -\
  3165. nonlocal_immobility_intervals.indices.start_index)) * 100
  3167. HSEs_contain_nonlocal_5[animal_idx, sess_idx] = np.percentile(perc_SWRs_contain_nonlocal, 5)
  3168. HSEs_contain_nonlocal_95[animal_idx, sess_idx] = np.percentile(perc_SWRs_contain_nonlocal, 95)
  3169. HSEs_contain_nonlocal_median[animal_idx, sess_idx] = np.percentile(perc_SWRs_contain_nonlocal, 50)
  3170. nonlocal_contain_HSEs_5[animal_idx, sess_idx] = np.percentile(perc_nonlocal_contain_SWRs, 5)
  3171. nonlocal_contain_HSEs_95[animal_idx, sess_idx] = np.percentile(perc_nonlocal_contain_SWRs, 95)
  3172. nonlocal_contain_HSEs_median[animal_idx, sess_idx] = np.percentile(perc_nonlocal_contain_SWRs, 50)
  3174. # %%
  3175. # % NonLocal Decode Overlapping with SWRs
  3176. n_nonnan = np.count_nonzero(~np.isnan(hse_overlap))
  3177. n_more_than_95 = np.nansum(nonlocal_contain_HSEs_95 < hse_overlap)
  3178. print(f"{n_more_than_95} of {n_nonnan} sessions ({n_more_than_95/n_nonnan*100:0.2f}%) had HSEs during non-local content than if SWRs were randomly allocated across immobility.")
  3180. # % SWRs with NonLocal Decode
  3181. n_nonnan = np.count_nonzero(~np.isnan(HSE_data_by_session['non_local']))
  3182. n_more_than_95 = np.nansum(HSEs_contain_nonlocal_95 < HSE_data_by_session['non_local'])
  3183. print(f"{n_more_than_95} of {n_nonnan} sessions ({n_more_than_95/n_nonnan*100:0.2f}%) had more non-local content during HSEs than if SWRs were randomly allocated across immobility.")
  3185. # %%
  3186. f, ax = plt.subplots(1, figsize=(4,4))
  3187. sns.kdeplot(nonlocal_contain_HSEs_95.flatten(), bw_adjust=0.75, color='red')
  3188. ax.set(xlabel='% HSEs with NonLocal Decode', ylabel='Proportion', xlim=[0,100])
  3189. f.savefig(join(figure_path,'hse_overlap_with_nl_over_events_shuffle.pdf'), format='pdf', transparent=True,
  3190. dpi=300, bbox_inches='tight')
  3192. f, ax = plt.subplots(1, figsize=(4,4))
  3193. sns.kdeplot(HSEs_contain_nonlocal_95.flatten(), bw_adjust=0.75, color='red')
  3194. ax.set(xlabel='% NonLocal Decode Overlapping with HSEs', ylabel='Proportion', xlim=[0,100])
  3195. f.savefig(join(figure_path,'nl_overlap_with_hse_over_events_shuffle.pdf'), format='pdf', transparent=True,
  3196. dpi=300, bbox_inches='tight')
  3198. # %% [markdown]
  3199. # ### Non-local content that overlaps with SWRs vs not
  3201. # %%
  3202. sessions = pd.read_csv('Z:/WT_Sequences/all_sessions_Xmaze_CA1.csv')
  3204. animal = []
  3205. session = []
  3206. nl_ripple_overlap = []
  3207. seq_stats = None
  3209. for _, row in sessions.iterrows():
  3210. if not row['Recording_Error'] and not row['Position_Error'] and \
  3211. row['Task'] == 'X Maze' and int(row['Session'][-2:])<=20:
  3213. animal_idx = animal_dict[row['Animal']]
  3214. sess_idx = int(row['Session'][-2:])-1
  3216. # load intervals
  3217. seq = Sequences(join(row['Base_Directory'],
  3218. 'Preprocessed_Data/Sequences',
  3219. row['File']+'_MEC_nonlocal_immobility_sequences.txt'))
  3220. intervals = Sequence_Intervals(join(row['Base_Directory'],
  3221. 'Preprocessed_Data/Sequences/Intervals',
  3222. row['File']+'_MEC_nonlocal_immobility_intervals.txt'))
  3223. SWR_intervals = Sequence_Intervals(join(row['Base_Directory'],
  3224. 'Preprocessed_Data/Sequences/Intervals',
  3225. row['File']+'_CA1_SWR_intervals.txt'))
  3226. nl_ripple_overlap_session = np.full(len(seq.stats), False)
  3228. # iterate over nonlocal intervals, then SWRs
  3229. if len(SWR_intervals.intervals>0):
  3230. for nl in range(len(seq.stats)):
  3231. # find overlaps
  3232. for ripple in SWR_intervals.intervals.itertuples():
  3233. if ripple.start_time<intervals.intervals.iloc[nl].end_time and \
  3234. ripple.end_time>intervals.intervals.iloc[nl].start_time:
  3235. nl_ripple_overlap_session[nl] = True
  3236. break
  3237. elif ripple.end_time>end_time:
  3238. break
  3240. # build arrays over all sequences
  3241. animal_idx = animal_dict[row['Animal']]
  3242. animal.extend([animal_idx] * len(seq.stats))
  3243. session_idx = int(row['Session'][-2:])-1
  3244. session.extend([session_idx] * len(seq.stats))
  3245. nl_ripple_overlap.extend(nl_ripple_overlap_session)
  3246. if seq_stats is not None:
  3247. seq_stats = pd.concat([seq_stats, seq.stats], axis=0, ignore_index=True)
  3248. else:
  3249. seq_stats = seq.stats
  3251. # build full df
  3252. ripple_overlap_by_seq = pd.DataFrame({'Animal': animal, 'Session': session, 'Ripple_overlap': nl_ripple_overlap})
  3253. ripple_overlap_by_seq = pd.concat([ripple_overlap_by_seq, seq_stats], axis=1)
  3255. # %%
  3256. # match durations
  3257. ripple_equivalent = []
  3258. yesripple_data_by_seq = ripple_overlap_by_seq.loc[ripple_overlap_by_seq.Ripple_overlap]
  3259. noripple_data_by_seq = ripple_overlap_by_seq.loc[~ripple_overlap_by_seq.Ripple_overlap]
  3261. for a in noripple_data_by_seq.Animal.unique():
  3262. start_idx = np.min(np.where(yesripple_data_by_seq.Animal==a))
  3263. for __, dur in noripple_data_by_seq.loc[noripple_data_by_seq.Animal==a].Duration.items():
  3264. idx = np.abs(yesripple_data_by_seq.loc[yesripple_data_by_seq.Animal==a].Duration - dur).argmin()
  3265. ripple_equivalent.append(start_idx + idx)
  3267. # %%
  3268. nl_ripple_overlap = []
  3269. nl_ripple_overlap.extend([False] * len(noripple_data_by_seq))
  3270. nl_ripple_overlap.extend([True] * len(yesripple_data_by_seq.iloc[ripple_equivalent]))
  3271. ripple_overlap_equivalent_by_seq = pd.concat([noripple_data_by_seq, yesripple_data_by_seq.iloc[ripple_equivalent]], ignore_index=True)
  3272. ripple_overlap_equivalent_by_seq.Ripple_overlap = nl_ripple_overlap
  3274. plot_seq_metrics_comparison_over_events(ripple_overlap_equivalent_by_seq, 'FR', 'Ripple_overlap',
  3275. figure_path, 'Nonlocal_SWR_vs_not_fr_matched', 'FR (Hz)', [0,50])
  3276. plot_seq_metrics_comparison_over_events(ripple_overlap_equivalent_by_seq, 'Percent_units', 'Ripple_overlap',
  3277. figure_path, 'Nonlocal_SWR_vs_not_percent_units_matched', '% Active Units', [0,100])
  3278. plot_seq_metrics_comparison_over_events(ripple_overlap_equivalent_by_seq, 'Spatial_information', 'Ripple_overlap',
  3279. figure_path, 'Nonlocal_SWR_vs_not_spatial_info_matched', 'Spatial information (bit/s)', [0,1.5])
  3280. plot_seq_metrics_comparison_over_events(ripple_overlap_equivalent_by_seq, 'Duration', 'Ripple_overlap',
  3281. figure_path, 'Nonlocal_SWR_vs_not_duration_matched', 'Duration (s)', [0,1])
  3282. plot_seq_metrics_comparison_over_events(ripple_overlap_equivalent_by_seq, 'Max_distance_from_animal', 'Ripple_overlap',
  3283. figure_path, 'Nonlocal_SWR_vs_not_lookahead_matched', 'Lookahead distance (cm)', [0,130])
  3284. plot_seq_metrics_comparison_over_events(ripple_overlap_equivalent_by_seq, 'Spatial_coverage_percent', 'Ripple_overlap',
  3285. figure_path, 'Nonlocal_SWR_vs_not_spatial_coverage_matched', '% Spatial coverage', [0,100])
  3287. # %% [markdown]
  3288. # ### SWRs that contain nonlocal content in MEC vs not
  3290. # %%
  3291. # match durations
  3292. nonlocal_SWR_equivalent = []
  3293. local_SWR_data_by_seq = SWR_data_by_seq.loc[~SWR_data_by_seq.Non_local]
  3294. nonlocal_SWR_data_by_seq = SWR_data_by_seq.loc[SWR_data_by_seq.Non_local]
  3296. for a in local_SWR_data_by_seq.Animal.unique():
  3297. start_idx = np.min(np.where(nonlocal_SWR_data_by_seq.Animal==a))
  3298. for __, dur in local_SWR_data_by_seq.loc[local_SWR_data_by_seq.Animal==a].Duration.items():
  3299. idx = np.abs(nonlocal_SWR_data_by_seq.loc[nonlocal_SWR_data_by_seq.Animal==a].Duration - dur).argmin()
  3300. nonlocal_SWR_equivalent.append(start_idx + idx)
  3302. # %%
  3303. non_local = []
  3304. non_local.extend([False] * len(local_SWR_data_by_seq))
  3305. non_local.extend([True] * len(nonlocal_SWR_data_by_seq.iloc[nonlocal_SWR_equivalent]))
  3306. SWR_data_equivalent_by_seq = pd.concat([local_SWR_data_by_seq, nonlocal_SWR_data_by_seq.iloc[nonlocal_SWR_equivalent]], ignore_index=True)
  3307. SWR_data_equivalent_by_seq.Non_local = non_local
  3309. plot_seq_metrics_comparison_over_events(SWR_data_equivalent_by_seq, 'FR', 'Non_local',
  3310. figure_path, 'SWR_local_vs_nonlocal_fr_matched', 'FR (Hz)', [0,50])
  3311. plot_seq_metrics_comparison_over_events(SWR_data_equivalent_by_seq, 'Percent_units', 'Non_local',
  3312. figure_path, 'SWR_local_vs_nonlocal_percent_units_matched', '% Active Units', [0,20])
  3313. plot_seq_metrics_comparison_over_events(SWR_data_equivalent_by_seq, 'Spatial_information', 'Non_local',
  3314. figure_path, 'SWR_local_vs_nonlocal_spatial_info_matched', 'Spatial information (bit/s)', [0,1.5])
  3315. plot_seq_metrics_comparison_over_events(SWR_data_equivalent_by_seq, 'Duration', 'Non_local',
  3316. figure_path, 'SWR_local_vs_nonlocal_duration_matched', 'Duration (s)', [0,0.2])
  3317. plot_seq_metrics_comparison_over_events(SWR_data_equivalent_by_seq, 'Animal_spatial_bin', 'Non_local',
  3318. figure_path, 'SWR_local_vs_nonlocal_animal_location_matched', 'Animal linearized position (cm)', [0,123])
  3320. # %% [markdown]
  3321. # # Figure S6 - part 2
  3323. # %% [markdown]
  3324. # ### example raw & filtered traces
  3326. # %%
  3327. sessions = pd.read_csv(r'Z:/WT_Sequences/all_sessions_Xmaze.csv')
  3328. row = sessions.iloc[243]
  3330. nonlocal_intervals = Sequence_Intervals(join(row['Base_Directory'], 'Preprocessed_Data/Sequences/Intervals',
  3331. row['File']+'_MEC_nonlocal_immobility_intervals.txt'))
  3332. local_intervals = Sequence_Intervals(join(row['Base_Directory'], 'Preprocessed_Data/Sequences/Intervals',
  3333. row['File']+'_MEC_local_immobility_intervals.txt'))
  3334. move_intervals = Sequence_Intervals(join(row['Base_Directory'], 'Preprocessed_Data/Sequences/Intervals',
  3335. row['File']+'_movement_intervals.txt'))
  3336. lfp = LFP(join(row['Base_Directory'], 'Preprocessed_Data/LFP', row['File']+'_MEC_theta_channel_lfp.txt'))
  3338. lfp_move = deepcopy(lfp)
  3339. lfp_move.subset_by_time(move_intervals.intervals.values)
  3341. lfp_nonlocal = deepcopy(lfp)
  3342. lfp_nonlocal.subset_by_time(nonlocal_intervals.intervals.values)
  3344. lfp_local = deepcopy(lfp)
  3345. lfp_local.subset_by_time(local_intervals.intervals.values)
  3347. # %%
  3348. print(lfp_move.name)
  3349. f, ax = plt.subplots(3)
  3350. high_theta_lfp = lf.bandpass_filter(lfp_local.lfp, low=lf.HIGH_THETA[0], high=lf.HIGH_THETA[1], pass_band=lf.HIGH_THETA[2])
  3351. low_theta_lfp = lf.bandpass_filter(lfp_local.lfp, low=lf.LOW_THETA[0], high=lf.LOW_THETA[1], pass_band=lf.LOW_THETA[2])
  3352. ax[0].plot(low_theta_lfp.squeeze(), color='b')
  3353. ax[0].plot(high_theta_lfp.squeeze(), color='r')
  3354. ax[0].plot(lfp_local.lfp)
  3355. ax[0].set(xlim=[1958,1958+625], ylim=[-150,150])
  3357. high_theta_lfp = lf.bandpass_filter(lfp_nonlocal.lfp, low=lf.HIGH_THETA[0], high=lf.HIGH_THETA[1], pass_band=lf.HIGH_THETA[2])
  3358. low_theta_lfp = lf.bandpass_filter(lfp_nonlocal.lfp, low=lf.LOW_THETA[0], high=lf.LOW_THETA[1], pass_band=lf.LOW_THETA[2])
  3359. ax[1].plot(low_theta_lfp.squeeze(), color='b')
  3360. ax[1].plot(high_theta_lfp.squeeze(), color='r')
  3361. ax[1].plot(lfp_nonlocal.lfp)
  3362. ax[1].set(xlim=[37798,37798+625], ylim=[-150,150])
  3364. high_theta_lfp = lf.bandpass_filter(lfp_move.lfp, low=lf.HIGH_THETA[0], high=lf.HIGH_THETA[1], pass_band=lf.HIGH_THETA[2])
  3365. low_theta_lfp = lf.bandpass_filter(lfp_move.lfp, low=lf.LOW_THETA[0], high=lf.LOW_THETA[1], pass_band=lf.LOW_THETA[2])
  3366. ax[2].plot(low_theta_lfp.squeeze(), color='b')
  3367. ax[2].plot(high_theta_lfp.squeeze(), color='r')
  3368. ax[2].plot(lfp_move.lfp)
  3369. ax[2].set(xlim=[625,2*625], ylim=[-150,150])
  3370. f.savefig(join(figure_path,'TopHat_DY06_theta_traces.pdf'), format='pdf', transparent=True,
  3371. dpi=300, bbox_inches='tight')
  3373. # %% [markdown]
  3374. # ### PSD in theta range
  3376. # %%
  3377. # PSD over movement vs immobiilty times
  3378. sessions = pd.read_csv(r'Z:/WT_Sequences/all_sessions_Xmaze.csv')
  3379. n_sessions = 20
  3380. n_animals = len(animal_dict)
  3381. rate = 625
  3382. n_freqs = 129
  3383. immob_psd = np.full((n_animals, n_sessions, n_freqs), np.NaN)
  3384. move_psd = np.full((n_animals, n_sessions, n_freqs), np.NaN)
  3386. start_time = time.time()
  3387. for _, row in sessions.iterrows():
  3388. if not row['Recording_Error'] and not row['Position_Error'] and int(row['Session'][-2:])<=20:
  3389. animal_idx = animal_dict[row['Animal']]
  3390. sess_idx = int(row['Session'][-2:])-1
  3392. immob_intervals = Sequence_Intervals(join(row['Base_Directory'], 'Preprocessed_Data/Sequences/Intervals',
  3393. row['File']+'_immobility_intervals.txt'))
  3394. move_intervals = Sequence_Intervals(join(row['Base_Directory'], 'Preprocessed_Data/Sequences/Intervals',
  3395. row['File']+'_movement_intervals.txt'))
  3396. lfp = LFP(join(row['Base_Directory'], 'Preprocessed_Data/LFP', row['File']+'_MEC_theta_channel_lfp.txt'))
  3398. lfp_immob = deepcopy(lfp)
  3399. lfp_immob.subset_by_time(immob_intervals.intervals.values)
  3400. freqs, immob_psd[animal_idx, sess_idx, :] = signal.welch(np.squeeze(lfp_immob.lfp), rate)
  3402. lfp_move = deepcopy(lfp)
  3403. lfp_move.subset_by_time(move_intervals.intervals.values)
  3404. _, move_psd[animal_idx, sess_idx, :] = signal.welch(np.squeeze(lfp_move.lfp), rate)
  3405. print(f"{row['File']} {time.time()-start_time}")
  3407. # %%
  3408. np.save(join(figure_path,'immob_psd.npy'), immob_psd)
  3409. np.save(join(figure_path,'move_psd.npy'), move_psd)
  3410. nonlocal_psd = np.load(join(figure_path,'nonlocal_psd.npy'))
  3411. local_psd = np.load(join(figure_path,'local_psd.npy'))
  3412. for a in range(n_animals):
  3413. for s in range(n_sessions):
  3414. immob_psd[a][s] /= np.nansum(immob_psd[a][s])
  3415. move_psd[a][s] /= np.nansum(move_psd[a][s])
  3416. nonlocal_psd[a][s] /= np.nansum(nonlocal_psd[a][s])
  3417. local_psd[a][s] /= np.nansum(local_psd[a][s])
  3418. freqs = freqs + (freqs[1]-freqs[0])/2 # plot at bin centers, not edges
  3420. # %%
  3421. # sum & plot over sessions
  3422. immob_psd_summed = np.nanmean(immob_psd, axis=1)
  3423. move_psd_summed = np.nanmean(move_psd, axis=1)
  3424. nonlocal_psd_summed = np.nanmean(nonlocal_psd, axis=1)
  3425. local_psd_summed = np.nanmean(local_psd, axis=1)
  3427. # comparison: mean + sem over animals
  3428. f, ax = plt.subplots()
  3429. ax.plot(freqs, np.nanmean(local_psd_summed, axis=0), color='black')
  3430. sem = stats.sem(local_psd_summed, nan_policy='omit')
  3431. ax.fill_between(freqs, np.nanmean(local_psd_summed, axis=0)-sem, \
  3432. np.nanmean(local_psd_summed, axis=0)+sem, color='grey', alpha=0.5, rasterized=True)
  3433. ax.plot(freqs, np.nanmean(nonlocal_psd_summed, axis=0), color='#1C968B')
  3434. sem = stats.sem(nonlocal_psd_summed, nan_policy='omit')
  3435. ax.fill_between(freqs, np.nanmean(nonlocal_psd_summed, axis=0)-sem, \
  3436. np.nanmean(nonlocal_psd_summed, axis=0)+sem, color='#1C968B', alpha=0.5, rasterized=True)
  3437. ax.plot(freqs, np.nanmean(move_psd_summed, axis=0), color='magenta')
  3438. sem = stats.sem(move_psd_summed, nan_policy='omit')
  3439. ax.fill_between(freqs, np.nanmean(move_psd_summed, axis=0)-sem, \
  3440. np.nanmean(move_psd_summed, axis=0)+sem, color='magenta', alpha=0.5, rasterized=True)
  3441. ax.set(xlabel='Frequency', ylabel='PSD')
  3442. f.savefig(join(figure_path,'local_nonlocal_move_psd.pdf'), format='pdf')
  3443. ax.set(xlim=[0,20])
  3444. f.savefig(join(figure_path,'local_nonlocal_move_psd_inset.pdf'), format='pdf')
  3446. # %% [markdown]
  3447. # ### Instantaneous frequency
  3449. # %%
  3450. # inst freq
  3451. sessions = pd.read_csv(r'Z:/WT_Sequences/all_sessions_Xmaze.csv')
  3452. n_sessions = 20
  3453. n_animals = len(animal_dict)
  3454. nonlocal_freqs, local_freqs, move_freqs = \
  3455. [np.full((n_animals, n_sessions), np.NaN) for _ in range(3)]
  3457. start_time = time.time()
  3458. for _, row in sessions.iterrows():
  3459. if not row['Recording_Error'] and not row['Position_Error'] and int(row['Session'][-2:])<=20:
  3460. animal_idx = animal_dict[row['Animal']]
  3461. sess_idx = int(row['Session'][-2:])-1
  3463. nonlocal_intervals = Sequence_Intervals(join(row['Base_Directory'], 'Preprocessed_Data/Sequences/Intervals',
  3464. row['File']+'_MEC_nonlocal_immobility_intervals.txt'))
  3465. local_intervals = Sequence_Intervals(join(row['Base_Directory'], 'Preprocessed_Data/Sequences/Intervals',
  3466. row['File']+'_MEC_local_immobility_intervals.txt'))
  3467. move_intervals = Sequence_Intervals(join(row['Base_Directory'], 'Preprocessed_Data/Sequences/Intervals',
  3468. row['File']+'_movement_intervals.txt'))
  3469. lfp = LFP(join(row['Base_Directory'], 'Preprocessed_Data/LFP', row['File']+'_MEC_theta_channel_lfp.txt'))
  3471. lfp_nonlocal = deepcopy(lfp)
  3472. lfp_nonlocal.subset_by_time(nonlocal_intervals.intervals.values)
  3473. _, _, inst_freqs = lf.hilbert_envelope_phase_freq(lf.bandpass_filter(lfp_nonlocal.lfp, low=lf.THETA[0], high=lf.THETA[1], pass_band=lf.THETA[2]))
  3474. inst_freqs = inst_freqs.squeeze()
  3475. nonlocal_freqs[animal_idx, sess_idx] = np.mean(inst_freqs[inst_freqs<20])
  3477. lfp_local = deepcopy(lfp)
  3478. lfp_local.subset_by_time(local_intervals.intervals.values)
  3479. _, _, inst_freqs = lf.hilbert_envelope_phase_freq(lf.bandpass_filter(lfp_local.lfp, low=lf.THETA[0], high=lf.THETA[1], pass_band=lf.THETA[2]))
  3480. inst_freqs = inst_freqs.squeeze()
  3481. local_freqs[animal_idx, sess_idx] = np.mean(inst_freqs[inst_freqs<20])
  3483. lfp_move = deepcopy(lfp)
  3484. lfp_move.subset_by_time(move_intervals.intervals.values)
  3485. _, _, inst_freqs = lf.hilbert_envelope_phase_freq(lf.bandpass_filter(lfp_move.lfp, low=lf.THETA[0], high=lf.THETA[1], pass_band=lf.THETA[2]))
  3486. inst_freqs = inst_freqs.squeeze()
  3487. move_freqs[animal_idx, sess_idx] = np.mean(inst_freqs[inst_freqs<20])
  3489. print(f"{row['File']} {time.time()-start_time}")
  3491. # %%
  3492. delta_freqs = local_freqs - nonlocal_freqs
  3493. print(f"Instantaneous frequency decreases by {np.mean(np.nanmean(delta_freqs, axis=1)):.02f} +/- {stats.sem(np.nanmean(delta_freqs, axis=1)):.02f} Hz")
  3494. print(f"Instantaneous frequencies: {np.mean(np.nanmean(local_freqs, axis=1)):.02f} +/- {stats.sem(np.nanmean(local_freqs, axis=1)):.02f} Hz\n" + \
  3495. f"{np.mean(np.nanmean(nonlocal_freqs, axis=1)):.02f} +/- {stats.sem(np.nanmean(nonlocal_freqs, axis=1)):.02f} Hz\n" + \
  3496. f"{np.mean(np.nanmean(move_freqs, axis=1)):.02f} +/- {stats.sem(np.nanmean(move_freqs, axis=1)):.02f} Hz")
  3498. # %%
  3499. data_by_session = pd.DataFrame({'Animal': np.repeat(np.arange(n_animals),n_sessions),
  3500. 'Session': np.tile(np.arange(n_sessions),n_animals),
  3501. 'Feature_l': local_freqs.flatten(),
  3502. 'Feature_nl': nonlocal_freqs.flatten()})
  3503. plot_compare_paired_metrics('Feature_l', 'Feature_nl', data_by_session, \
  3504. figure_path, 'theta_inst_freq_l_vs_nl', 'Instantaneous Frequency (Hz)', [6,9])
  3505. data_by_session = pd.DataFrame({'Animal': np.repeat(np.arange(n_animals),n_sessions),
  3506. 'Session': np.tile(np.arange(n_sessions),n_animals),
  3507. 'Feature_l': nonlocal_freqs.flatten(),
  3508. 'Feature_nl': move_freqs.flatten()})
  3509. plot_compare_paired_metrics('Feature_l', 'Feature_nl', data_by_session, \
  3510. figure_path, 'theta_inst_freq_nl_vs_move', 'Instantaneous Frequency (Hz)', [6,9])
  3512. # %% [markdown]
  3513. # ### Power in theta bands
  3515. # %%
  3516. # Z-score of low and high theta bands
  3517. sessions = pd.read_csv(r'Z:/WT_Sequences/all_sessions_Xmaze.csv')
  3518. n_sessions = 20
  3519. n_animals = len(animal_dict)
  3520. window = 0.5
  3521. nl_oldmethod_theta, nl_newmethod_theta, l_oldmethod_theta, l_newmethod_theta = \
  3522. [np.full((n_animals, n_sessions), np.NaN) for _ in range(4)]
  3524. start_time = time.time()
  3525. for _, row in sessions.iterrows():
  3526. if not row['Recording_Error'] and not row['Position_Error'] and int(row['Session'][-2:])<=20:
  3527. animal_idx = animal_dict[row['Animal']]
  3528. sess_idx = int(row['Session'][-2:])-1
  3530. nonlocal_intervals = Sequence_Intervals(join(row['Base_Directory'], 'Preprocessed_Data/Sequences/Intervals',
  3531. row['File']+'_MEC_nonlocal_immobility_intervals.txt'))
  3532. local_intervals = Sequence_Intervals(join(row['Base_Directory'], 'Preprocessed_Data/Sequences/Intervals',
  3533. row['File']+'_MEC_local_immobility_intervals.txt'))
  3534. immob_intervals = Sequence_Intervals(join(row['Base_Directory'], 'Preprocessed_Data/Sequences/Intervals',
  3535. row['File']+'_immobility_intervals.txt'))
  3536. lfp = LFP(join(row['Base_Directory'], 'Preprocessed_Data/LFP', row['File']+'_MEC_theta_channel_lfp.txt'))
  3538. lfp_nonlocal = deepcopy(lfp)
  3539. lfp_nonlocal.subset_by_time(nonlocal_intervals.intervals.values)
  3540. lfp_local = deepcopy(lfp)
  3541. lfp_local.subset_by_time(local_intervals.intervals.values)
  3542. lfp_immob = deepcopy(lfp)
  3543. lfp_immob.subset_by_time(immob_intervals.intervals.values)
  3545. oldmethod_immob_theta, oldmethod_theta_mean = lf.multitaper_filtered_power(lfp_immob.lfp, low=lf.NONLOCAL[0], high=lf.NONLOCAL[1], window=1)
  3546. oldmethod_theta_std = np.nanstd(oldmethod_immob_theta)
  3547. _, nl_oldmethod_theta_power = lf.multitaper_filtered_power(lfp_nonlocal.lfp, low=lf.NONLOCAL[0], high=lf.NONLOCAL[1], window=1)
  3548. nl_oldmethod_theta[animal_idx, sess_idx] = (nl_oldmethod_theta_power-oldmethod_theta_mean)/oldmethod_theta_std
  3549. _, l_oldmethod_theta_power = lf.multitaper_filtered_power(lfp_local.lfp, low=lf.NONLOCAL[0], high=lf.NONLOCAL[1], window=1)
  3550. l_oldmethod_theta[animal_idx, sess_idx] = (l_oldmethod_theta_power-oldmethod_theta_mean)/oldmethod_theta_std
  3552. nl_newmethod_theta_power, l_newmethod_theta_power = \
  3553. [[] for _ in range(2)]
  3555. # subset on intervals that are at least the size of the window
  3556. intervals = nonlocal_intervals.intervals.values
  3557. idx = np.where((intervals[:,1]-intervals[:,0])>window)[0]
  3558. for _, i in enumerate(idx):
  3559. power, _ = lf.multitaper_filtered_power(lfp.lfp[np.where((lfp.timestamps >= intervals[i, 0]) \
  3560. & (lfp.timestamps <= intervals[i, 1]))[0]], \
  3561. low=lf.NONLOCAL[0], high=lf.NONLOCAL[1], window=window)
  3562. nl_newmethod_theta_power.extend(power)
  3564. # repeat for local
  3565. intervals = local_intervals.intervals.values
  3566. idx = np.where((intervals[:,1]-intervals[:,0])>window)[0]
  3567. for _, i in enumerate(idx):
  3568. power, _ = lf.multitaper_filtered_power(lfp.lfp[np.where((lfp.timestamps >= intervals[i, 0]) \
  3569. & (lfp.timestamps <= intervals[i, 1]))[0]], \
  3570. low=lf.NONLOCAL[0], high=lf.NONLOCAL[1], window=window)
  3571. l_newmethod_theta_power.extend(power)
  3573. # z-score and add to per-session calculation
  3574. newmethod_theta_mean = np.mean(nl_newmethod_theta_power+l_newmethod_theta_power)
  3575. newmethod_theta_std = np.std(nl_newmethod_theta_power+l_newmethod_theta_power)
  3576. nl_newmethod_theta[animal_idx, sess_idx] = (np.mean(nl_newmethod_theta_power)-theta_mean)/theta_std
  3577. l_newmethod_theta[animal_idx, sess_idx] = (np.mean(l_newmethod_theta_power)-theta_mean)/theta_std
  3579. print(f"{row['File']} {time.time()-start_time}")
  3581. # %%
  3582. print(f"Type II theta power changes by {np.mean(np.nanmean(nl_oldmethod_theta-l_oldmethod_theta, axis=1)):.02f}" +\
  3583. f" +/- {stats.sem(np.nanmean(nl_oldmethod_theta-l_oldmethod_theta, axis=1)):.04f}")
  3584. print(f"Z-score Power: {np.mean(np.nanmean(l_oldmethod_theta, axis=1)):.02f} +/- {stats.sem(np.nanmean(l_oldmethod_theta, axis=1)):.04f}\n" + \
  3585. f"{np.mean(np.nanmean(nl_oldmethod_theta, axis=1)):.02f} +/- {stats.sem(np.nanmean(nl_oldmethod_theta, axis=1)):.04f}")
  3587. # %%
  3588. data_by_session = pd.DataFrame({'Animal': np.repeat(np.arange(n_animals),n_sessions),
  3589. 'Session': np.tile(np.arange(n_sessions),n_animals),
  3590. 'Feature_l': l_oldmethod_theta.flatten(),
  3591. 'Feature_nl': nl_oldmethod_theta.flatten()})
  3592. plot_compare_paired_metrics('Feature_l', 'Feature_nl', data_by_session, \
  3593. figure_path, 'low_theta_zscore_l_vs_nl', '3-7 Hz Power (Z-score)', [-1,1])
  3594. data_by_session = pd.DataFrame({'Animal': np.repeat(np.arange(n_animals),n_sessions),
  3595. 'Session': np.tile(np.arange(n_sessions),n_animals),
  3596. 'Feature_l': l_newmethod_theta.flatten(),
  3597. 'Feature_nl': nl_newmethod_theta.flatten()})
  3598. plot_compare_paired_metrics('Feature_l', 'Feature_nl', data_by_session, \
  3599. figure_path, 'high_theta_zscore_l_vs_nl', '3-7 Hz Power (Z-score)', [-1,1])
  3601. # %% [markdown]
  3602. # ### Spike-field coherence in theta band
  3604. # %%
  3605. # calculate
  3606. sessions = pd.read_csv(r'Z:/WT_Sequences/all_sessions_Xmaze.csv')
  3607. n_sessions = 20
  3608. n_animals = len(animal_dict)
  3609. nl_low_theta_coh, nl_high_theta_coh, l_low_theta_coh, l_high_theta_coh, \
  3610. delta_low_theta_coh, delta_high_theta_coh = \
  3611. [np.full((n_animals, n_sessions), np.NaN) for _ in range(6)]
  3613. start_time = time.time()
  3614. for _, row in sessions.iterrows():
  3615. if not row['Recording_Error'] and not row['Position_Error'] and int(row['Session'][-2:])<=20:
  3616. animal_idx = animal_dict[row['Animal']]
  3617. sess_idx = int(row['Session'][-2:])-1
  3619. # load intervals
  3620. nonlocal_intervals = Sequence_Intervals(join(row['Base_Directory'], 'Preprocessed_Data/Sequences/Intervals',
  3621. row['File']+'_MEC_nonlocal_immobility_intervals.txt'))
  3622. local_intervals = Sequence_Intervals(join(row['Base_Directory'], 'Preprocessed_Data/Sequences/Intervals',
  3623. row['File']+'_MEC_local_immobility_intervals.txt'))
  3625. # load LFP and subset on intervals
  3626. lfp = LFP(join(row['Base_Directory'], 'Preprocessed_Data/LFP', row['File']+'_MEC_theta_channel_lfp.txt'))
  3627. lfp_nonlocal = deepcopy(lfp)
  3628. lfp_nonlocal.subset_by_time(nonlocal_intervals.intervals.values)
  3629. lfp_local = deepcopy(lfp)
  3630. lfp_local.subset_by_time(local_intervals.intervals.values)
  3632. # load spikes
  3633. if ('2021_pilot' in row['Base_Directory']) or ('2022_winter' in row['Base_Directory']):
  3634. electrodes = Electrodes(join(row['Base_Directory'], 'Preprocessed_Data/Electrodes',
  3635. row['File']+'_electrodes.txt'))
  3636. spikes = Spikes(join(row['Base_Directory'], 'Preprocessed_Data/Spikes/g1', row['File']+'_imec0_spikes.txt'))
  3637. channels = electrodes.subset_by_location(regions=['Entorhinal area medial part dorsal zone'])
  3638. else:
  3639. if (row['Animal']=='Lamarr') and (sess_idx>=7):
  3640. electrodes = Electrodes(join(row['Base_Directory'], 'Preprocessed_Data/Electrodes',
  3641. row['File']+'_imec0_electrodes.txt'))
  3642. spikes = Spikes(join(row['Base_Directory'], 'Preprocessed_Data/Spikes',
  3643. row['File']+'_imec0_spikes.txt'))
  3644. else:
  3645. electrodes = Electrodes(join(row['Base_Directory'], 'Preprocessed_Data/Electrodes',
  3646. row['File']+'_imec1_electrodes.txt'))
  3647. spikes = Spikes(join(row['Base_Directory'], 'Preprocessed_Data/Spikes',
  3648. row['File']+'_imec1_spikes.txt'))
  3649. channels = electrodes.subset_by_location(regions=['ENTm1','ENTm2','ENTm3','ENTm4','ENTm5','ENTm6'])
  3650. spikes.subset_by_channel(channels)
  3652. # bin spikes and subset on intervals
  3653. spikes_nonlocal_binned, _ = sp.calc_binned_spikes(spikes.spikes, \
  3654. intervals=np.append(lfp_nonlocal.timestamps, lfp_nonlocal.timestamps[-1]+1/lfp.rate))
  3655. spikes_nonlocal_fr = np.mean(spikes_nonlocal_binned, axis=1)
  3656. spikes_local_binned, _ = sp.calc_binned_spikes(spikes.spikes, \
  3657. intervals=np.append(lfp_local.timestamps, lfp_local.timestamps[-1]+1/lfp.rate))
  3658. spikes_local_fr = np.mean(spikes_local_binned, axis=1)
  3660. nl_low_theta_coh_units, nl_high_theta_coh_units, l_low_theta_coh_units, l_high_theta_coh_units = \
  3661. [np.full(len(spikes.spikes), np.NaN) for _ in range(4)]
  3663. # calculate rate-adjusted coherence per cell
  3664. for u in range(len(spikes.spikes)):
  3665. _, nl_low_theta_coh_units[u] = lf.multitaper_filtered_spike_coherence(spikes_nonlocal_binned[u][:, np.newaxis], lfp_nonlocal.lfp, \
  3666. spikes_nonlocal_fr[u], spikes_local_fr[u], \
  3667. low=lf.LOW_THETA[0], high=lf.LOW_THETA[1], window=1)
  3668. _, nl_high_theta_coh_units[u] = lf.multitaper_filtered_spike_coherence(spikes_nonlocal_binned[u][:, np.newaxis], lfp_nonlocal.lfp, \
  3669. spikes_nonlocal_fr[u], spikes_local_fr[u], \
  3670. low=lf.HIGH_THETA[0], high=lf.HIGH_THETA[1], window=1)
  3671. _, l_low_theta_coh_units[u] = lf.multitaper_filtered_spike_coherence(spikes_local_binned[u][:, np.newaxis], lfp_local.lfp, \
  3672. spikes_local_fr[u], spikes_nonlocal_fr[u], \
  3673. low=lf.LOW_THETA[0], high=lf.LOW_THETA[1], window=1)
  3674. _, l_high_theta_coh_units[u] = lf.multitaper_filtered_spike_coherence(spikes_local_binned[u][:, np.newaxis], lfp_local.lfp, \
  3675. spikes_local_fr[u], spikes_nonlocal_fr[u], \
  3676. low=lf.HIGH_THETA[0], high=lf.HIGH_THETA[1], window=1)
  3678. # average across the session
  3679. nl_low_theta_coh[animal_idx][sess_idx] = np.nanmean(nl_low_theta_coh_units)
  3680. nl_high_theta_coh[animal_idx][sess_idx] = np.nanmean(nl_high_theta_coh_units)
  3681. l_low_theta_coh[animal_idx][sess_idx] = np.nanmean(l_low_theta_coh_units)
  3682. l_high_theta_coh[animal_idx][sess_idx] = np.nanmean(l_high_theta_coh_units)
  3683. delta_low_theta_coh[animal_idx][sess_idx] = np.nanmean(nl_low_theta_coh_units-l_low_theta_coh_units)
  3684. delta_high_theta_coh[animal_idx][sess_idx] = np.nanmean(nl_high_theta_coh_units-l_high_theta_coh_units)
  3686. print(f"{row['File']} {time.time()-start_time}")
  3688. # %%
  3689. print(f"Type II theta coherence changes by {np.mean(np.nanmean(delta_low_theta_coh, axis=1)):.02f} +/- {stats.sem(np.nanmean(delta_low_theta_coh, axis=1)):.04f}")
  3690. print(f"Coherence: {np.mean(np.nanmean(l_low_theta_coh, axis=1)):.02f} +/- {stats.sem(np.nanmean(l_low_theta_coh, axis=1)):.04f}\n" + \
  3691. f"{np.mean(np.nanmean(nl_low_theta_coh, axis=1)):.02f} +/- {stats.sem(np.nanmean(nl_low_theta_coh, axis=1)):.04f}")
  3693. # %%
  3694. # analyze & plot
  3695. data_by_session = pd.DataFrame({'Animal': np.repeat(np.arange(n_animals),n_sessions),
  3696. 'Session': np.tile(np.arange(n_sessions),n_animals),
  3697. 'Feature_l': nl_low_theta_coh.flatten(),
  3698. 'Feature_nl': l_low_theta_coh.flatten()})
  3699. plot_compare_paired_metrics('Feature_l', 'Feature_nl', data_by_session, \
  3700. figure_path, 'low_theta_spike_field_coh_l_vs_nl', 'Coherence', [0,1])
  3701. data_by_session = pd.DataFrame({'Animal': np.repeat(np.arange(n_animals),n_sessions),
  3702. 'Session': np.tile(np.arange(n_sessions),n_animals),
  3703. 'Feature_l': nl_high_theta_coh.flatten(),
  3704. 'Feature_nl': l_high_theta_coh.flatten()})
  3705. plot_compare_paired_metrics('Feature_l', 'Feature_nl', data_by_session, \
  3706. figure_path, 'high_theta_spike_field_coh_l_vs_nl', 'Coherence', [0,1])
  3708. # %% [markdown]
  3709. # # Figure 4 &

Figures_published_edition.ipynb at commit dfbbde1, under GPL-3.0 · at the source

Overview

Authors: Emily A Aery Jones1,2, Isabel I C Low1,3, Frances S Cho1, Lisa M Giocomo1
ORCID iDs: Emily A Aery Jones, Lisa M Giocomo
  1. Department of Neurobiology, Stanford University School of Medicine, Stanford, CA USA
  2. Present Address: Department of Neurobiology, University of Maryland School of Medicine, Baltimore, MD USA
  3. Present Address: Zuckerman Mind Brain Behavior Institute, Columbia University, New York, NY USA
Institutions: University of Maryland, Baltimore (United States); Stanford Medicine (United States); Stanford University (United States); Mortimer B. Zuckerman Mind Brain Behavior Institute (United States); Columbia University (United States)
Journal: Nature neuroscience, volume 29, issue 5, pages 1181-1190
Dates: received 10 August 2024; accepted 28 January 2026; published online 1 April 2026; in print 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41593-026-02232-0 · PMID 41922514 · PMCID PMC13107481 · OpenAlex W7147527953
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: mouse (organism), systems (subfield)
Methods: Spectral & time-frequency, Connectivity, Statistics, Machine learning, Preprocessing, Single-unit activity, calcium imaging
Keywords: Learning and memory, Neural circuits
MeSH: CA1 Region, Hippocampal*, Entorhinal Cortex*, Neurons*, Space Perception*, Action Potentials, Animals, Male, Mice, Mice, Inbred C57BL, Reward (* major topic)
Topic: Neural dynamics and brain function (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: NIMH NIH HHS (R01 MH106475); NINDS NIH HHS (U19 NS118284, F32 NS138225, K99 NS134734); Howard Hughes Medical Institute
Citations: cited by 1 paper (Europe PMC); 54 references in the paper
Research resources: RRID:IMSR_JAX:000664

Abstract

Neurons can collectively represent the current sensory experience during exploration or remote experiences during immobility. Remote representations can reflect learned associations and support learning. Neurons in medial entorhinal cortex (MEC) represent the animal’s current location during movement, but little is known about MEC representations during immobility. We recorded hundreds of neurons simultaneously in MEC and CA1 as mice learned to associate pairs of rewarded locations. During immobility, the MEC neural population frequently represented positions far from the animal’s location (‘nonlocal coding’). Cells with spatial firing fields at remote locations drove nonlocal coding, even as cells representing the current position remained active. While MEC nonlocal coding has been reported during sharp-wave ripples in CA1, we observed nonlocal coding more often outside of ripples and saw less CA1–MEC coordination during nonlocal coding. Further, nonlocal coding preferentially represented remote task-relevant locations at appropriate times. Together, this work suggests that MEC nonlocal coding could strengthen associations between locations independently from CA1.

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

Repositories

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

emilyasterjones/X_maze

License: MIT
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: b5692c46163cf889a074870ca44af2317f76132e, 24 October 2025
Languages: Python (2)
Size: 46 files, 2 scripts
Software Heritage: not archived
Found in: the text, “Electrophysiology”
Holds: README, license file
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (2 files), OpenCV (1 file), Pillow (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
4 files
At the source: github.com/emilyasterjones/X_maze

petersaj/AP_histology

License: GPL-3.0
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: 67d71af75657dbd14bd97604514337f4822b2a98, 23 September 2026
Languages: MATLAB (31)
Size: 35 files, 31 scripts
Software Heritage: not archived
Found in: the text, “Histology”
Holds: README, license file
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Image Processing Toolbox (9 files), Statistics and Machine Learning Toolbox (3 files)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
33 files
At the source: github.com/petersaj/AP_histology

emilyasterjones/ecephys_spike_sorting

License: BSD-2-Clause
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: b6ef3f0de4f16bf1bba2852edab4988d6345c521, 4 June 2024
Languages: Python (102), C++ (14), MATLAB (9), JavaScript (5), Shell (3), Jupyter (2), C/C++ (2)
Size: 237 files, 137 scripts
Software Heritage: not archived
Found in: the text, “Neural data preprocessing”
Holds: README, license file, environment (Pipfile, Pipfile.lock, setup.cfg, setup.py), tests, continuous integration, documentation, 2 notebooks
Not found: CITATION.cff
Tools: NumPy (51 files), pandas (17 files), Matplotlib (12 files), SciPy (11 files), Kilosort (4 files), scikit-learn (4 files), Phy (3 files), xarray (3 files), h5py (1 file), Pillow (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
139 files
At the source: github.com/emilyasterjones/ecephys_spike_sorting/tree/master

jenniferColonell/ecephys_spike_sorting

License: BSD-2-Clause
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: d61fba8e4387780ea4277c78eb891aacf7a4ee6b, 22 April 2026
Languages: Python (94), C++ (14), JavaScript (5), MATLAB (4), Shell (2), C/C++ (2)
Size: 220 files, 121 scripts
Software Heritage: not archived
Found in: the text, “Neural data preprocessing”
Holds: README, license file, environment (Pipfile, Pipfile.lock, setup.cfg, setup.py), tests, continuous integration, documentation
Not found: CITATION.cff
Tools: NumPy (46 files), pandas (13 files), SciPy (12 files), Matplotlib (9 files), Kilosort (5 files), scikit-learn (4 files), xarray (3 files), h5py (1 file), Phy (1 file), Pillow (1 file), PyTorch (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
123 files
At the source: github.com/jenniferColonell/ecephys_spike_sorting

emilyasterjones/bombcell

License: GPL-3.0
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: b7127c41f35ec6e57973b3aa9b61ec79159e86f2, 11 March 2024
Languages: MATLAB (137), Python (8), C (3), R (1)
Size: 206 files, 149 scripts
Software Heritage: archived
Found in: the text, “Neural data preprocessing”
Holds: README, license file
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Statistics and Machine Learning Toolbox (23 files), Image Processing Toolbox (5 files), Signal Processing Toolbox (5 files), Phy (5 files), NumPy (4 files), SciPy (2 files), Optimization Toolbox (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
151 files
At the source: github.com/emilyasterjones/bombcell

LorenFrankLab/track_linearization

License: MIT
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: 424087881215e3e30bfbf00926571d9e5bef1239, 16 October 2025
Languages: Python (13), Jupyter (3)
Size: 30 files, 16 scripts
Software Heritage: not archived
Found in: the text, “Position tracking”
Holds: README, license file, environment (environment.yml, postBuild, pyproject.toml), tests, continuous integration, 3 notebooks
Not found: CITATION.cff, documentation
Tools: NumPy (14 files), NetworkX (9 files), Matplotlib (7 files), pandas (5 files), SciPy (4 files), Numba (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
18 files
At the source: github.com/LorenFrankLab/track_linearization

Eden-Kramer-Lab/ripple_detection

License: MIT
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: 68e1784b183f161a568d3154325cf770730bfb83, 26 September 2026
Languages: Python (38), Jupyter (4)
Size: 160 files, 42 scripts
Software Heritage: not archived
Found in: the text, “LFP analysis”
Holds: README, license file, CITATION.cff, environment (environment.yml, pyproject.toml, uv.lock), tests, continuous integration, documentation, 4 notebooks
Tools: NumPy (35 files), pandas (24 files), SciPy (12 files), Matplotlib (4 files)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
44 files
At the source: github.com/Eden-Kramer-Lab/ripple_detection

emilyasterjones/AeryJones_2025

License: GPL-3.0
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: dfbbde116929d420523421d7bfe5d441eb91568a, 11 June 2025
Languages: Python (18), Jupyter (6)
Size: 141 files, 24 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, environment (setup.py), 6 notebooks
Not found: CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (19 files), pandas (16 files), SciPy (12 files), Matplotlib (10 files), seaborn (7 files), Neurodata Without Borders (PyNWB, MatNWB) (4 files), SpikeInterface (4 files), CircStat (2 files), DeepLabCut (2 files), OpenCV (2 files), Pingouin (2 files), NetworkX (1 file), Phy (1 file), Pillow (1 file), scikit-image (1 file), scikit-learn (1 file), statsmodels (1 file), xarray (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
26 files
At the source: github.com/emilyasterjones/AeryJones_2025

Code availability

All original code can be found at https://github.com/emilyasterjones/AeryJones_2025.

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

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:

  • 8 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 522 scripts, each with its path and the digest of its content;
  • 19 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 0.1–300-Hz filtered LFP, isolated unit spike times, electrode site locations, trial data, mouse position and head direction, subject metadata and session metadata are available at https://dandiarchive.org/dandiset/001701/0.260120.0303. Source data are provided with this paper.

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 2, 28 September 2026

  • Publisher: n/a → Nature Portfolio

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 4 authors, 2 keywords, 10 MeSH terms, 3 funders, 52 references, 1 RRID.

Cite

This paper

Aery Jones, E. A., Low, I. I. C., Cho, F. S., & Giocomo, L. M. (2026). Entorhinal cortex represents task-relevant remote locations independently of CA1. Nature neuroscience, 29(5), 1181-1190. https://doi.org/10.1038/s41593-026-02232-0

BibTeX

@article{aeryjones2026entorhinal,
author = {Aery Jones, Emily A and Low, Isabel I C and Cho, Frances S and Giocomo, Lisa M},
title = {{Entorhinal cortex represents task-relevant remote locations independently of CA1}},
journal = {Nature neuroscience},
year = {2026},
month = apr,
volume = {29},
number = {5},
pages = {1181--1190},
publisher = {Nature Portfolio},
issn = {1097-6256},
doi = {10.1038/s41593-026-02232-0},
url = {https://doi.org/10.1038/s41593-026-02232-0},
pmid = {41922514},
pmcid = {PMC13107481}
}

RIS

TY - JOUR
AU - Aery Jones, Emily A
AU - Low, Isabel I C
AU - Cho, Frances S
AU - Giocomo, Lisa M
TI - Entorhinal cortex represents task-relevant remote locations independently of CA1
T2 - Nature neuroscience
J2 - Nat Neurosci
PY - 2026
DA - 2026/04/01
VL - 29
IS - 5
SP - 1181
EP - 1190
SN - 1097-6256
PB - Nature Portfolio
DO - 10.1038/s41593-026-02232-0
UR - https://doi.org/10.1038/s41593-026-02232-0
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41593-026-02232-0",
"type": "article-journal",
"title": "Entorhinal cortex represents task-relevant remote locations independently of CA1",
"container-title": "Nature neuroscience",
"author": [
{
"family": "Aery Jones",
"given": "Emily A"
},
{
"family": "Low",
"given": "Isabel I C"
},
{
"family": "Cho",
"given": "Frances S"
},
{
"family": "Giocomo",
"given": "Lisa M"
}
],
"container-title-short": "Nat Neurosci",
"volume": "29",
"issue": "5",
"page": "1181-1190",
"DOI": "10.1038/s41593-026-02232-0",
"PMID": "41922514",
"PMCID": "PMC13107481",
"ISSN": "1097-6256",
"publisher": "Nature Portfolio",
"URL": "https://doi.org/10.1038/s41593-026-02232-0",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
1
]
]
}
}

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.1016/j.patter.2026.101590 [code]
Density-based longitudinal neuron tracking in high-density electrophysiological recordings.
Journal: Patterns (New York, N.Y.)
In common: Phy, Kilosort, SpikeInterface, 14 other tools, 3 references
[2] doi:10.7554/elife.110588 [code]
Opening the black box toward a modular approach to spike sorting.
Journal: eLife
In common: Kilosort, SpikeInterface, Neurodata Without Borders (PyNWB, MatNWB), 12 other tools, mouse, 4 references
[3] doi:10.1038/s41593-026-02362-5 [code]
Replay of procedural memory is independent of the hippocampus.
Journal: Nature neuroscience
In common: Pingouin, h5py, Pillow, 9 other tools, mouse, 9 references
[4] doi:10.1038/s41467-026-75347-4 [code]
Sleep reveals dynamics integrating and segregating movement and stimulus representations in V1.
Journal: Nature communications
In common: Phy, Neurodata Without Borders (PyNWB, MatNWB), CircStat, 13 other tools, systems, mouse, 1 reference
[5] doi:10.1016/j.celrep.2026.117646 [code]
Medial entorhinal-hippocampal desynchronization parallels the emergence of memory impairment in a mouse model of Alzheimer's disease pathology.
Journal: Cell reports
In common: Neurodata Without Borders (PyNWB, MatNWB), CircStat, Optimization Toolbox, 9 other tools, systems, mouse, 5 references
[6] doi: [code]
Naturalistic behavior and self-generated neural activity predictive of self-correction
Journal: bioRxiv : the preprint server for biology
In common: SpikeInterface, Neurodata Without Borders (PyNWB, MatNWB), xarray, 11 other tools, 2 references
[7] doi:10.1038/s41467-026-76581-6 [code]
Thalamocortical bursts encode reward contingencies and drive associative learning.
Journal: Nature communications
In common: Kilosort, xarray, CircStat, 11 other tools, systems, mouse, 2 references
[8] doi:10.7554/elife.109717 [code]
Retrosplenial cortex enables context-dependent goal-directed sensorimotor transformation.
Journal: eLife
In common: Kilosort, Neurodata Without Borders (PyNWB, MatNWB), Numba, 12 other tools, systems, mouse, 1 reference
[9] 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, Numba, NetworkX, 12 other tools, systems, mouse, 2 references
[10] doi:10.1002/hipo.70131 [code]
Decoding Medial Entorhinal Cortical Dynamics Produces Planning-Like Alternations in Hippocampal theta Sequences.
Journal: Hippocampus
In common: xarray, Numba, OpenCV, 6 other tools, systems, 6 references

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.