OSCR

A critical initialization for biological neural networks.

Code ↔ Paper

13 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 13 matches
  1. [1] § Methods › Data acquisition › Neuropixel recordings and processing ↔ fig2/analysis.py, lines 269–325 · score 0.85 · hippocampal formation, subcortical areas, firing rate, visual cortex, cortices, striatum
  2. [2] § Intrinsic structure in neural recordings ↔ fig2/fig2.py, lines 10–72 · score 0.82 · 2385–10344 regions, 2961–8566 ROIs, 1716–2914, 2385 regions, 2961 ROIs, Neuropixels
  3. [3] § Intrinsic structure in neural recordings ↔ fig2/fig2.py, lines 10–72 · score 0.82 · 2385–10344 regions, 2961–8566 ROIs, 1716–2914, 2385 regions, 2961 ROIs, Neuropixels
  4. [4] § Computations with symmetric dynamics ↔ fig5/fig5.ipynb, lines 61–163 · score 0.74 · zero shot working, persistent activity, working memory, symmetric dynamics, subspace, training
  5. [5] § Computations with symmetric dynamics ↔ fig5/fig5.ipynb, lines 61–163 · score 0.69 · zero shot working, working memory, symmetric dynamics, binary, training, models
  6. [6] § Methods › Data analysis › Estimating rotational components from data ↔ fig3/other_datasets.py, lines 405–474 · score 0.65 · linear track, virtual reality, visual cortex, reward, maze, PSTHs
  7. [7] § Methods › Data analysis › Estimating rotational components from data ↔ fig3/other_datasets.py, lines 405–474 · score 0.65 · linear track, virtual reality, visual cortex, reward, maze, PSTHs
  8. [8] § Methods ↔ fig_utils.py, lines 1–27 · score 0.62 · Howard Hughes Medical
  9. [9] § Methods › Simulations of dynamical systems › Sparse/varied connectivity ↔ fig4/sparse_clustered_local_sims.py, lines 136–193 · score 0.58 · exponential decay, locally connected, Bernoulli, distances, global, Sparse
  10. [10] § Methods › Simulations of dynamical systems ↔ simulations.py, lines 241–319 · score 0.54 · shot noise, relu, SNR, Poisson, decay, simulate
  11. [11] § Intrinsic structure in neural recordings ↔ fig2/fig2.py, lines 413–470 · score 0.54 · brainwide ephys, power law exponents, duration, SVCA2, CA1, neurons
  12. [12] § Intrinsic structure in neural recordings ↔ fig2/fig2.py, lines 413–470 · score 0.54 · brainwide ephys, power law exponents, duration, SVCA2, CA1, neurons
  13. [13] § Methods › Simulations of dynamical systems › Sparse/varied connectivity ↔ fig4/sparse_clustered_local_sims.py, lines 136–193 · score 0.53 · global connection, local connection, Bernoulli, clustered, Sparse, symmetric

Paper

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

The paper is loaded when this pane is shown.

The authors' code

