OSCR

Active dissociation of intracortical spiking and high gamma activity.

Code ↔ Paper

6 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 6 matches
  1. [1] § Methods › ONF performance metrics ↔ behavior_analysis.ipynb, lines 194–260 · score 0.62 · entry angle, HG target trials, spike target trials, monkey CE, metrics, cursor
  2. [2] § Methods › Factor analysis ↔ subspace_analysis.ipynb, lines 139–170 · score 0.62 · linear regression model, concatenated spike, scored, CE HGA, spike rates, weight
  3. [3] § Methods › Neural signal recording ↔ neural_signal_analysis.ipynb, lines 52–114 · score 0.60 · low spike rate, shunted electrodes, neural, channel, signals, binned
  4. [4] § Distributed spikes do not leak into HGA ↔ subspace_analysis.ipynb, lines 197–228 · score 0.60 · weighted sum, randomly selected electrodes, DWSS, spike rate, correlated, CE
  5. [5] § Distributed spikes do not leak into HGA ↔ subspace_analysis.ipynb, lines 197–228 · score 0.56 · random DWSS, weighted sum, randomly selected electrodes, Correlation, CE, HGA
  6. [6] § Methods › cPCA ↔ contrastive/__init__.py, lines 12–20 · score 0.55 · eigenvalue decomposition, cPCA, foreground, variance, background

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 · 515 lines · 20 KB · CC-BY-4.0 · 3 matches

  1. # %% [markdown]
  2. # # This Jupyter notebook file produces the subspace analysis results in the paper
  4. # %%
  5. import numpy as np
  6. import scipy
  7. import mat73
  8. import sklearn
  9. from utility_functions import *
  10. import pandas as pd
  11. import matplotlib.pyplot as plt
  12. import matplotlib.patches as Patches
  13. from itertools import chain
  14. from sklearn.preprocessing import StandardScaler
  15. from sklearn.decomposition import FactorAnalysis,PCA
  16. from sklearn.linear_model import LinearRegression
  17. import matplotlib.ticker as mticker
  18. from contrastive import CPCA
  20. plt.rcParams["axes.linewidth"] = 2
  21. plt.rcParams["axes.edgecolor"] = 'k'
  22. plt.rcParams['grid.alpha'] = 0
  24. # %% [markdown]
  25. # # Load the data and meta data
  26. # Here we need to load the entire trials. So `data_type` should be "trials"
  28. # %%
  29. data_type='trials' # this could be 'trials', 'last4s', 'first4s'.
  30. monkey='C' # monkey should be one of 'C', 'J' or 'M'
  31. CE=63 # channel name. NOTE: this is not the channel number or the channel index
  33. # load the grid layout of the monkey
  34. grid=get_grid(monkey)
  35. # load the spike channel number
  36. spike_channel=get_channel(monkey,CE)
  37. # load the location of the CE on the grid
  38. loc_x,loc_y=get_CE_loc(grid,spike_channel)
  39. # load the shunted electrodes for this monkey. This is the python index
  40. shunted_electrodes=load_shunted_electrodes(monkey,spike_channel)
  41. # finally, load the raw data
  42. raw_data=load_ONF_data(data_type, monkey, CE)
  43. # We only need the "learned" files.
  44. raw_data=raw_data[raw_data['file_types']=='learned']
  46. # %% [markdown]
  47. # Optional: Visualize the data table. This may take some time.
  49. # %%
  50. raw_data
  52. # %% [markdown]
  53. # # Pre-process the spike rate and HGA data
  54. #
  55. # 1. include only the successful trials.
  56. # 2. remove shunted channels
  57. # 3. remove NAN and empty channels
  58. # 4. remove channels with low firing rate
  60. # %%
  61. # Here we concatenate all the successful trials together for every file.
  62. spike_rate_all = list(chain.from_iterable([[d[0] for d in data] for data in raw_data['spike_rate'].to_list()]))
  63. HGA_all = list(chain.from_iterable([[d[0] for d in data] for data in raw_data['HGA'].to_list()]))
  65. # We concatenate the trial type
  66. trial_types_all=np.concatenate(raw_data['trial_types'].to_list())
  68. # We concatenate the target type
  69. target_types_all=np.concatenate(raw_data['target_types'].to_list())
  71. # We take out only the successful trials
  72. spike_rate_successful=[a for a,t in zip(spike_rate_all,trial_types_all) if t==1 ]
  73. HGA_successful=[a for a,t in zip(HGA_all,trial_types_all) if t==1]
  75. # concatenate the signals
  76. spike_rate_flat=np.concatenate(spike_rate_successful,axis=-1)
  77. HGA_flat=np.concatenate(HGA_successful,axis=-1)
  78. n_time=spike_rate_flat.shape[0]
  80. # save the indices of the leftover electrodes
  81. leftover_electrode_list=np.arange(0,spike_rate_flat.shape[0])
  83. print(f"Before preprocessing, signals have the shape (chan x (time x n_trials)): {spike_rate_flat.shape}")
  84. # remove shunted electrodes
  85. spike_rate_flat=np.delete(spike_rate_flat,shunted_electrodes,axis=0)
  86. HGA_flat=np.delete(HGA_flat,shunted_electrodes,axis=0)
  87. leftover_electrode_list=np.delete(leftover_electrode_list,shunted_electrodes)
  88. print(f"Removed the shunted electrodes, signals have the shape (chan x (time x n_trials)): {spike_rate_flat.shape}")
  90. # check NAN electrodes in the data
  91. has_nan_spike_rate=np.sum(np.isnan(spike_rate_flat),axis=1)
  92. has_nan_HGA=np.sum(np.isnan(HGA_flat),axis=1)
  93. has_nan=np.where(has_nan_spike_rate+has_nan_HGA)
  95. # remove the electrodes with NAN values
  96. spike_rate_flat=np.delete(spike_rate_flat,has_nan,axis=0)
  97. HGA_flat=np.delete(HGA_flat,has_nan,axis=0)
  98. leftover_electrode_list=np.delete(leftover_electrode_list,has_nan,axis=0)
  99. print(f"Removed the electrodes with NAN values, signals have the shape (chan x (time x n_trials)): {spike_rate_flat.shape}")
  101. # check the electrodes with constant low spike rate
  102. low_spike_rate=5 # spikes/s
  103. has_low_spike_rate=np.sum(spike_rate_flat,axis=1)<n_time*low_spike_rate/20 # 20 Hz
  104. has_low_spike_rate=np.where(has_low_spike_rate==spike_rate_flat.shape[-1])
  106. # remove the electrodes with low spike rate
  107. spike_rate_flat=np.delete(spike_rate_flat,has_low_spike_rate,axis=0)
  108. HGA_flat=np.delete(HGA_flat,has_low_spike_rate,axis=0)
  109. leftover_electrode_list=np.delete(leftover_electrode_list,has_low_spike_rate)
  110. print(f"Removed the electrodes with low spike rate, signals have the shape (chan x (time x n_trials)): {spike_rate_flat.shape}")
  112. # get the number of clean electrodes
  113. n_clean_electrodes=len(leftover_electrode_list)
  114. CE_index=np.where(leftover_electrode_list==spike_channel-1)[0]
  116. # %% [markdown]
  117. # ### Optional: check the correlation between CE neural siginals and velocities
  119. # %%
  120. # get the velocity
  121. velocity=list(chain.from_iterable(raw_data['velocity'].to_list()))
  122. velocity=np.concatenate([v[0] for v, t in zip(velocity,trial_types_all) if t==1.0],axis=0)
  124. # plot the correlation of velocity x and spike rate, velocity y and HGA
  125. fig,ax=plt.subplots(1,2, figsize=(5,3))
  126. ax[0].scatter(np.squeeze(velocity[:,0]),np.squeeze(spike_rate_flat[CE_index,:]),color='b',alpha=0.2,s=2,marker='.')
  127. ax[0].set_xlabel('velocity x')
  128. ax[0].set_ylabel('spike rate (Hz)')
  130. ax[1].scatter(np.squeeze(velocity[:,1]),np.squeeze(HGA_flat[CE_index,:]),color='b',alpha=0.2,s=2,marker='.')
  131. ax[1].set_xlabel('velocity y')
  132. ax[1].set_ylabel('HGA amplitude')
  134. # title
  135. plt.suptitle("Correlation between the neural signals and velocity")
  136. plt.tight_layout()
  137. plt.show()
  139. # %% [markdown]
  140. # # Comparing the first factor to CE HGA and spike rate
  142. # %%
  143. # factor analysis on the mean centered, concatenated spike rate
  144. n_factors=10
  146. # zscore the spike rate data
  147. spike_rate_zscore=scipy.stats.zscore(spike_rate_flat,axis=None)
  149. # run factor analysis
  150. fa=FactorAnalysis(n_components=n_factors,rotation='varimax')
  151. factors=fa.fit_transform(spike_rate_zscore.T)
  152. weights=fa.components_
  154. # calculate the R between the factors and CE neural signals
  156. R_factors_spike_rate=[]
  157. R_factors_HGA=[]
  158. for i in range(n_factors):
  159. # fit a linear regression model on factors and spike rate
  160. lr_factors_spike_rate=LinearRegression().fit(factors[:,:i+1],spike_rate_zscore[CE_index,:].T)
  161. # calculate the R2
  162. r2_factors_spike_rate=lr_factors_spike_rate.score(factors[:,:i+1],spike_rate_zscore[CE_index,:].T)
  163. R_factors_spike_rate.append(np.sqrt(r2_factors_spike_rate))
  166. # fit a linear regression model on factors and HGA
  167. lr_factors_HGA=LinearRegression().fit(factors[:,:i+1],HGA_flat[CE_index,:].T)
  168. # calculate the R2
  169. r2_factors_HGA=lr_factors_HGA.score(factors[:,:i+1],HGA_flat[CE_index,:].T)
  170. R_factors_HGA.append(np.sqrt(r2_factors_HGA))
  172. # %% [markdown]
  173. # plot R and the factor weights
  175. # %%
  176. fig,ax = plt.subplots(1,2, figsize=(7,3))
  177. im=plot_on_grid(ax[0],
  178. leftover_electrode_list,grid,
  179. np.abs(weights[0]),(loc_x,loc_y),
  180. cmap=mpb.cm.viridis
  181. )
  182. im.set_clim(0,1)
  183. cb=plt.colorbar(im,shrink=0.8)
  184. cb.set_label('First factor weights',fontsize=12)
  186. l1=ax[1].plot(np.arange(1,n_factors+1),R_factors_spike_rate,lw=3, label='spike rate',color='magenta')
  187. l2=ax[1].plot(np.arange(1,n_factors+1),R_factors_HGA,lw=3,label='HGA',color='darkcyan')
  188. ax[1].set_xticks(np.arange(1,n_factors+1,3))
  189. ax[1].set_yticks(np.arange(0,1,0.2))
  190. ax[1].set_xlabel('number of factors')
  191. ax[1].set_ylabel('R')
  192. ax[1].legend(title='Factors fitted to')
  194. plt.tight_layout()
  195. plt.show()
  197. # %% [markdown]
  198. # ### here we calculate the correlation between DWSS and HGA
  200. # %%
  201. # calculate the weighted grid for CE
  203. weighted_grid=calculate_weighted_grid((loc_x,loc_y),power=2)
  204. weighted_grid=get_the_weights(weighted_grid,
  205. leftover_electrode_list,
  206. grid.flatten())
  208. # calculate the weighted sum of spike using the weights_grid
  209. corr_HGA_DWSS=np.corrcoef(
  210. np.vstack([HGA_flat[CE_index],np.average(spike_rate_zscore,axis=0,weights=weighted_grid)]))[0,1]
  212. # calculate the weighted grid for 20 random electrodes
  213. n_random_pairs=20
  214. random_pairs=draw_random_pairs(n_random_pairs,except_loc=(loc_x,loc_y))
  215. random_weighted_grids=[get_the_weights(calculate_weighted_grid(r_pairs,power=2,emit_start_location=False),leftover_electrode_list,grid.flatten()) for r_pairs in random_pairs]
  217. # calculate the weighted sum of spike using random weighted grids.
  218. corr_HGA_random_DWSS=[np.corrcoef(np.vstack([HGA_flat[CE_index],np.average(spike_rate_zscore,axis=0,weights=r_grid)]))[0,1] for r_grid in random_weighted_grids]
  220. # calculate the DWSS and HGA on non CE electrodes
  221. # get the randomly selected electrodes
  222. random_electrodes=np.random.choice(np.delete(leftover_electrode_list,np.where(leftover_electrode_list==spike_channel-1)),n_random_pairs,replace=False)+1
  223. random_electrodes_loc=[np.where(grid==r) for r in random_electrodes]
  225. # calculate the DWSS
  226. random_weighted_grids=[get_the_weights(calculate_weighted_grid(r,power=2),leftover_electrode_list,grid.flatten()) for r in random_electrodes_loc]
  227. # calculate the correlation between DWSS and HGA on those random electrodes
  228. corr_random_HGA_DWSS=[np.corrcoef(np.squeeze(HGA_flat[np.where(leftover_electrode_list==random_electrodes[n]-1)[0]]), np.average(spike_rate_zscore,axis=0,weights=random_weighted_grids[n]))[0,1] for n in range(n_random_pairs)]
  230. # %% [markdown]
  231. # plotting
  233. # %%
  234. fig,ax=plt.subplots(1,1,figsize=(3,3))
  235. im=plot_on_grid(ax,leftover_electrode_list,grid,np.sqrt(weighted_grid),(loc_x,loc_y),cmap=mpb.cm.Blues)
  236. im.set_clim(0,0.003)
  237. ax.set_title("distance weights",fontsize=10)
  238. ax.spines[['left','bottom','top','right']].set_linewidth(3)
  239. cb=fig.colorbar(im,shrink=0.8,ticks=[0,0.002],format=mticker.FixedFormatter(['0', '3']))
  240. cb.set_label("weights",size=10)
  241. cb.ax.set_yticklabels(['0', '3']) # vertically oriented colorbar
  242. cb.ax.tick_params(labelsize=10)
  245. fig,ax=plt.subplots(1,1,figsize=(3,2))
  246. ax.bar([0,1,2],
  247. np.abs(np.squeeze([
  248. corr_HGA_DWSS,
  249. np.mean(corr_HGA_random_DWSS),
  250. np.mean(corr_random_HGA_DWSS)]))
  251. ,color=['crimson','grey',"lightblue"],yerr=(0,np.std(corr_HGA_random_DWSS),np.std(corr_random_HGA_DWSS)),error_kw=dict(lw=3))
  252. # plot the raw data
  253. x = np.random.normal(1, 0.05, size=len(corr_HGA_random_DWSS)) # Jitter
  254. ax.scatter(x, np.array(corr_HGA_random_DWSS), color='k', alpha=0.5, edgecolor='black', linewidth=0)
  255. x = np.random.normal(2, 0.05, size=len(corr_random_HGA_DWSS)) # Jitter
  256. ax.scatter(x, np.array(corr_random_HGA_DWSS), color='k', alpha=0.5, edgecolor='black', linewidth=0)
  258. ax.set_xticks([0,1,2],["DWSS","random","non-CE"],rotation=0,fontsize=10)
  260. ax.set_ylabel("R",fontsize=10)
  261. ax.set_ylim(-0.2,0.8)
  262. ax.tick_params(axis='y', labelsize=10)
  263. ax.tick_params(axis='x', labelsize=10)
  264. ax.spines[['top','right']].set_visible(False)
  265. ax.spines[['left','bottom']].set_linewidth(3)
  267. plt.show()
  270. # %% [markdown]
  271. # # PCA and cPCA analysis on the correlation and subspace angles
  272. # 1. define SP+ and HG+ NAS
  273. # 2. perform PCA on the data separately
  274. # 3. perform cPCA on the data separately
  275. # 4. calculate the angles between them
  277. # %%
  278. # calculate the CE spike rate and HGA
  279. CE_spike_rate=StandardScaler().fit_transform(spike_rate_zscore[CE_index].T)
  280. CE_HGA=StandardScaler().fit_transform(HGA_flat[CE_index].T)
  282. # calculated the sp+ and hg+ indices
  283. NAS=np.squeeze(np.arctan2(CE_HGA, CE_spike_rate)*180/np.pi)
  284. sp_ind=(NAS>=-15)&(NAS<=15)
  285. hg_ind=(NAS>=75)&(NAS<=105)
  287. # %% [markdown]
  288. # Quick check of the dimensionality of the two data
  290. # %%
  291. n_PCs=10
  292. pca=PCA(n_components=n_PCs)
  293. plt.figure()
  294. plt.plot(np.cumsum(pca.fit(spike_rate_zscore[:,sp_ind].T).explained_variance_ratio_),lw=3,color='magenta')
  295. plt.plot(np.cumsum(pca.fit(spike_rate_zscore[:,hg_ind].T).explained_variance_ratio_),lw=3,color='darkblue')
  296. plt.ylabel("cumulative VAF")
  297. plt.xlabel("number of PCs")
  298. plt.legend(['SP+','HG+'])
  300. # %%
  301. # first we need to fine the best alpha for cPCA
  302. cpca=CPCA(n_components=n_PCs,standardize=False)
  304. # specify the range the step for the alphas to try
  305. n_alphas=15
  306. max_log_alpha=2
  308. # prepare the data
  309. sp_data=scipy.stats.zscore(spike_rate_zscore[:,sp_ind],axis=None)
  310. hg_data=scipy.stats.zscore(spike_rate_zscore[:,hg_ind],axis=None)
  312. # Use the plots to determin the best alpha
  313. projected_fg,best_alphas=cpca.fit_transform(hg_data.T,
  314. sp_data.T,
  315. n_alphas=n_alphas, max_log_alpha=max_log_alpha, alpha_selection='all',return_alphas=True)
  316. projected_bg=cpca.transform(sp_data.T,n_alphas=n_alphas, max_log_alpha=max_log_alpha, alpha_selection='all')
  318. # %% [markdown]
  319. # - Plotting the background and foreground data with different alpha
  320. # - The peak of the ratio of variances in the second plot is where the best alpha value is
  322. # %%
  323. plt.figure(figsize=[3*n_alphas+1,3])
  324. subsample=0.2 # only plot number of subsample x data
  325. all_ratio=[]
  326. for j, pack in enumerate(zip(projected_fg,projected_bg)):
  327. # retrieve and subsample the background and foreground data
  328. fg,bg=pack
  329. subsample_bg=np.random.choice(len(bg),int(len(bg)*subsample),replace=False)
  330. subsample_fg=np.random.choice(len(fg),int(len(fg)*subsample),replace=False)
  332. # calculate the variance
  333. fg_variance=np.var(fg)
  334. bg_variance=np.var(bg)
  336. plt.subplot(1,n_alphas+1,j+1)
  337. plt.scatter(bg[subsample_bg,0],bg[subsample_bg,1], color='k', alpha=0.5, label='background',s=3)
  338. plt.scatter(fg[subsample_fg,0],fg[subsample_fg,1], color='r', alpha=0.5, label='foreground',s=3)
  339. plt.title('Alpha='+str(np.round(best_alphas[j],2)))
  340. plt.xlabel(f"background {bg_variance:0.2f}\nforeground {fg_variance:0.2f}\nratio {fg_variance/bg_variance:0.2f}")
  341. all_ratio.append(fg_variance/bg_variance)
  342. plt.legend()
  344. # plot the variance ratios for each alpha we selected
  345. plt.figure()
  346. plt.plot(best_alphas,all_ratio,lw=3)
  347. plt.xscale("log")
  348. plt.yticks(np.linspace(1,19,6))
  349. plt.xlabel('alpha value')
  350. plt.ylabel('ratio of variances')
  351. plt.ylim(0,5)
  353. print(f"The best alpha value for this dataset it: {best_alphas[np.argmax(all_ratio)]}")
  354. plt.show()
  356. # %% [markdown]
  357. # ## Now we are using the best alphao to produce the cPCs and correlate them with CE spike rate and CE HGA
  359. # %%
  360. alpha_selected=best_alphas[np.argmax(all_ratio)]
  362. ######################## fg: SP+, bg: HG+ ################################
  363. fg_cov,bg_cov=get_cov(sp_data.T,hg_data.T)
  364. v_top_SP=get_cpca_loadings(fg_cov,bg_cov,alpha_selected,n_PCs)
  366. projected_bg_SP=hg_data.T @ v_top_SP
  367. projected_fg_SP=sp_data.T @ v_top_SP
  369. # fit the unique component one by one and see which component(s) contribute the most
  370. cumsum_R2_SP_CPCA_FR=[]
  371. cumsum_R2_SP_CPCA_HG=[]
  372. for i in range(n_PCs):
  373. R2=(LinearRegression(fit_intercept=False)
  374. .fit(spike_rate_zscore.T@v_top_SP[:,:i+1],np.squeeze(spike_rate_zscore[CE_index]))
  375. .score(spike_rate_zscore.T@v_top_SP[:,:i+1],np.squeeze(spike_rate_zscore[CE_index])))
  376. cumsum_R2_SP_CPCA_FR.append(np.maximum(R2,0))
  378. R2=(LinearRegression(fit_intercept=False)
  379. .fit(spike_rate_zscore.T@v_top_SP[:,:i+1],np.squeeze(HGA_flat[CE_index]))
  380. .score(spike_rate_zscore.T@v_top_SP[:,:i+1],np.squeeze(HGA_flat[CE_index])))
  381. cumsum_R2_SP_CPCA_HG.append(np.maximum(R2,0))
  384. ######################### fg: SP+, bg: HG+ ################################
  385. fg_cov,bg_cov=get_cov(hg_data.T,sp_data.T)
  386. v_top_HG=get_cpca_loadings(fg_cov,bg_cov,alpha_selected,n_PCs)
  388. projected_bg_HG=sp_data.T @ v_top_HG
  389. projected_fg_HG=hg_data.T @ v_top_HG
  391. # fit the unique component one by one and see which component(s) contribute the most
  392. cumsum_R2_HG_CPCA_FR=[]
  393. cumsum_R2_HG_CPCA_HG=[]
  394. for i in range(n_PCs):
  395. R2=(LinearRegression(fit_intercept=False)
  396. .fit(spike_rate_zscore.T@v_top_HG[:,:i+1],np.squeeze(spike_rate_zscore[CE_index,:]))
  397. .score(spike_rate_zscore.T@v_top_HG[:,:i+1],np.squeeze(spike_rate_zscore[CE_index,:])))
  398. cumsum_R2_HG_CPCA_FR.append(np.maximum(R2,0))
  400. R2=(LinearRegression(fit_intercept=False)
  401. .fit(spike_rate_zscore.T@v_top_HG[:,:i+1],np.squeeze(HGA_flat[CE_index,:]))
  402. .score(spike_rate_zscore.T@v_top_HG[:,:i+1],np.squeeze(HGA_flat[CE_index,:])))
  403. cumsum_R2_HG_CPCA_HG.append(np.maximum(R2,0))
  405. # %% [markdown]
  406. # Now we perform PCA on SP+ data and HG+ data separately
  408. # %%
  409. # same analysis using PCA (this equals to fg with alpha=0)
  410. pca_projection_HG=PCA(n_PCs).fit_transform(hg_data.T)
  411. pca_projection_SP=PCA(n_PCs).fit_transform(sp_data.T)
  412. pca_loading_HG=PCA(n_PCs).fit(hg_data.T).components_.T
  413. pca_loading_SP=PCA(n_PCs).fit(sp_data.T).components_.T
  414. shuffle_pca_loading_HG=np.copy(pca_loading_HG)
  415. shuffle_pca_loading_SP=np.copy(pca_loading_SP)
  416. np.random.shuffle(shuffle_pca_loading_HG)
  417. np.random.shuffle(shuffle_pca_loading_SP)
  420. # %% [markdown]
  421. # Prepare the shuffled dataset as control
  423. # %%
  424. sp_fr_shuffle = np.copy(sp_data)
  425. np.random.shuffle(sp_fr_shuffle)
  426. hg_fr_shuffle = np.copy(hg_data)
  427. np.random.shuffle(hg_fr_shuffle)
  428. fg_cov_shuffle,bg_cov_shuffle=get_cov(sp_fr_shuffle.T,hg_fr_shuffle.T)
  429. v_top_SP_shuffle=get_cpca_loadings(fg_cov_shuffle,bg_cov_shuffle,alpha_selected,n_PCs)
  431. np.random.shuffle(sp_fr_shuffle)
  432. np.random.shuffle(hg_fr_shuffle)
  433. bg_cov_shuffle,fg_cov_shuffle=get_cov(sp_fr_shuffle.T,hg_fr_shuffle.T)
  434. v_top_HG_shuffle=get_cpca_loadings(fg_cov_shuffle,bg_cov_shuffle,alpha_selected,n_PCs)
  436. # %% [markdown]
  437. # Calculate the subspace angles between different PCs and cPCs
  439. # %%
  440. subspace_angles_cPCs=scipy.linalg.subspace_angles(v_top_HG, v_top_SP)
  441. subspace_angles_cPCs_shuffle=scipy.linalg.subspace_angles(v_top_HG_shuffle, v_top_SP_shuffle)
  443. subspace_angles_PCs=scipy.linalg.subspace_angles(pca_loading_HG, pca_loading_SP)
  444. subspace_angles_PCs_shuffle=scipy.linalg.subspace_angles(shuffle_pca_loading_HG, shuffle_pca_loading_SP)
  446. subspace_angles_cPCA_PCA_HG=scipy.linalg.subspace_angles(v_top_HG, pca_loading_HG)
  447. subspace_angles_cPCA_PCA_SP=scipy.linalg.subspace_angles(v_top_SP, pca_loading_SP)
  448. subspace_angles_shuffle_cPCA_PCA_HG=scipy.linalg.subspace_angles(v_top_HG_shuffle, shuffle_pca_loading_HG)
  449. subspace_angles_shuffle_cPCA_PCA_SP=scipy.linalg.subspace_angles(v_top_SP_shuffle, shuffle_pca_loading_SP)
  451. # %% [markdown]
  452. # Here we demostrated the correlation (R) between the spike rate, HGA and SP+/HG+ cPCs
  454. # %%
  455. fig,ax=plt.subplots(1,1,figsize=(4,3))
  456. ax.plot(np.linspace(1,3*n_PCs,n_PCs),np.sqrt(cumsum_R2_HG_CPCA_HG),color='tab:blue',lw=3)
  457. ax.plot(np.linspace(1,3*n_PCs,n_PCs),np.sqrt(cumsum_R2_HG_CPCA_FR),color='grey',lw=3)
  459. ax.set_xticks(np.linspace(2,3*n_PCs+0.5,5),np.int16(np.linspace(1,n_PCs,5)),fontsize=14)
  460. ax.set_yticks(np.linspace(0,1,3),np.linspace(0,1,3),fontsize=14)
  461. ax.legend(["HGA","Spike rate"],frameon=False,loc="upper right",title="HG+ cPCs fitted to",title_fontsize=14,fontsize=10)
  462. ax.set_xlabel("num of HG+ cPCs added",fontsize=14)
  463. ax.set_ylabel("R",fontsize=14)
  464. ax.spines['top'].set_visible(False)
  465. ax.spines['right'].set_visible(False)
  466. plt.tight_layout()
  469. fig,ax=plt.subplots(1,1,figsize=(4,3))
  470. ax.plot(np.linspace(1,3*n_PCs,n_PCs),np.sqrt(cumsum_R2_SP_CPCA_HG),color='grey',lw=2)
  471. ax.plot(np.linspace(1,3*n_PCs,n_PCs),np.sqrt(cumsum_R2_SP_CPCA_FR),color='crimson',lw=2)
  473. ax.set_xticks(np.linspace(2,3*n_PCs+0.5,5),np.int16(np.linspace(1,n_PCs,5)),fontsize=14)
  474. ax.set_yticks(np.linspace(0,1,3),np.linspace(0,1,3),fontsize=15)
  475. ax.legend(["HGA","Spike rate"],frameon=False,loc="best",title="SP+ cPCs fitted to",title_fontsize=14,fontsize=10)
  476. ax.set_xlabel("num of SP+ cPCs added",fontsize=14)
  477. ax.set_ylabel("R",fontsize=14)
  478. ax.spines['top'].set_visible(False)
  479. ax.spines['right'].set_visible(False)
  480. plt.tight_layout()
  482. plt.show()
  484. # %% [markdown]
  485. # Here we demonstrate the cPCA weights and that HG+ cPCs have more activation in other electrodes, compared to SP+ cPCs
  487. # %%
  488. ### another figure: plot the first cPC on the grid
  489. first_cPC_weights_SP=np.abs(v_top_SP[:,0])/np.abs(v_top_SP[:,0]).max()
  490. first_cPC_weights_HG=np.abs(v_top_HG[:,0])/np.abs(v_top_HG[:,0]).max()
  492. fig,ax=plt.subplots(1,2,layout='constrained',figsize=(7,4))
  493. im1=plot_on_grid(ax[0],leftover_electrode_list,grid,first_cPC_weights_SP,(loc_x,loc_y),
  494. cmap=mpb.cm.viridis,show_chan_indx=False,CE_color='r')
  495. im2=plot_on_grid(ax[1],leftover_electrode_list,grid,first_cPC_weights_HG,(loc_x,loc_y),
  496. cmap=mpb.cm.viridis,show_chan_indx=False,CE_color='r')
  497. #clim_min=np.min([v_top_SP[:,0],v_top_HG[:,0]])
  498. clim=np.max([first_cPC_weights_SP,first_cPC_weights_HG])
  499. #clim=np.maximum(clim_max,np.abs(clim_min))
  500. im1.set_clim([0,clim-0.3])
  501. im2.set_clim([0,clim-0.1])
  502. cb=fig.colorbar(im2, ax=ax[:], location='right',shrink=0.6,pad=0.1)
  503. cb.set_label("cPC weights",fontsize=12)
  504. cb.ax.tick_params(labelsize=12)
  505. ax[0].set_title("SP+ 1st cPC",fontsize=12)
  506. ax[1].set_title("HG+ 1st cPC",fontsize=12)
  508. ### some quantification: save the first cPC weights
  509. fig,ax=plt.subplots(1,1,figsize=(4,4))
  510. ax.violinplot([first_cPC_weights_SP,first_cPC_weights_HG])
  511. plt.show()
  513. # %% [markdown]
  514. # ---
  515. # Code developed by: [Tianhao Lei](https://github.com/caraido)

subspace_analysis.ipynb at commit 8ca016b, under CC-BY-4.0 · at the source

Overview

Authors: Tianhao Lei1, Michael R Scheid1, Robert D Flint1, Joshua I Glaser1,2,3, Marc W Slutzky1,4,5,6
ORCID iDs: Tianhao Lei, Michael R Scheid, Marc W Slutzky
  1. Department of Neurology, Northwestern University Feinberg School of Medicine, Chicago, IL USA
  2. Department of Computer Science, Northwestern University, Evanston, IL USA
  3. National Institute for Theory and Mathematics in Biology, Chicago, IL USA
  4. Department of Neuroscience, Northwestern University Feinberg School of Medicine, Chicago, IL USA
  5. Department of Physical Medicine and Rehabilitation, Northwestern University Feinberg School of Medicine, Chicago, IL USA
  6. Department of Biomedical Engineering, Northwestern University, Evanston, IL USA
Institutions: Northwestern University (United States)
Journal: Nature, volume 654, issue 8119, pages 744-750
Dates: received 6 September 2024; accepted 26 February 2026; published online 1 April 2026; in print 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1038/s41586-026-10331-y · PMID 41922776 · PMCID PMC13067879 · OpenAlex W7147612737
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: extracellular electrophysiology (units, LFP) (modality), non-human primate (organism)
Methods: Spectral & time-frequency, Connectivity, Statistics, Preprocessing, Graphs
Keywords: Neuronal physiology, Network models, Biophysical models, Neurophysiology, Extracellular recording
MeSH: Action Potentials*, Gamma Rhythm*, Macaca mulatta*, Neurons*, Animals, Brain-Computer Interfaces, Male, Models, Neurological, Synaptic Potentials (* major topic)
Topic: Neural dynamics and brain function (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: NINDS NIH HHS (RF1 NS125026, T32 NS047987, K08 NS060223, R01 NS112942, R01 NS094748, R00 NS119787, R01 NS099210); NIBIB NIH HHS (T32 EB009406)
Citations: cited by 8 papers (Europe PMC); 65 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repositories

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

abidlabs/contrastive

License: MIT
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: f0d3c330009cfdc76675f95859e70dc96adc3689, 7 October 2025
Languages: Jupyter (13), Python (6)
Size: 49 files, 19 scripts
Software Heritage: 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), Matplotlib (15 files), scikit-learn (12 files), SciPy (7 files), Pillow (3 files), h5py (2 files), TensorFlow (2 files)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
21 files
At the source: github.com/abidlabs/contrastive