Python · 973 lines · 43 KB · GPL-3.0 · 2 matches

  1. from matplotlib.patches import Ellipse
  2. from fig_utils import *
  3. from scipy.stats import zscore
  4. from powerlaw import fit_powerlaw_exp
  5. from scipy.stats import ttest_rel
  6. from matplotlib import gridspec
  7. from scipy.interpolate import interp1d
  10. def fig2(dat):
  11. areas_all = dat["areas_all"]
  12. evals_all = dat["evals_svca2_all"]
  13. evals_shuff_all = dat["evals_shuff_all"]
  14. evals_run = dat["evals_run_all"]
  16. fig = plt.figure(figsize=(14, 7), dpi=300)
  17. yratio = 14/7
  18. grid = plt.GridSpec(8, 7, wspace=0.4, hspace=0.2, figure=fig,
  19. bottom=0.04, top=0.94, left=0.02, right=1.)
  20. titles = ["cortical 2P imaging", "CA1 2P imaging",
  21. "brainwide Neuropixels"]
  22. sneur = ["2,385 - 10,344 ROIs", "2,961 - 8,566 ROIs", "1,716 - 2,914 units"]
  23. colors = dcolors[:3].copy()
  24. n_dset = len(areas_all)
  25. ids = np.zeros(n_dset, "int")
  26. for d in range(len(areas)):
  27. ids[np.array(areas_all)==areas[d]] = d
  28. ls = ["-", "-", "-", "--", "-."]
  29. print(areas_all, ids)
  30. dy = [0, 2, 1]
  31. alphas = np.zeros(n_dset)
  32. alphas_shuff = np.zeros(n_dset)
  33. ymax = 500
  34. norms = np.zeros(n_dset)
  36. if "imexs" in dat:
  37. # example images/traces (extracted from suite2p folders from example data)
  38. # (if available make example panels)
  39. imexs = dat["imexs"]
  40. masksexs = dat["masksexs"]
  42. for d in range(3):
  43. ax = plt.subplot(grid[:2, d])
  44. pos = ax.get_position().bounds
  45. ax.set_position([pos[0]+0.04*d + 0.01, pos[1], pos[2]*(1.3 +0.4*(d<2)), pos[3]*1.])
  46. if d==2:
  47. im = plt.imread("allenprobes.png")
  48. #print(im.shape, imexs[1].shape)
  49. ax.imshow(im[58:-100])
  50. else:
  51. ax.imshow(imexs[d], vmin=0, vmax=0.9,
  52. cmap="gray", aspect=0.75/0.5 if d==1 else 1)
  53. masks0 = masksexs[d].copy()
  54. yo, xo = np.nonzero(masks0[:,:,-1]>0)
  55. masks0[yo,xo,-1] = 0.25
  56. ax.imshow(masks0, aspect=0.75/0.5 if d==1 else 1)
  57. if d==0:
  58. ax.set_ylim([380, 380+90])
  59. ax.set_xlim(50, 50 + 130)
  60. elif d==1:
  61. #ax.set_ylim([350, 430])
  62. #ax.set_xlim([100, 100 + 80*0.75/0.5])
  63. ax.set_ylim([50, 50+90])
  64. ax.set_xlim([250, 250+130*0.75/0.5])
  65. ax.set_title(titles[d], color=colors[d], fontstyle="italic",
  66. y=1.075, loc="left", x=-0.05)#, fontweight="bold")
  67. ax.text(0.5, 1.025, sneur[d], transform=ax.transAxes,
  68. fontsize="small", ha="center")
  69. ax.axis("off")
  70. il = d
  71. transl = mtransforms.ScaledTranslation(-25 / 72, 12 / 72, fig.dpi_scale_trans)
  72. il = plot_label(ltr, il, ax, transl)
  74. for j in range(2):
  75. Xexs = dat["Xexs"] if j==0 else dat["Xemb_ex"]
  76. Xexs = [Xexs[0], Xexs[1]] if j==0 else Xexs
  77. if j==0:
  78. titles = ['Rastermap of 2p data - V1', 'CA1']
  79. tmins = [5800, 42000//7+200] # 5000
  80. tlen = 22*60*3//7
  81. else:
  82. titles = ['Simulated 2p data - max eigenvalue = 0.998', 'max eigenvalue = 0.975']
  83. tmins = [5960*2, 0]
  84. tlen = 22*60*3//7
  85. for d in range(2):
  86. ax = plt.subplot(grid[2+3*j:5+3*j, 2*d:2*d+2])
  87. pos = ax.get_position().bounds
  88. ax.set_position([pos[0]-d*0.015, pos[1]-j*0.02, pos[2]+0.005, pos[3]*0.9])
  89. #Xe = zscore(Xexs[d][:, tmins[d]:tmins[d]+tlen].copy(), axis=1)
  90. im = ax.imshow(zscore(Xexs[d][:, tmins[d]:tmins[d]+tlen], axis=1),
  91. aspect="auto", cmap="gray_r", vmin=0., vmax=1.5)
  92. ax.plot([0, 22*30/7], -Xexs[d].shape[0]*0.04*np.ones(2), color="k", lw=1.5)
  93. ax.plot(-0.025*tlen*np.ones(2), [0, 50], color="k", lw=1.5)
  94. if d==0:
  95. if j==0:
  96. ax.text((22*30/2)/7, -50, "30 sec.", ha="center", va="center", fontsize="small")
  97. ax.text(-0.05*tlen, 0, "1000 neurons", ha="center", va="bottom", fontsize="small", rotation=90)
  98. cax = ax.inset_axes([0.85, -0.01, 0.15, 0.025])
  99. cb = plt.colorbar(im, cax=cax, orientation="horizontal")
  100. cb.set_ticks([0, 1.0])
  101. cb.set_ticklabels(["0", "1.0"], fontsize="small")
  102. ax.text(0.83, -0.005, "z-scored activity", fontsize="small",
  103. ha="right", va="top", transform=ax.transAxes)
  104. ax.set_ylim([-Xexs[d].shape[0]*0.05, Xexs[d].shape[0]+0.5])
  105. ax.set_xlim([-0.028*tlen, tlen])
  106. ax.axis("off")
  107. ax.set_title(titles[d], fontsize='medium')
  108. if d==0:
  109. il = 3 + 3*j
  110. transl = mtransforms.ScaledTranslation(-10 / 72, 0 / 72, fig.dpi_scale_trans)
  111. il = plot_label(ltr, il, ax, transl)
  113. grid1 = gridspec.GridSpecFromSubplotSpec(3, 1, subplot_spec=grid[:5, 4], wspace=0., hspace=0.3)
  114. for d in range(3):
  115. ax = fig.add_subplot(grid1[d,0])
  116. pos = ax.get_position().bounds
  117. ax.set_position([pos[0]+0.01, pos[1]+(pos[3]-pos[2]*yratio)/2+0.01,
  118. pos[2]*0.8, pos[2]*yratio*0.8])
  119. pose = ax.get_position().bounds
  120. ix = np.nonzero(ids==d)[0]
  121. if d==0:
  122. ix = np.hstack((ix, np.nonzero(ids==4)[0]))
  123. ix = np.hstack((ix, np.nonzero(ids==3)[0]))
  124. lns = []
  125. for i in ix:
  126. ss = evals_all[i][:1000].copy()
  127. alphas[i], yp = fit_powerlaw_exp(evals_all[i],
  128. np.arange(10, ymax))
  129. ss /= yp[0]
  130. norms[i] = yp[0] #yp[10] * 10
  131. ax.loglog(np.arange(1, min(len(evals_all[i])+1, 1001)),
  132. ss, color=colors[d], lw=0.5, alpha=0.75, ls=ls[ids[i]])
  133. alphas_shuff[i] = fit_powerlaw_exp(evals_shuff_all[i],
  134. np.arange(10, ymax))[0]
  135. ln = ax.plot([], [], color=colors[d], lw=1, ls=ls[ids[i]])
  136. lns.append(ln)
  137. if d==0:
  138. ax.legend([lns[0][0], lns[-5][0], lns[-1][0]], ["V1", "sensori-\nmotor", "PPC"],
  139. frameon=False, loc="upper left", bbox_to_anchor=(0.3, 1.3),
  140. handlelength=1.2)
  142. ax.set_ylim(0.003, 3)
  143. ax.set_xlim(1, 1000)
  144. ax.set_xticks([1, 10, 100, 1000])
  145. ax.set_xticklabels(["1", "10", "100", "1,000"])
  146. ax.set_yticks([0.01, 0.1 ,1])
  147. ax.set_yticklabels(["0.01", "0.1", "1"])
  148. aexp = [-0.69, -1.254]
  149. yc = 1
  150. ax.fill_between([1, 1000], [yc, yc * 1000**aexp[0]], [yc, yc * 1000**aexp[1]],
  151. color="k", lw=0, alpha=0.1)
  152. if d==0:
  153. il = 4
  154. transl_e = mtransforms.ScaledTranslation(-50 / 72, 12 / 72, fig.dpi_scale_trans)
  155. il = plot_label(ltr, il, ax, transl_e)
  156. ax.set_ylabel("normalized variance")
  158. axin = ax.inset_axes([1.2, 0.5, 0.7, 0.7])
  159. axin.fill_between([1, 1000], [yc, yc * 1000**aexp[0]], [yc, yc * 1000**aexp[1]],
  160. color="k", lw=0, alpha=0.1)
  161. axin.text(0.27, 0.37, "symmetric", transform=axin.transAxes, fontsize="small",
  162. color=0.*np.ones(3), rotation=-(0.38)*90, fontstyle="italic")
  163. axin.text(0.08, 0.01, "non-symmetric", transform=axin.transAxes, fontsize="small",
  164. color=0.*np.ones(3), rotation=-(0.58)*90, fontstyle="italic")
  165. axin.set_xlim(1, 1000)
  166. axin.set_ylim(0.001, 1)
  167. axin.set_xscale("log")
  168. axin.set_yscale("log")
  169. axin.set_xticks([]); axin.set_yticks([])
  170. axin.minorticks_off()
  171. elif d==2:
  172. ax.set_xlabel("PC index")
  174. ax = plt.subplot(grid[:5, -2:])
  175. pos = ax.get_position().bounds
  176. ax.set_position([pos[0]+0.03, pos[1]+(pos[3]-pos[2]*yratio)/2-0.03,
  177. pos[2]*0.8, pos[2]*yratio*0.8])
  178. yh = [0, 2, 1]
  179. for d in range(3):
  180. ix = np.nonzero(ids==d)[0]
  181. emax = min(1000, np.array([len(evals_all[i]) for i in ix]).min())
  182. ev_all = np.array([evals_all[i][:emax] / norms[i] for i in ix])
  183. ev_mean = ev_all.mean(axis=0)
  184. ev_std = ev_all.std(axis=0) #/ np.sqrt(len(ix)-1)
  185. ax.loglog(np.arange(1, len(ev_mean)+1), ev_mean, color=colors[d], lw=1)
  186. ax.fill_between(np.arange(1, len(ev_mean)+1), ev_mean-ev_std, ev_mean+ev_std,
  187. color=colors[d], alpha=0.5, lw=0)
  188. ax.text(0.05, 0.1 + 0.12*yh[d], r"$\alpha$" + f" = {alphas[ix].mean():.2f}",
  189. transform=ax.transAxes, color=colors[d], fontsize="large")
  190. ax.set_ylim(0.003, 3)
  191. ax.set_xlim(1, 1000)
  192. ax.set_xticks([1, 10, 100, 1000])
  193. ax.set_xticklabels(["1", "10", "100", "1,000"])
  194. ax.set_yticks([0.01, 0.1 ,1])
  195. ax.set_yticklabels(["0.01", "0.1", "1"])
  196. ax.set_xlabel("PC index")
  197. ax.set_ylabel("normalized variance")
  198. ax.set_title("average", y=1, loc='left')
  199. il = 5
  200. transl = mtransforms.ScaledTranslation(-40 / 72, 0 / 72, fig.dpi_scale_trans)
  201. il = plot_label(ltr, il, ax, transl)
  203. ax = ax.inset_axes([0.7, 0.9, 0.3, 0.4])
  204. for d in range(3):
  205. ix = np.nonzero(ids==d)[0]
  206. if d==0:
  207. ix = np.hstack((ix, np.nonzero(ids==4)[0]))
  208. ix = np.hstack((ix, np.nonzero(ids==3)[0]))
  209. ash = np.stack((alphas[ix], alphas_shuff[ix]), axis=0)
  210. p = ttest_rel(ash[0], ash[1]).pvalue
  211. for i in ix:
  212. ax.plot(np.arange(0,2) + dy[d] * 1.5, (alphas[i], alphas_shuff[i]),
  213. color=colors[d], lw=1., alpha=0.75, ls=ls[ids[i]])
  214. star = "***" if p < 0.001 else "**" if p < 0.01 else "*" if p < 0.05 else "n.s."
  215. print(p)
  216. ax.text(dy[d]*1.5 + 0.5 + (dy[d]-1)*0., 0.9, f"{star}", ha="center",
  217. va="center", color=colors[d], fontsize="small" if p>=0.05 else "medium")
  218. ax.set_ylim([0., 1.0])
  219. ax.set_yticks([0., 0.5, 1.0])
  220. ax.set_xticks([0, 1])
  221. ax.set_xticklabels(["original", "shuffled"], rotation=90, ha="center", va="top")
  222. ax.set_ylabel("power-law\nexponent ($\\alpha$)")
  223. #il += 1
  224. #il = plot_label(ltr, il, ax, transl)
  226. evals_all = dat['evals_sub'].copy()
  227. alphas_sim = np.zeros((2, evals_all.shape[1])) * np.nan
  228. ax = plt.subplot(grid[-2:, 4])
  229. pos = ax.get_position().bounds
  230. ypos = pos[1]+(pos[3]-pos[2]*yratio)/2+0.03
  231. ax.set_position([pose[0], ypos, pos[2]*0.9, pos[2]*yratio*0.9])
  232. for d in range(2):
  233. for i in range(len(evals_all[0])):
  234. ss = evals_all[d, i].copy()
  235. alphas_sim[d, i], yp = fit_powerlaw_exp(ss, np.arange(10, ymax))
  236. ss /= yp[0]
  237. ax.loglog(np.arange(1, min(len(ss)+1, 1001)),
  238. ss[:1000], color=colors[d], lw=0.5, alpha=0.75)
  239. ax.text(0.35, 0.65 + 0.15*(1-d), r"$\alpha$" + f" = {alphas_sim[d].mean():.2f}",
  240. transform=ax.transAxes, color=colors[d], fontsize="medium")
  242. # ax.fill_between(np.arange(1, len(ev_mean)+1), ev_mean-ev_std, ev_mean+ev_std,
  243. # color=colors[d], alpha=0.5, lw=0)
  244. print(alphas_sim.mean(axis=-1))
  245. ax.set_title('simulated 2p data', fontsize='medium')
  246. ax.set_ylim(0.003, 3)
  247. ax.set_xlim(1, 1000)
  248. ax.set_xticks([1, 10, 100, 1000])
  249. ax.set_xticklabels(["1", "10", "100", "1,000"])
  250. ax.set_yticks([0.01, 0.1 ,1])
  251. ax.set_yticklabels(["0.01", "0.1", "1"])
  252. ax.set_xlabel("PC index")
  253. ax.set_ylabel("normalized variance")
  254. il += 1
  255. il = plot_label(ltr, il, ax, transl_e)
  257. transl = mtransforms.ScaledTranslation(-50 / 72, 0 / 72, fig.dpi_scale_trans)
  258. evals_enorms = dat['evals_enorms']
  259. enorms = dat['enorms']
  260. nnorm = len(enorms)
  261. n_sim = len(evals_enorms[0])
  262. alphas_enorms = np.zeros((nnorm, n_sim)) * np.nan
  263. ecolor = ['g', [0, 1, 0]]
  264. for d in range(nnorm):
  265. for i in range(len(evals_enorms[0])):
  266. alphas_enorms[d, i], yp = fit_powerlaw_exp(evals_enorms[d, i], np.arange(10, ymax))
  267. ax = plt.subplot(grid[-2:, -2])
  268. pos = ax.get_position().bounds
  269. ax.errorbar(enorms, alphas_enorms.mean(axis=1), alphas_enorms.std(axis=1)/(n_sim-1)**0.5, color=ecolor[0])
  270. ax1 = ax.twinx()
  271. ax1.set_position([pos[0] + pos[2]*0.25, ypos,
  272. pos[2]*0.5, pos[2]*yratio])
  273. ax1.plot(enorms, 0.5/(1-enorms) * 0.02, color=ecolor[1])
  274. ax1.spines['right'].set_visible(True)
  275. ax1.spines['right'].set_color(ecolor[1])
  276. ax1.tick_params(axis='y', colors=ecolor[1])
  277. ax1.set_ylabel('timescale (sec.)', color=ecolor[1], rotation=-90, va='bottom')
  278. ax1.set_yticks([0, 2, 4])
  279. ax.spines['left'].set_color(ecolor[0])
  280. ax.tick_params(axis='y', colors=ecolor[0])
  281. ax.set_xticks([0.975, 0.985, 0.998])
  282. ax.set_xlim([0.975, 0.998])
  283. ax.set_yticks([0.5, 0.6, 0.7])
  284. ax.set_ylabel('power-law exponent ($\\alpha$)', color=ecolor[0])
  285. ax.set_position([pos[0] + pos[2]*0.6, ypos,
  286. pos[2]*1.25, pos[2]*yratio])
  287. il = plot_label(ltr, il, ax, transl)
  289. return fig, alphas
  292. def suppfig_svca2(dat, dat_subsample, dat_janelia, dat_london):
  293. fig = plt.figure(figsize=(14, 10))
  294. yratio = 14 / 10
  295. grid = plt.GridSpec(5, 5, wspace=0.4, hspace=0.6, figure=fig,
  296. bottom=0.05, top=0.92, left=0.05, right=0.95)
  297. ix = 5
  298. il = 0
  299. dy = 0.
  301. shapes = dat['shapes']
  302. titles = ["eigenspectrum - direct", "SVCA", "SVCA2"]
  303. ticks = ["direct", "SVCA", "SVCA2"]
  305. evals_all_list = [dat["evals_all"], dat["evals_svca_all"], dat["evals_svca2_all"]]
  306. areas_all = dat["areas_all"]
  307. n_dset = len(areas_all)
  308. ids = np.zeros(n_dset, "int")
  309. alphas_all = np.zeros((n_dset, 3))
  310. for d in range(5):
  311. ids[np.array(areas_all)==areas[d]] = d
  313. colors = dcolors[:3].copy()
  314. colors = np.vstack((colors, colors[0], colors[0]))
  315. ls = ["-", "-", "-", "--", "-."]
  316. anames = ['2P cortex', '2P CA1', 'brainwide\nephys']
  318. grid1 = gridspec.GridSpecFromSubplotSpec(3, 3, subplot_spec=grid[:3, :3], wspace=0.2, hspace=0.4)
  320. for d in range(3):
  321. for k in range(3):
  322. ax = plt.subplot(grid1[d, k])
  323. pos = ax.get_position().bounds
  324. ax.set_position([pos[0] + 0.02*(1-k), pos[1], pos[3]/yratio, pos[3]])
  325. ix = (ids==0) + (ids==3) + (ids==4) if d==0 else ids==d
  326. ix = np.nonzero(ix)[0]
  327. for i in ix:
  328. evals = evals_all_list[k][i].copy()
  329. alphas_all[i, k], yp = fit_powerlaw_exp(evals, np.arange(10, 500))
  330. evals /= yp[0]
  331. ax.loglog(np.arange(1, len(evals)+1), evals, color=colors[ids[i]],
  332. lw=1, alpha=0.25 if d!=3 else 0.75, ls=ls[ids[i]], zorder=-30*(ids[i]==1) + 20*(ids[i]!=1))
  334. ax.minorticks_on()
  335. ax.set_ylim(0.001, 3)
  336. ax.set_xlim(1, 3000)
  337. ax.set_yticks([0.01, 0.1, 1])
  338. ax.set_yticklabels(["0.01", "0.1", "1"], fontsize='small')
  339. ax.set_xticks([1, 10, 100, 1000])
  340. ax.set_xticklabels(['1', '10 ', '100 ', ' 1,000'], fontsize='small')
  341. ax.text(0.3, 0.85, "$\\alpha = $" + f"{alphas_all[ix, k].mean():.2f}",
  342. color=colors[d], transform=ax.transAxes, fontweight="bold")
  343. ax.xaxis.set_minor_locator(matplotlib.ticker.LogLocator(base=10, subs=np.arange(2,10), numticks=10))
  344. if k==0:
  345. ax.text(0.05, 0.05, anames[d], color=colors[d],
  346. transform=ax.transAxes, fontweight='bold', fontstyle='italic')
  347. transl = mtransforms.ScaledTranslation(-40 / 72, 5 / 72, fig.dpi_scale_trans)
  348. il = plot_label(ltr, il, ax, transl)
  349. if d==0:
  350. if k==0:
  351. ax.set_xlabel("PC index")
  352. ax.set_ylabel("normalized variance")
  353. #ax.set_title("neural recordings", loc="left", x=-0.1, fontweight="bold")
  355. ax.set_title(titles[k], y=1.1, fontstyle='italic')
  358. ax = plt.subplot(grid[:3, 3])
  359. pos = ax.get_position().bounds
  360. ax.set_position([pos[0]- 0.03, pos[1]+0.25*pos[3], pos[2], pos[3]*0.8])
  361. xp = np.arange(3)*np.ones((len(evals_all_list[0]),1))
  362. xp += np.random.randn(*xp.shape)*0.05
  363. cols = np.tile(colors[ids][:,np.newaxis], (1,3)).reshape(-1,3)
  364. ax.scatter(xp.flatten(), alphas_all.flatten(), color=cols, s=10)
  365. for d in range(3):
  366. ix = (ids==0) + (ids==4) + (ids==5) if d==0 else ids==d
  367. ax.scatter(np.arange(3), alphas_all[ix].mean(axis=0), color=colors[d],
  368. s=400, marker="_")
  369. ax.set_ylabel("power-law exponent ($\\alpha$)")
  370. ax.set_xticks(np.arange(3))
  371. ax.set_xticklabels(["direct", "SVCA", "SVCA2"], rotation=45, ha='right')
  372. transl = mtransforms.ScaledTranslation(-50 / 72, -5 / 72, fig.dpi_scale_trans)
  373. il = plot_label(ltr, il, ax, transl)
  374. ax.set_ylim([0.2, 1.4])
  376. ax = plt.subplot(grid[:3, 4])
  377. cols = [colors[0], colors[0], [0.8, 0, 0], [0.5, 0, 0.5], [1, 0, 1]]
  378. pos = ax.get_position().bounds
  379. ax.set_position([pos[0]-0.01, pos[1]+0.25*pos[3], pos[2]*1.25, pos[3]*0.8])
  380. ix = (ids==0) + (ids==4) + (ids==5)
  381. for j in range(5):
  382. if j<2:
  383. alphas = alphas_all[ix, (1-j)+1]
  384. elif j==2:
  385. alphas = dat_subsample['alphas_all']
  386. elif j==3:
  387. alphas = dat_janelia['alphas_all']
  388. elif j==4:
  389. alphas = dat_london['alphas_all']
  390. xp = j*np.ones(len(alphas))
  391. xp += np.random.randn(*xp.shape)*(0.05 + 0.05*(j==4))
  392. ax.scatter(xp, alphas, color=cols[j], s=10)
  393. ax.scatter(j, alphas.mean(), color=cols[j], s=500, marker="_")
  394. ax.set_ylabel("power-law exponent ($\\alpha$)")
  395. ax.set_ylim([0.2, 1.4])
  396. ax.set_xticks(np.arange(5))
  397. ax.set_xticklabels(["SVCA2, GCaMP8s 22Hz", "SVCA, GCaMP8s 22Hz",
  398. "SVCA, GCaMP8s 3Hz", "SVCA, GCaMP6s 3Hz\n(Janelia)",
  399. "SVCA, GCaMP6s 3Hz\n(London)"], rotation=45, ha='right')
  400. for idx, lbl in enumerate(ax.get_xticklabels()):
  401. lbl.set_color(cols[idx])
  402. transl = mtransforms.ScaledTranslation(-50 / 72, -5 / 72, fig.dpi_scale_trans)
  403. il = plot_label(ltr, il, ax, transl)
  406. shapes = dat['shapes']
  407. areas_all = dat['areas_all']
  408. n_dset = len(shapes)
  409. ids = np.zeros(n_dset, "int")
  410. for d in range(5):
  411. ids[np.array(areas_all)==areas[d]] = d
  413. anames = ['2P cortex', '2P CA1', 'brainwide ephys']
  414. scolors = ['g', 'y', 'b']
  415. bsize = 4000
  416. ikeeps = np.arange(50, bsize+1, 200)
  417. ikeeps[-1] = 4000
  418. ineurs = np.arange(0.01, 1.05, 0.05)
  419. ineurs[-1] = 1.0
  420. nneurons = ineurs * np.array(shapes)[:,:1]
  421. ntimes = ikeeps / bsize * np.array(shapes)[:,1:] / (22*60)
  422. xticks = [[[0, 2500, 5000], [0, 50, 100]],
  423. [[0, 2500, 5000], [0, 50, 100]],
  424. [[0, 1000, 2000], [0, 15, 30]]]
  425. fstr = ['', 'svca_', 'svca2_']
  426. grid1 = gridspec.GridSpecFromSubplotSpec(1, 6, subplot_spec=grid[-2:, :],
  427. wspace=0.4, hspace=0.4)
  428. axs = [[plt.subplot(grid1[-2:, 2*d]) for d in range(3)],
  429. [plt.subplot(grid1[-2:, 2*d+1]) for d in range(3)]]
  430. for k in range(3):
  431. evals_all = [dat[f'evals_{fstr[k]}all_times'].copy(), dat[f'evals_{fstr[k]}all_neurons'].copy()]
  432. for j in range(2):
  433. nvar = len(evals_all[j][0])
  434. alphas_all = np.zeros((n_dset, nvar))*np.nan
  435. for i in range(len(evals_all[j])):
  436. for n in range(len(evals_all[j][i])):
  437. evals = evals_all[j][i][n].copy()
  438. ymax = min(len(evals)//2, 500) #500 if len(evals) > 500 else min(len(evals), 100)
  439. if ymax > 15:
  440. alphas_all[i,n], yp = fit_powerlaw_exp(evals, np.arange(10, ymax))
  442. for d in range(3):
  443. ax = axs[j][d]
  444. ix = (ids==0) + (ids==4) + (ids==5) if d==0 else ids==d
  445. nn = nneurons[ix].mean(axis=0) if j==1 else ntimes[ix].mean(axis=0)
  446. ax.errorbar(nn, np.nanmean(alphas_all[ix], axis=0),
  447. np.nanstd(alphas_all[ix], axis=0) / ((~np.isnan(alphas_all[ix])).sum(axis=0)-1)**0.5,
  448. color=scolors[k], lw=1)
  449. ax.set_ylim([0.2, 1.5])
  450. ax.tick_params(axis='both', labelsize='small')
  451. ax.set_xticks(xticks[d][1-j])
  452. ax.set_yticks([0.5, 1, 1.5])
  453. ax.set_xlabel('# of neurons' if j==1 else 'duration (min.)')
  454. if k==0 and j==0:
  455. ax.set_ylabel('power-law exponent ($\\alpha$)')
  456. transl = mtransforms.ScaledTranslation(-55 / 72, -5 / 72, fig.dpi_scale_trans)
  457. il = plot_label(ltr, il, ax, transl)
  458. ax.set_title(anames[d], loc='left', color=colors[d],
  459. fontstyle='italic', fontweight='bold')
  461. if k==0:
  462. pos = ax.get_position().bounds
  463. ax.set_position([pos[0]+0.015*(j==0)*(k==0), pos[1], pos[2], pos[3]*0.9])
  465. if j==0 and d==0:
  466. ax.text(1, 0.9-0.1*k, ['direct', 'SVCA', 'SVCA2'][k], color=scolors[k],
  467. ha='right', transform=ax.transAxes)
  470. return fig
  473. def suppfig_ephys(dat, dat_areas, dat_tbins):
  474. fig = plt.figure(figsize=(8, 9))
  475. yratio = 8 / 9
  476. grid = plt.GridSpec(7, 4, wspace=0.45, hspace=0.6, figure=fig,
  477. bottom=0.07, top=0.98, left=0.1, right=0.95)
  478. ix = 5
  479. il = 0
  480. dy = 0.
  482. shapes = dat['shapes']
  484. areas_all = dat["areas_all"]
  485. n_dset = len(areas_all)
  486. ids = np.zeros(n_dset, "int")
  487. alphas_all = np.zeros((n_dset, 3))
  488. for d in range(5):
  489. ids[np.array(areas_all)==areas[d]] = d
  492. ax = plt.subplot(grid[:3, :])
  493. pos = ax.get_position().bounds
  494. ax.set_position([pos[0]-0.05, pos[1]-0.2*pos[3], pos[2]+0.065, pos[3]*1.08])
  495. tmin = 1120#42000//7 - 0
  496. tlen = 22*60*3//7
  497. Xex = zscore(dat['Xexs'][2][:, tmin:tmin+tlen], axis=1)
  498. im = ax.imshow(Xex, aspect="auto", cmap="gray_r", vmin=0., vmax=1.25)
  499. ax.plot([0, 22*30/7], -Xex.shape[0]*0.03*np.ones(2), color="k", lw=1.5)
  500. ax.plot(-0.015*tlen*np.ones(2), [0, 50], color="k", lw=1.5)
  501. ax.text((22*30/2)/7, -0.07*Xex.shape[0], "30 sec.", ha="center", va="center", fontsize="small")
  502. ax.text(-0.03*tlen, 0, "500 neurons", ha="center", va="bottom", fontsize="small", rotation=90)
  503. cax = ax.inset_axes([0.85, -0.01, 0.15, 0.025])
  504. cb = plt.colorbar(im, cax=cax, orientation="horizontal")
  505. cb.set_ticks([0, 1.0])
  506. cb.set_ticklabels(["0", "1.0"], fontsize="small")
  507. ax.text(0.83, -0.005, "z-scored activity", fontsize="small",
  508. ha="right", va="top", transform=ax.transAxes)
  509. ax.set_ylim([-Xex.shape[0]*0.05, Xex.shape[0]+0.5])
  510. ax.set_xlim([-0.028*tlen, tlen])
  511. ax.axis("off")
  512. ax.set_title('Rastermap of brainwide ephys activity')
  513. transl = mtransforms.ScaledTranslation(-20 / 72, 0 / 72, fig.dpi_scale_trans)
  514. il = plot_label(ltr, il, ax, transl)
  516. colors = dcolors[:3].copy()
  518. acolors = ['g', 'm', [0.25, 1, 0.25], 'r', [0.5, 0, 0]]
  520. ineurs = np.arange(0.01, 1.05, 0.05)
  521. ineurs[-1] = 1.0
  522. nneurons = ineurs * np.array(shapes)[:,:1]
  523. evals_all = dat[f'evals_svca2_all_neurons'].copy()
  524. nvar = len(evals_all[0])
  525. alphas_all = np.zeros((n_dset, nvar))*np.nan
  526. for i in range(len(evals_all)):
  527. for n in range(len(evals_all[i])):
  528. evals = evals_all[i][n].copy()
  529. ymax = min(len(evals)//2, 500) #500 if len(evals) > 500 else min(len(evals), 100)
  530. if ymax > 15:
  531. alphas_all[i,n], yp = fit_powerlaw_exp(evals, np.arange(10, ymax))
  533. nh = [4, 2, 3, 1, 0]
  534. area_names = dat_areas['area_groups'].keys()
  535. evals_all = dat_areas['evals_areas'].copy()
  536. nareas = len(area_names)
  537. alphas_areas = np.zeros((nareas, 3)) * np.nan
  538. nn_areas = np.stack([dat_areas['nneurons'][a] for a in area_names], axis=0)
  539. for i in range(3):
  540. ax = plt.subplot(grid[3:5, i+1])
  541. pos = ax.get_position().bounds
  542. ax.set_position([pos[0] + (2-i)*0.03 + 0.01, pos[1]-0.02, pos[2]*0.95, pos[2]*yratio*0.95])
  543. for n, area in enumerate(area_names):
  544. evals = evals_all[area][i].copy()
  545. if len(evals) > 0:
  546. ymax = min(len(evals)//2, 500) #500 if len(evals) > 500 else min(len(evals), 100)
  547. alphas_areas[n, i], yp = fit_powerlaw_exp(evals, np.arange(10, ymax))
  548. evals /= yp[0]
  549. ax.loglog(np.arange(1, len(evals)+1), evals, color=acolors[n], lw=1)
  550. nstr = '# of neurons ' if nh[n]==0 else ''
  551. nhn = nh[n] if i<2 else nh[n] - 1*(nh[n]>1)
  552. ax.text(1, 0.95-0.12*nhn, f'{nstr}= {int(nn_areas[n, i]):,d}', color=acolors[n],
  553. ha='right', transform=ax.transAxes, fontsize='small')
  554. ax.set_ylim(0.003, 3)
  555. ax.set_xlim(1, 1000)
  556. ax.set_yticks([0.01, 0.1, 1])
  557. ax.set_yticklabels(["0.01", "0.1", "1"], fontsize='small')
  558. ax.set_xticks([1, 10, 100, 1000])
  559. ax.set_xticklabels(['1', '10', '100', '1,000'], fontsize='small')
  560. ax.set_title(f'mouse {i+1}', fontsize='medium', y=1.1)
  561. ax.set_xlabel("PC index")
  562. if i==0:
  563. ax.set_ylabel("normalized variance")
  564. transl = mtransforms.ScaledTranslation(-50 / 72, 10 / 72, fig.dpi_scale_trans)
  565. il = 2
  566. il = plot_label(ltr, il, ax, transl)
  567. il -= 2
  570. d = 2
  571. ax = plt.subplot(grid[3:5, 0])
  572. pos = ax.get_position().bounds
  573. ax.set_position([pos[0], pos[1]-0.02, pos[2], pos[2]*yratio*1.1])
  574. ix = ids==d
  575. nn = nneurons[ix].mean(axis=0)
  576. ax.errorbar(nn, np.nanmean(alphas_all[ix], axis=0),
  577. np.nanstd(alphas_all[ix], axis=0) / ((~np.isnan(alphas_all[ix])).sum(axis=0)-1)**0.5,
  578. color=colors[d], lw=1)
  579. for n, area in enumerate(area_names):
  580. ax.scatter(nn_areas[n], alphas_areas[n], color=acolors[n],
  581. s=30, marker='x', zorder=30, alpha=0.9)
  582. ax.text(1.25, 1-0.12*nh[n], area, color=acolors[n], ha='right',
  583. transform=ax.transAxes, fontsize='small', va='bottom')
  584. ax.text(1.25, 0.02, u'\u2014 random\nsubsets', color=dcolors[d], ha='right',
  585. transform=ax.transAxes, fontweight='bold')
  586. ax.set_ylim([0.4, 1.8])
  587. ax.set_xscale('log')
  588. ax.set_xticks([100, 1000])
  589. ax.set_xticklabels(['100', '1,000'])
  590. ax.set_xlabel('# of neurons')
  591. ax.set_ylabel('power-law exponent ($\\alpha$)')
  594. transl = mtransforms.ScaledTranslation(-55 / 72, 0 / 72, fig.dpi_scale_trans)
  595. il = plot_label(ltr, il, ax, transl)
  596. il += 1
  598. tbins = dat_tbins['tbins']
  599. acolors = plt.get_cmap('YlOrBr_r')(np.linspace(0, 0.6, len(tbins)))
  600. evals_all = dat_tbins['evals_tbins']
  601. alphas_tbins = np.zeros((len(tbins), 3)) * np.nan
  602. for i in range(3):
  603. ax = plt.subplot(grid[5:7, i+1])
  604. pos = ax.get_position().bounds
  605. ax.set_position([pos[0] + (2-i)*0.03 + 0.01, pos[1], pos[2]*0.95, pos[2]*yratio*0.95])
  606. for n, tbin in enumerate(tbins):
  607. evals = evals_all[i][n].copy()
  608. if len(evals) > 0:
  609. ymax = min(len(evals)//2, 500) #500 if len(evals) > 500 else min(len(evals), 100)
  610. alphas_tbins[n, i], yp = fit_powerlaw_exp(evals, np.arange(10, ymax))
  611. evals /= yp[0]
  612. ax.loglog(np.arange(1, len(evals)+1), evals, color=acolors[n], lw=1)
  613. ax.set_ylim(0.003, 3)
  614. ax.set_xlim(1, 1000)
  615. ax.set_yticks([0.01, 0.1, 1])
  616. ax.set_yticklabels(["0.01", "0.1", "1"], fontsize='small')
  617. ax.set_xticks([1, 10, 100, 1000])
  618. ax.set_xticklabels(['1', '10', '100', '1,000'], fontsize='small')
  619. ax.set_title(f'mouse {i+1}', fontsize='medium')
  620. ax.set_xlabel("PC index")
  621. if i==0:
  622. ax.set_ylabel("normalized variance")
  623. transl = mtransforms.ScaledTranslation(-50 / 72, 10 / 72, fig.dpi_scale_trans)
  624. il += 1
  625. il = plot_label(ltr, il, ax, transl)
  626. il -= 2
  628. ax = plt.subplot(grid[5:7, 0])
  629. pos = ax.get_position().bounds
  630. ax.set_position([pos[0], pos[1], pos[2], pos[2]*yratio*1.1])
  631. ax.errorbar(tbins*1000, alphas_tbins.mean(axis=-1),
  632. alphas_tbins.std(axis=-1)/(2**0.5), color='k')
  633. for n, tbin in enumerate(tbins):
  634. ax.scatter(tbin*1000, alphas_tbins[n].mean(), s=30,
  635. color=acolors[n], zorder=30)
  636. ax.set_xscale('log')
  637. ax.set_xlabel('time bin (ms)')
  638. ax.set_ylabel('power-law exponent ($\\alpha$)')
  639. ax.set_xticks([10, 100])
  640. ax.set_xticklabels(['10', '100'])
  641. ax.set_ylim([0.65, 0.85])
  642. ax.set_yticks([0.7, 0.8])
  644. transl = mtransforms.ScaledTranslation(-58 / 72, 0 / 72, fig.dpi_scale_trans)
  645. il = plot_label(ltr, il, ax, transl)
  647. return fig
  650. def suppfig_running(dat):
  651. colors = [[0.7, 0, 0], [1, 0.5, 0.5]]
  652. fig = plt.figure(figsize=(14, 3), dpi=150)
  653. yratio = 14 / 3
  654. grid = plt.GridSpec(1, 5, wspace=0.4, hspace=0.2, figure=fig,
  655. bottom=0.18, top=0.9, left=0.05, right=0.95)
  656. transl = mtransforms.ScaledTranslation(-45 / 72, 5 / 72, fig.dpi_scale_trans)
  657. il = 0
  658. evals_run_all = dat['evals_run_all'].copy()
  659. nrun = evals_run_all.shape[1]
  660. alphas_all = np.zeros((2, nrun))
  661. ax = plt.subplot(grid[0, 0])
  662. pos = ax.get_position().bounds
  663. ax.set_position([pos[0], pos[1], pos[2], pos[2]*yratio])
  664. rstr = ['running', 'not running']
  665. for j in range(2):
  666. for i in range(nrun):
  667. alphas_all[j, i], yp = fit_powerlaw_exp(evals_run_all[j, i], np.arange(10, 500))
  668. evals_run_all[j, i] /= yp[0]
  669. evals = np.nanmean(evals_run_all[j], axis=0)
  670. evals_std = np.nanstd(evals_run_all[j], axis=0) / (nrun-1)**0.5
  671. ax.loglog(np.arange(1, len(evals)+1), evals, lw=1, color=colors[j])
  672. ax.fill_between(np.arange(1, len(evals)+1), evals - evals_std, evals + evals_std,
  673. edgecolor='none', facecolor=colors[j], alpha=0.25)
  674. alpha = fit_powerlaw_exp(evals, np.arange(10, 500))[0]
  675. ax.text(1, 1-0.1*j, f'{rstr[j]}, $\\alpha$={alpha:.2f}', ha='right',
  676. transform=ax.transAxes, color=colors[j])
  677. ax.set_yscale('log')
  678. ax.set_xscale('log')
  679. ax.set_ylim(0.003, 0.003*700)
  680. ax.set_xlim(1, 700)
  681. ax.set_yticks([0.01, 0.1, 1])
  682. ax.set_yticklabels(["0.01", "0.1", "1"], fontsize='small')
  683. ax.set_xticks([1, 10, 100])
  684. ax.set_xticklabels(['1', '10', '100'], fontsize='small')
  685. ax.set_xlabel("PC index")
  686. ax.set_ylabel("normalized variance")
  687. il = plot_label(ltr, il, ax, transl)
  689. ax = plt.subplot(grid[0, 1])
  690. pos = ax.get_position().bounds
  691. ax.set_position([pos[0] + pos[2]*0.1, pos[1], pos[2]*0.5, pos[2]*yratio])
  692. areas = dat['areas']
  693. ls = {'V1': '-', 'sensorimotor': '--', 'PPC': '-.'}
  694. for i in range(nrun):
  695. ax.plot(np.arange(0,2), alphas_all[:,i],
  696. color='k', lw=1., alpha=0.75, ls=ls[areas[i]])
  697. p = ttest_rel(alphas_all[0], alphas_all[1]).pvalue
  698. star = "***" if p < 0.001 else "**" if p < 0.01 else "*" if p < 0.05 else "n.s."
  699. print(p)
  700. ax.text(0.5, 1., f"{star}", ha="center",
  701. va="center", color='k', fontsize="small" if p>=0.05 else "medium")
  702. ax.set_ylim([0., 1.0])
  703. ax.set_yticks([0., 0.5, 1.0])
  704. ax.set_xticks([0, 1])
  705. ax.set_xticklabels(["running", "not\nrunning"], rotation=0, ha="center", va="top")
  706. ax.set_ylabel("power-law exponent ($\\alpha$)")
  707. lns = []
  708. area_names = ['V1', 'sensorimotor', 'PPC']
  709. for area in area_names:
  710. ln = ax.plot([], [], color='k', lw=1, ls=ls[area])
  711. lns.append(ln[0])
  712. ax.legend(lns, area_names,
  713. frameon=False, loc="lower left", #bbox_to_anchor=(0.2, 1.1),
  714. handlelength=1.2)
  715. il = plot_label(ltr, il, ax, transl)
  717. iex = 2
  718. for j in range(2):
  719. ax = plt.subplot(grid[0, 2+j])
  720. pos = ax.get_position().bounds
  721. ax.set_position([pos[0], pos[1], pos[2], pos[2]*yratio])
  722. e = dat['dmd_evals_all'][j, iex].copy()
  723. ix = np.abs(e)>.25
  724. iang = np.angle(e[ix]) / (2*np.pi)
  725. iabs = -np.log10(np.abs(e[ix]))
  726. irot = iang/iabs
  727. ixx = e[ix].imag>=0
  728. mu = irot[ixx].mean()
  729. sd = irot[ixx].std()
  730. m = np.percentile(irot[ixx], [5, 25, 75, 95])
  731. med = np.median(irot[ixx])
  732. print(m, med)
  733. ax.scatter(e.real, e.imag, s=10, color=colors[j])
  734. ax.set_ylim([-1, 1])
  735. ax.set_xlim([0.25, 1])
  736. ax.set_xlabel('real part')
  737. ax.set_ylabel('imaginary part')
  738. ax.text(0.5, 0.8, rstr[j], color=colors[j], ha='center', transform=ax.transAxes)
  739. if j==0:
  740. il = plot_label(ltr, il, ax, transl)
  741. ax.set_title('Eigenvalues of DMD matrix (dt=0.23s), example mouse', loc='left', y=1.02, fontsize='medium')
  743. ax = plt.subplot(grid[0, -1])
  744. pos = ax.get_position().bounds
  745. ax.set_position([pos[0], pos[1], pos[2], pos[2]*yratio])
  746. igood = np.nonzero(~np.isnan(dat['dmd_evals_all'][0,:,0]))[0]
  747. for i, iex in enumerate(igood):
  748. for j in range(2):
  749. e = dat['dmd_evals_all'][j, iex].copy()
  750. ix = np.abs(e)>.25
  751. iang = np.angle(e[ix]) / (2*np.pi)
  752. iabs = -np.log10(np.abs(e[ix]))
  753. irot = iang/iabs
  754. ixx = e[ix].imag>=0
  755. mu = irot[ixx].mean()
  756. sd = irot[ixx].std()
  757. m = np.percentile(irot[ixx], [5, 25, 50, 75, 95])
  758. med = np.median(irot[ixx])
  759. dy = 0.3
  760. yy = -(4*i - j)
  761. print(ixx.sum())
  762. ax.plot([m[0], m[-1]], -(4*i - j) * np.ones(2), lw=1, color=np.minimum(1, np.array(colors[j])+0.25))
  763. ax.plot([m[1], m[-2]], -(4*i - j) * np.ones(2), lw=4, color=np.minimum(1, np.array(colors[j])+0.25))
  764. ax.plot([m[2], m[2]], [yy-dy, yy+dy], lw=1, color=colors[j])
  765. if j==0:
  766. ax.text(-0.1, -4*i, f'mouse {len(igood) - i}', rotation=90, ha='center', va='center', fontsize='small')
  767. ax.set_xlim([-0.1, 1])
  768. ax.set_ylim([3, -9])
  769. ax.spines['left'].set_visible(False)
  770. ax.set_yticks([])
  771. ax.set_xlabel('rotations per 10-fold attenuation')
  772. il = plot_label(ltr, il, ax, transl)
  774. return fig
  779. # def suppfig_svca(dat):
  780. # fig = plt.figure(figsize=(14, 14), dpi=150)
  781. # yratio = 14 / 14
  782. # grid = plt.GridSpec(5, 6, wspace=0.4, hspace=0.2, figure=fig,
  783. # bottom=0.05, top=1, left=0.08, right=0.95)
  784. # ix = 5
  785. # transl = mtransforms.ScaledTranslation(-45 / 72, ix / 72, fig.dpi_scale_trans)
  786. # il = 0
  787. # dy = 0.
  789. # titles = ["eigenspectrum - direct", "SVCA", "SVCA2"]
  790. # ticks = ["direct", "SVCA", "SVCA2"]
  792. # evals_all_list = [dat["evals_all"], dat["evals_svca_all"], dat["evals_svca2_all"]]
  793. # areas_all = dat["areas_all"]
  794. # ids = np.zeros(len(areas_all), "int")
  795. # alphas_all = np.zeros((len(evals_all_list[0]), 3))
  796. # for d in range(5):
  797. # ids[np.array(areas_all)==areas[d]] = d
  799. # colors = dcolors[:3].copy()
  800. # colors.append(colors[0])
  801. # colors.append(colors[0])
  802. # ls = ["-", "-", "-", "--", "-."]
  804. # for d in range(5):
  805. # ax = plt.subplot(grid[-1, d])
  806. # pos = ax.get_position().bounds
  807. # ax.set_position([pos[0], pos[1], pos[2], pos[2]*yratio])
  808. # for i in range(len(evals_all_list[d])):
  809. # evals = evals_all_list[d][i].copy()
  810. # alphas_all[i, d], yp = fit_powerlaw_exp(evals, np.arange(10, 500))
  811. # evals /= yp[0]
  812. # ax.loglog(np.arange(1, len(evals)+1), evals, color=dcolors[ids[i]],
  813. # lw=1, alpha=1, ls=ls[d])
  814. # ax.set_ylim(0.001, 3)
  815. # ax.set_xlim(1, 3000)
  816. # ax.set_yticks([0.01, 0.1, 1])
  817. # ax.set_yticklabels(["0.01", "0.1", "1"])
  818. # ax.set_xticks([1, 10, 100, 1000])
  819. # ax.set_xticklabels(["1", "10", "100", "1,000"])
  820. # for k in range(3):
  821. # ax.text(0.05, 0.32 - 0.12*k, "$\\alpha = $" + f"{alphas_all[ids==k, d].mean():.2f}",
  822. # color=dcolors[k], transform=ax.transAxes, fontweight="bold")
  823. # ax.set_xlabel("PC index")
  824. # if d==0:
  825. # ax.set_ylabel("normalized variance")
  826. # il = plot_label(ltr, il, ax, transl)
  827. # ax.set_title("neural recordings", loc="left", x=-0.1, fontweight="bold")
  829. # ax = plt.subplot(grid[-1, 3])
  830. # pos = ax.get_position().bounds
  831. # ax.set_position([pos[0] + 0.025, pos[1]+dy-0.02, pos[2], pos[2]*yratio + 0.04])
  832. # xp = np.arange(3)*np.ones((len(evals_all_list[0]),1))
  833. # xp += np.random.randn(*xp.shape)*0.05
  834. # cols = np.tile(dcolors[ids][:,np.newaxis], (1,3)).reshape(-1,3)
  835. # ax.scatter(xp.flatten(), alphas_all.flatten(), color=cols, s=10)
  836. # for k in range(3):
  837. # ax.scatter(np.arange(3), alphas_all[ids==k].mean(axis=0), color=dcolors[k],
  838. # s=400, marker="_")
  839. # ax.set_ylabel("power-law exponent ($\\alpha$)")
  840. # ax.set_xticks(np.arange(3))
  841. # ax.set_xticklabels(["direct", "SVCA", "SVCA2"], rotation=0)
  842. # il = plot_label(ltr, il, ax, transl)
  844. # return fig
  846. # def suppfig_svca2_sizes(dat_sim, alphas=None):
  847. # evals_svca2_all = dat_sim["evals_svca2_all"]
  848. # nonsyms = dat_sim["nonsyms"]
  849. # nneurons = dat_sim["nneurons"]
  850. # ntimes = dat_sim["ntimes"].astype('float32') / 23 # convert to seconds
  851. # noise_levels = dat_sim["noise_levels"]
  852. # n_sim = len(evals_svca2_all)
  854. # if alphas is None:
  855. # alphas = np.zeros((evals_svca2_all.shape[:-1]))
  856. # for i in range(n_sim):
  857. # for ni, nonsym in enumerate(nonsyms):
  858. # for nl, noise_level in enumerate(noise_levels):
  859. # for ii, nneur in enumerate(nneurons):
  860. # for jj, ntime in enumerate(ntimes):
  861. # evals = evals_svca2_all[i,ni,nl,ii,jj].copy()
  862. # ymax = min((~np.isnan(evals)).sum(), int(nneur*0.4), int(ntime*0.4))
  863. # evals = evals_svca2_all[i,ni,nl,ii,jj].copy()
  864. # yrange = np.arange(10, min(500, (~np.isnan(evals)).sum()-1))
  865. # alpha = fit_powerlaw_exp(evals, yrange)[0]
  866. # alphas[i,ni,nl,ii,jj] = alpha
  869. # fig = plt.figure(figsize=(14, 9), dpi=150)
  870. # yratio = 14 / 9
  871. # grid = plt.GridSpec(4, 7, wspace=0.6, hspace=0.7, figure=fig,
  872. # bottom=0.05, top=0.95, left=0.06, right=0.98)
  873. # grid1 = gridspec.GridSpecFromSubplotSpec(4, 6, subplot_spec=grid[:, 2:], wspace=0.3, hspace=0.7)
  874. # il = 0
  876. # n_sim = len(evals_svca2_all)
  878. # for ni, nonsym in enumerate(nonsyms):
  879. # for j in range(2):
  880. # if j==0:
  881. # lcolors = plt.get_cmap("Oranges")(np.linspace(0.5, 1, len(nneurons)//2))
  882. # else:
  883. # lcolors = plt.get_cmap("Greens")(np.linspace(0.5, 1, len(ntimes)//2))
  885. # ax = plt.subplot(grid[ni, j])
  886. # pos = ax.get_position().bounds
  887. # ax.set_position([pos[0], pos[1], pos[2], pos[2]*yratio])
  888. # if ni==0:
  889. # transl = mtransforms.ScaledTranslation(-40 / 72, -4 / 72, fig.dpi_scale_trans)
  890. # else:
  891. # transl = mtransforms.ScaledTranslation(-40 / 72, -10 / 72, fig.dpi_scale_trans)
  892. # il = plot_label(ltr, il, ax, transl)
  894. # if j==0:
  895. # if ni==0:
  896. # x0, y0 = -0.65, 1.3
  897. # else:
  898. # x0, y0 = -0.65, 1.2
  899. # ax.text(x0, y0, ["symmetric", "1/3 non-symmetric", "2/3 non-symmetric", "non-symmetric"][ni] + " connectivity",
  900. # transform=ax.transAxes, fontsize="medium", fontweight="bold", fontstyle='italic')
  901. # if ni == 0:
  902. # ax.text(0.7, 1.05, '# neurons = ' if j==0 else 'time =',
  903. # color='k', transform=ax.transAxes, ha='right')
  905. # for ii, nn in enumerate(nneurons) if j ==0 else enumerate(ntimes):
  906. # if ii%2==1:
  907. # continue
  908. # evals = evals_svca2_all[0,ni,0,ii,-1].copy() if j==0 else evals_svca2_all[0,ni,0,-1,ii].copy()
  909. # ymax = min((~np.isnan(evals)).sum(), int(nn*0.4))
  910. # evals = evals[:ymax]
  911. # alpha = plot_spectrum(ax, evals, color=lcolors[ii//2], lw=1, plot_fit=False)
  912. # if ii%4==0 and ni==0:
  913. # if j==0:
  914. # ax.text(0.75, 0.65 + 0.1*ii/4, f"{nn:.0f}",
  915. # color=lcolors[ii//2], transform=ax.transAxes, fontsize="small")
  916. # else:
  917. # ax.text(0.75, 1.05 - 0.1*ii/4, f"{nn/60:.1f} min",
  918. # color=lcolors[ii//2], transform=ax.transAxes, fontsize="small")
  919. # if ni==0:
  920. # if j==0:
  921. # ax.set_xlabel("PC index")
  922. # ax.set_ylabel("normalized\nvariance")
  923. # else:
  924. # ax.set_xlabel("PC index")
  926. # for j in range(len(noise_levels)):
  927. # ax = plt.subplot(grid1[ni, j])
  928. # pos = ax.get_position().bounds
  929. # ax.set_position([pos[0]+0.008*(len(noise_levels) - j), pos[1], pos[2], pos[2]*yratio])
  930. # im = ax.imshow(alphas.mean(axis=0)[ni, j].T, vmin=0.5, vmax=2, cmap='viridis',
  931. # aspect='auto')
  932. # #ax.invert_yaxis()
  934. # xticks = np.array(1 * 2**np.arange(0, 9, 2))
  935. # f = interp1d(ntimes, np.arange(len(ntimes)))
  936. # ax.set_xticks(f(xticks * 60))
  937. # yticks = np.array(200 * 2**np.arange(0, 6, 1))
  938. # yticks = np.hstack((yticks, 10000))
  939. # f = interp1d(nneurons, np.arange(len(nneurons)))
  940. # ax.set_yticks(f(yticks))
  941. # ax.tick_params(labelsize='small')
  942. # #ax.axis('square')
  943. # ax.set_title(['low', 'medium', 'high'][j%3], fontsize='small')
  944. # if j==1:
  945. # ax.text(0.5, 1.2, 'Gaussian noise + smoothing', fontsize='medium',
  946. # ha='center', transform=ax.transAxes, fontstyle='italic')
  947. # elif j==4:
  948. # ax.text(0.5, 1.2, 'Poisson noise', fontsize='medium',
  949. # ha='center', transform=ax.transAxes, fontstyle='italic')
  951. # if ni==0:
  952. # if j==1:
  953. # axin = ax.inset_axes([0.05, -0.25, 0.8, 0.1]);
  954. # cb = plt.colorbar(im, cax=axin, orientation='horizontal')
  955. # cb.ax.tick_params(labelsize='small')
  956. # axin.text(1.05, 0.5, 'power-law exponent ($\\alpha$)', fontsize='small',
  957. # ha='left', transform=axin.transAxes, va='center')
  958. # cb.ax.set_xticks([0.5, 1, 1.5, 2])
  959. # #cb.ax.yaxis.label.set_size('medium')
  960. # elif j==0:
  961. # ax.set_ylabel('# of neurons')
  962. # ax.set_xlabel('time (min)')
  964. # if j==0:
  965. # transl = mtransforms.ScaledTranslation(-50 / 72, -5 / 72, fig.dpi_scale_trans)
  966. # il = plot_label(ltr, il, ax, transl)
  967. # ax.set_yticklabels(yticks)
  968. # ax.set_xticklabels(xticks, rotation=30, ha='right', va='top')
  969. # else:
  970. # ax.set_yticklabels([])
  971. # ax.set_xticklabels([])
  973. # return fig

fig2.py at commit 2c15edf, under GPL-3.0 · at the source

Overview

Authors: Marius Pachitariu1, Lin Zhong1, Alexa Gracias1, Amanda Minisi1, Crystall Lopez1, Carsen Stringer1
ORCID iDs: Marius Pachitariu, Lin Zhong, Carsen Stringer
  1. HHMI Janelia Research Campus,Ashburn, VA USA
Institutions: Janelia Research Campus (United States)
Journal: Nature, volume 655, issue 8124, pages 990-996
Dates: received 29 January 2025; accepted 10 April 2026; published online 20 May 2026; in print 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41586-026-10528-1 · PMID 42162432 · PMCID PMC13391357 · OpenAlex W4406241453
Open access: hybrid, a free copy (OpenAlex)
Preprint: doi:10.1101/2025.01.10.632397
Status: code verified
Categories: computational modeling (no new data) (modality), mouse (organism)
Methods: Statistics, Smoothing, state filtering, decompositions, Machine learning, Preprocessing, Connectivity, fMRI & imaging, Single-unit activity, calcium imaging, Spectral & time-frequency
Keywords: Network models, Cellular neuroscience
MeSH: Models, Neurological*, Nerve Net*, Neurons*, Animals, CA1 Region, Hippocampal, Memory, Short-Term, Mice, Time Factors (* major topic)
Topic: Neural Networks and Applications (Artificial Intelligence, Computer Science), according to OpenAlex
Citations: cited by 5 papers (Europe PMC); 110 references in the paper

Abstract

Intrinsically generated, brainwide neural activity displays macroscopic coordination among large populations of neurons that persists beyond the biophysical timescales of individual neurons1–3. It is not well understood how these macroscopic behaviours arise from microscopic, short-lived interactions between pairs of neurons. Here we show that the eigenvalue spectrum and dynamical properties of large-scale neural recordings in mice are similar to those produced by linear dynamics governed by a random symmetric matrix that is critically normalized. An exception was population activity in hippocampal area CA1, which resembled an efficient, uncorrelated neural code that may be optimized for information storage capacity. High-dimensional, global activity modes emerged in critically normalized artificial networks and persisted under sparse, clustered or spatial connectivity. These dynamics were useful for solving time-dependent tasks such as a zero-shot working memory task.

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

mouseland/critical_init

License: GPL-3.0
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: 2c15edf4e770165fc3962dcc3920c8bcaf555bed, 4 June 2026
Languages: Python (13), Jupyter (7)
Size: 32 files, 20 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, environment (requirements.txt), 7 notebooks
Not found: CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (17 files), PyTorch (13 files), Matplotlib (10 files), SciPy (8 files), scikit-learn (2 files), Neurodata Without Borders (PyNWB, MatNWB) (1 file), Suite2p (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
22 files
At the source: github.com/mouseland/critical_init

Zenodo 19322086

License: GPL-3.0
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: the references
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (17 files), PyTorch (13 files), Matplotlib (10 files), SciPy (8 files), scikit-learn (2 files), Neurodata Without Borders (PyNWB, MatNWB) (1 file), Suite2p (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers (HTTP 200)
  • 28 September 2026: the link answers (HTTP 200)
22 files
At the source:

Code availability

Code to reproduce all the analyses and figures is available at GitHub (https://github.com/mouseland/critical_init)110.

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:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 40 scripts, each with its path and the digest of its content;
  • 13 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 new neural recordings from this study are available at Figshare (10.25378/janelia.27854448)108. Datasets in Fig. 4 are publicly available at Figshare (10.25378/janelia.23712957)75. The eight-probe Neuropixels recordings were published previously109; we have uploaded the Kilosort4 processing of these recordings to Figshare (10.25378/janelia.27854448)108. The datasets in Extended Data Fig. 8 are publicly available88,90,92,95,97,99.

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, 6 authors, 2 keywords, 8 MeSH terms, 84 references.

Cite

This paper

Pachitariu, M., Zhong, L., Gracias, A., Minisi, A., Lopez, C., & Stringer, C. (2026). A critical initialization for biological neural networks. Nature, 655(8124), 990-996. https://doi.org/10.1038/s41586-026-10528-1

BibTeX

@article{pachitariu2026critical,
author = {Pachitariu, Marius and Zhong, Lin and Gracias, Alexa and Minisi, Amanda and Lopez, Crystall and Stringer, Carsen},
title = {{A critical initialization for biological neural networks}},
journal = {Nature},
year = {2026},
month = may,
volume = {655},
number = {8124},
pages = {990--996},
publisher = {Nature Portfolio},
issn = {0028-0836},
doi = {10.1038/s41586-026-10528-1},
url = {https://doi.org/10.1038/s41586-026-10528-1},
pmid = {42162432},
pmcid = {PMC13391357}
}

RIS

TY - JOUR
AU - Pachitariu, Marius
AU - Zhong, Lin
AU - Gracias, Alexa
AU - Minisi, Amanda
AU - Lopez, Crystall
AU - Stringer, Carsen
TI - A critical initialization for biological neural networks
T2 - Nature
J2 - Nature
PY - 2026
DA - 2026/05/20
VL - 655
IS - 8124
SP - 990
EP - 996
SN - 0028-0836
PB - Nature Portfolio
DO - 10.1038/s41586-026-10528-1
UR - https://doi.org/10.1038/s41586-026-10528-1
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41586-026-10528-1",
"type": "article-journal",
"title": "A critical initialization for biological neural networks",
"container-title": "Nature",
"author": [
{
"family": "Pachitariu",
"given": "Marius"
},
{
"family": "Zhong",
"given": "Lin"
},
{
"family": "Gracias",
"given": "Alexa"
},
{
"family": "Minisi",
"given": "Amanda"
},
{
"family": "Lopez",
"given": "Crystall"
},
{
"family": "Stringer",
"given": "Carsen"
}
],
"container-title-short": "Nature",
"volume": "655",
"issue": "8124",
"page": "990-996",
"DOI": "10.1038/s41586-026-10528-1",
"PMID": "42162432",
"PMCID": "PMC13391357",
"ISSN": "0028-0836",
"publisher": "Nature Portfolio",
"URL": "https://doi.org/10.1038/s41586-026-10528-1",
"language": "en",
"issued": {
"date-parts": [
[
2026,
5,
20
]
]
}
}

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

Similar papers

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

[1] doi:10.1126/sciadv.aed6417 [code]
Intrinsic timing, not temporal prediction, underlies ramping dynamics in visual and parietal cortex during passive behavior.
Journal: Science advances
In common: scikit-learn, SciPy, Matplotlib, 1 other tool, mouse, 8 references
[2] 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: Neurodata Without Borders (PyNWB, MatNWB), PyTorch, scikit-learn, 3 other tools, mouse, 5 references
[3] doi:10.1016/j.crmeth.2026.101421 [code]
EthoPy provides an accessible platform for reproducible behavioral neuroscience.
Journal: Cell reports methods
In common: Neurodata Without Borders (PyNWB, MatNWB), scikit-learn, SciPy, 2 other tools, mouse, 6 references
[4] doi:10.1016/j.neuron.2026.07.016 [code]
Inferring brain-wide interactions using data-constrained recurrent neural network models.
Journal: Neuron
In common: Matplotlib, NumPy, computational modeling (no new data), mouse, 8 references
[5] doi:10.1038/s41467-026-71458-0 [code]
Early differential impact of MeCP2 mutations on functional networks in Rett syndrome patient-derived human cortical organoids.
Journal: Nature communications
In common: Suite2p, Neurodata Without Borders (PyNWB, MatNWB), PyTorch, 4 other tools, 1 reference
[6] doi:10.1038/s41467-026-71725-0 [code]
Interactions across hemispheres in prefrontal cortex reflect global cognitive processing.
Journal: Nature communications
In common: scikit-learn, SciPy, Matplotlib, 1 other tool, 6 references
[7] doi:10.3390/biomimetics11080569 [code]
Pretraining of Embodied Recurrent Networks Bridges the Gap Between Artificial and Cortical Neural Activities.
Journal: Biomimetics (Basel, Switzerland)
In common: PyTorch, scikit-learn, SciPy, 2 other tools, 4 references
[8] doi:10.1093/nsr/nwag420 [code]
Constructing mesoscale functionomics by neural dynamics subspace clustering.
Journal: National science review
In common: 8 references
[9] doi:10.1038/s41467-026-74227-1 [code]
Age-related changes in behavioural and neural variability in a decision-making task.
Journal: Nature communications
In common: scikit-learn, SciPy, Matplotlib, 1 other tool, mouse, 5 references
[10] doi:10.1038/s41592-026-03076-z [code]
Neuropixels Opto: combining high-resolution electrophysiology and optogenetics.
Journal: Nature methods
In common: Neurodata Without Borders (PyNWB, MatNWB), PyTorch, SciPy, 2 other tools, mouse, 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.