caraido/spike-highgamma

License: CC-BY-4.0
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: 8ca016b16e8b0452bec726fab1c0b118a77895dc, 4 May 2026
Languages: Jupyter (4), Python (1)
Size: 9 files, 5 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, environment (environment.yaml), 4 notebooks
Not found: license file, CITATION.cff, tests, continuous integration, documentation
Tools: Matplotlib (5 files), SciPy (5 files), NumPy (4 files), pandas (4 files), scikit-learn (4 files)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
6 files
At the source: github.com/caraido/spike-highgamma

Code availability statement

The paper has a code availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

Read it in the paper: doi.org/10.1038/s41586-026-10331-y.

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:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 24 scripts, each with its path and the digest of its content;
  • 6 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 statement

The paper has a data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

Read it in the paper: doi.org/10.1038/s41586-026-10331-y.

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, 5 authors, 5 keywords, 9 MeSH terms, 2 funders, 64 references.

Cite

This paper

Lei, T., Scheid, M. R., Flint, R. D., Glaser, J. I., & Slutzky, M. W. (2026). Active dissociation of intracortical spiking and high gamma activity. Nature, 654(8119), 744-750. https://doi.org/10.1038/s41586-026-10331-y

BibTeX

@article{lei2026active,
author = {Lei, Tianhao and Scheid, Michael R and Flint, Robert D and Glaser, Joshua I and Slutzky, Marc W},
title = {{Active dissociation of intracortical spiking and high gamma activity}},
journal = {Nature},
year = {2026},
month = apr,
volume = {654},
number = {8119},
pages = {744--750},
publisher = {Nature Portfolio},
issn = {0028-0836},
doi = {10.1038/s41586-026-10331-y},
url = {https://doi.org/10.1038/s41586-026-10331-y},
pmid = {41922776},
pmcid = {PMC13067879}
}

RIS

TY - JOUR
AU - Lei, Tianhao
AU - Scheid, Michael R
AU - Flint, Robert D
AU - Glaser, Joshua I
AU - Slutzky, Marc W
TI - Active dissociation of intracortical spiking and high gamma activity
T2 - Nature
J2 - Nature
PY - 2026
DA - 2026/04/01
VL - 654
IS - 8119
SP - 744
EP - 750
SN - 0028-0836
PB - Nature Portfolio
DO - 10.1038/s41586-026-10331-y
UR - https://doi.org/10.1038/s41586-026-10331-y
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41586-026-10331-y",
"type": "article-journal",
"title": "Active dissociation of intracortical spiking and high gamma activity",
"container-title": "Nature",
"author": [
{
"family": "Lei",
"given": "Tianhao"
},
{
"family": "Scheid",
"given": "Michael R"
},
{
"family": "Flint",
"given": "Robert D"
},
{
"family": "Glaser",
"given": "Joshua I"
},
{
"family": "Slutzky",
"given": "Marc W"
}
],
"container-title-short": "Nature",
"volume": "654",
"issue": "8119",
"page": "744-750",
"DOI": "10.1038/s41586-026-10331-y",
"PMID": "41922776",
"PMCID": "PMC13067879",
"ISSN": "0028-0836",
"publisher": "Nature Portfolio",
"URL": "https://doi.org/10.1038/s41586-026-10331-y",
"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.1073/pnas.2516293123
Distinct laminar origins of sensory-evoked high-gamma and low-frequency ECoG signals revealed by optogenetics.
Journal: Proceedings of the National Academy of Sciences of the United States of America
In common: 13 references
[2] doi:10.1371/journal.pbio.3003873 [code]
Distinct sources of decision-related signals in visual cortex are represented in different local field potential bands.
Journal: PLoS biology
In common: pandas, Matplotlib, NumPy, extracellular electrophysiology (units, LFP), non-human primate, 6 references
[3] doi:10.1038/s41467-026-75455-1 [code]
Shared latent representations of speech production for cross-patient speech decoding.
Journal: Nature communications
In common: TensorFlow, h5py, scikit-learn, 4 other tools, 4 references
[4] doi:10.1007/s12021-026-09807-z [code]
Optimal Size of Electrocorticography Grids for Classification of Hand Movements.
Journal: Neuroinformatics
In common: scikit-learn, pandas, SciPy, 2 other tools, 5 references
[5] doi:10.64898/2026.03.12.710517 [code]
Cortical excitability inversely modulates fMRI connectivity via low-frequency neuronal coupling
Journal: bioRxiv (preprint)
In common: h5py, pandas, SciPy, 2 other tools, 5 references
[6] doi:10.1038/s41467-026-72444-2 [code]
Broadband synergy versus oscillatory redundancy in the visual cortex.
Journal: Nature communications
In common: SciPy, NumPy, extracellular electrophysiology (units, LFP), non-human primate, 5 references
[7] doi:10.1002/advs.202519893 [code]
NeuroSuite for Long-Term Functional and Structural Studies of Air-Liquid Interface Cerebral Organoids.
Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany)
In common: Pillow, pandas, SciPy, 2 other tools, extracellular electrophysiology (units, LFP), 3 references
[8] doi:10.1038/s41598-026-52253-9 [code]
Sequential visual stimuli increase high frequency power in the visual cortex.
Journal: Scientific reports
In common: h5py, Pillow, scikit-learn, 4 other tools, 2 references
[9] doi:10.7554/elife.103046 [code]
Dichotomy between extracellular signatures of active dendritic chemical synapses and gap junctions.
Journal: eLife
In common: Matplotlib, NumPy, extracellular electrophysiology (units, LFP), 5 references
[10] doi:10.1038/s41467-026-76109-y [code]
Assistive algorithms influence neural representations in motor brain-computer interfaces.
Journal: Nature communications
In common: scikit-learn, pandas, SciPy, 2 other tools, non-human primate, 3 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.