OSCR

Redundant prefrontal hemispheres adapt storage strategy to working memory demands.

Code ↔ Paper

16 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 16 matches
  1. [1] § Methods › Neural data analysis › Gaussian mixture model ↔ Figure5/Figure5.ipynb, lines 1508–1543 · score 0.74 · Gaussian mixture, dual reactivations, reactivation strengths, covariance, classes, GMM
  2. [2] § Methods › Neural data analysis › Higher-capacity model ↔ Figure7/MultiitemModel_HigherCapacity.py, lines 224–273 · score 0.71 · weight profile, stimulus width, 6.1 deg, 14.4 deg, capacity, excitation
  3. [3] § Methods › Neural data analysis › Higher-capacity model ↔ Figure7/TwoAreaModelVaryingConnect_MultiitemParams.py, lines 216–265 · score 0.67 · weight profile, stimulus width, 6.1 deg, 14.4 deg, excitation, 50 deg
  4. [4] § Methods › Neural data analysis › Decoder-behavior error correlations ↔ Figure3/Figure3.ipynb, lines 1157–1222 · score 0.67 · circ_corrcc, circularly correlated, response errors, ipsilateral, contralateral, behavioral
  5. [5] § Methods › Behavioral data analysis › Serial dependence strength fit ↔ Figure5/Figure5.ipynb, lines 804–896 · score 0.66 · fit behavioral, decoder errors, dog1, serial dependence, OLS, neural
  6. [6] § Results › Each prefrontal hemisphere had private history effects, suggesting weak interhemispheric connectivity ↔ Figure5/Figure5.ipynb, lines 804–896 · score 0.65 · model fit, history drift, decoder error, dog1, serial dependence, OLS
  7. [7] § Methods › Neural data analysis › Single-trial decoders ↔ Figure3/Figure3_Figure5_SingleTrialDecoder.py, lines 1–24 · score 0.60 · single trial decoders, cross validation, Delay Decoder, reactivations, prediction
  8. [8] § Results › Both hemispheres equally reflected trial-by-trial behavioral variability with weakly correlated memory representations ↔ Figure3/Figure3_CorrelateEyeResponse.py, lines 103–191 · score 0.59 · eye movements, gaze errors, error correlations, monkeys, connections
  9. [9] § Methods › Behavioral data analysis › Determining the most-interfering non-target ↔ Supplement/Supplement_DoGSerialBias.py, lines 55–142 · score 0.55 · relative location, dog1, optimal, hyperparameter, BIC, sigma
  10. [10] § Methods › Neural data analysis › Cross-temporal decoders ↔ Figure2/Figure2_DecoderMatrix.py, lines 28–177 · score 0.55 · decoder matrices, linear regression, trained, model
  11. [11] § Methods › Neural data analysis › Uninstructed eye-movement controls ↔ Figure3/Figure3_CorrelateEyeResponse.py, lines 103–191 · score 0.54 · eye movements, error correlations, gaze
  12. [12] § Methods › Neural data analysis › Population decoding of stimulus location ↔ Figure2/Figure2_DecoderMatrix.py, lines 28–177 · score 0.54 · trained model, randomly shuffled, MSE, score, decoding
  13. [13] § Results › Each prefrontal hemisphere had private history effects, suggesting weak interhemispheric connectivity ↔ Figure5/Figure5.ipynb, lines 1342–1380 · score 0.53 · Gaussian Mixture, dual reactivations, class, GMM, diagonal, strength
  14. [14] § Methods › Neural data analysis › Bump-attractor model ↔ Figure7/TwoAreaModelVaryingConnect_MultiitemParams.py, lines 267–306 · score 0.52 · excitatory neurons, GABA, NMDA, AMPA, profile, std
  15. [15] § Methods › Behavioral data analysis › Serial dependence strength fit ↔ Supplement/Supplement_DoGSerialBias.py, lines 55–142 · score 0.52 · OLS model, dog1, derivative, serial, Gaussian, fit
  16. [16] § Results › Each hemisphere stably represented the full visual space with a contralateral preference ↔ Figure2/Figure2.ipynb, lines 549–638 · score 0.51 · contra ipsi, shorter stimulus, scored, shuffle, Pe, trained

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 · 1,797 lines · 82 KB · no license · 4 matches

  1. # %%
  2. %cd ..
  4. # %%
  5. import matplotlib.pyplot as plt
  6. plt.style.use("stylefile.mplstyle")
  9. # Define colors
  10. colors = {
  11. "Single": "#333333",
  12. "Left": "#246EB9",
  13. "Right": "#8B1E3F",
  14. "Ipsilateral": "#1B9E77",
  15. "Contralateral": "#D95F02",
  16. "Within": "#2B4162",
  17. "Across": "#E89D0B",
  18. 'SerialBiasWeak': '#254441',
  19. 'SerialBias': '#43AA8B',\
  20. 'ReactivationWeak': '#7A8A99',
  21. 'Reactivation': '#12719E'
  22. }
  24. from matplotlib.colors import LinearSegmentedColormap
  25. left_right_cmap = LinearSegmentedColormap.from_list("LeftRight",\
  26. [colors['Left'], "white", colors['Right']])
  27. ipsi_contra_cmap = LinearSegmentedColormap.from_list("IpsiContra",\
  28. [colors['Ipsilateral'], "white", colors['Contralateral']])
  29. within_cmap = LinearSegmentedColormap.from_list("Within",\
  30. ["white", colors['Within']])
  31. across_cmap = LinearSegmentedColormap.from_list("Across",\
  32. ["white", colors['Across']])
  34. # %%
  35. ## circstats
  36. def len2(x):
  37. if type(x) is not type([]):
  38. if type(x) is not type(array([])):
  39. return -1
  40. return len(x)
  42. def phase2(x):
  43. if not np.isnan(x):
  44. return phase(x)
  45. return nan
  47. def circdist(angles1,angles2):
  48. ''' calculates circular distance of angles [rad]'''
  49. if len2(angles2) < 0:
  50. if len2(angles1) > 0:
  51. angles2 = [angles2]*len(angles1)
  52. else:
  53. angles2 = [angles2]
  54. angles1 = [angles1]
  55. if len2(angles1) < 0:
  56. angles1 = [angles1]*len(angles2)
  57. return array(list(map(lambda a1,a2: phase2(np.exp(1j*a1)/np.exp(1j*a2)), angles1,angles2)))
  59. # %%
  60. def sig_bar(sigs,axis,y,ax,color):
  61. w=np.diff(axis)[0]
  62. for s in sigs:
  63. beg =axis[s]-w/2
  64. end = axis[s]+w/2
  65. ax.fill_between([beg,end],[y[0],y[0]],[y[1],y[1]],color=color)
  67. # %%
  68. def sign_rl(rel_loc):
  69. # returns sign of array unless value=0, returns 1
  70. sign = [np.sign(rel_loc[rl]) if rel_loc[rl]!=0.0 else 1 for rl in range(len(rel_loc))]
  71. return np.array(sign)
  73. # %%
  74. import numpy as np
  75. def plot_twolines_full(R=[],base=[], bins=0.2,labelR='',labelB='', errorbars='SEM', MeanType = np.nanmean,\
  76. borders=[], ylabel='decoding', yticks = False, titel='', significances=False,\
  77. currentTrial=False, end_border = 15, shorten_delay = 0,\
  78. refline =0, colors=[colors['Within'], colors['Across']]):
  80. offset_start = borders[1]-borders[0]
  81. x = (np.linspace(borders[1], borders[5]-shorten_delay, borders[5]-borders[1]-shorten_delay)-\
  82. (borders[2]-borders[1])-offset_start)*bins
  83. x3 = (np.linspace(borders[5]+shorten_delay, borders[8], borders[8]-borders[5]+1-shorten_delay)-\
  84. (borders[7]-borders[1])-offset_start)*bins
  85. if end_border < 8:# need end of one trial, start of next trial
  86. x3 = (np.linspace(borders[5]+shorten_delay, borders[end_border], borders[end_border]-borders[5]-shorten_delay)-\
  87. (borders[7]-borders[1])-offset_start)*bins
  89. width_radius = [len(x), len(x3)]
  90. # only add 3rd axis if we look at start of next trial
  91. if end_border > 12:
  92. end_prev = (np.linspace(borders[9], borders[10], borders[10]-borders[9]+1)-\
  93. (borders[13]-borders[1])-offset_start)*bins
  94. start_curr = (np.linspace(borders[12], borders[end_border], borders[end_border]-borders[12])-\
  95. (borders[13]-borders[1])-offset_start)*bins
  96. x5 = np.append(end_prev,start_curr)
  97. width_radius = [len(x), len(x3), len(x5)]# len(x3),
  98. elif end_border >=10:
  99. end_prev = (np.linspace(borders[9], borders[10], borders[10]-borders[9]+1)-offset_start)*bins
  100. x5 = end_prev
  101. width_radius = [len(x), len(x3), len(x5)]# len(x3),
  103. compare = 0.05#/len(np.append(np.append(x,x2),x5))
  104. # TTESTS:
  105. pvals = ttest_rel(R,base, axis=0, nan_policy='omit')[1]
  106. pvals1 = ttest_1samp(R, 0, axis=0, nan_policy='omit')[1]
  107. pvals2 = ttest_1samp(base, 0, axis=0, nan_policy='omit')[1]
  109. # compute mean and errorbars
  110. mean_R = MeanType(R, axis=0)
  111. mean_base = MeanType(base, axis=0)
  112. if MeanType == circmean:
  113. mean_R = MeanType(R, axis=0, low=-np.pi, high=np.pi)
  114. mean_base = MeanType(base, axis=0, low=-np.pi, high=np.pi)
  115. if errorbars=='CI':
  116. std_R = 2*sem(R, axis=0, nan_policy='omit')
  117. std_base = 2*sem(base, axis=0, nan_policy='omit')
  118. else:
  119. std_R = sem(R, axis=0, nan_policy='omit')
  120. std_base = sem(base, axis=0, nan_policy='omit')
  122. #f, (ax1,ax5) = plt.subplots(1, 2,sharey=True,figsize=(3.8,1.6), gridspec_kw={'width_ratios': width_radius})
  123. #f, (ax1,ax3,ax5) = plt.subplots(1, 3,sharey=True,figsize=(16,5), gridspec_kw={'width_ratios': width_radius})
  124. if end_border < 10:
  125. f, (ax1,ax3) = plt.subplots(1, 2,sharey=True,figsize=(3.,2), gridspec_kw={'width_ratios': width_radius})
  126. else:
  127. f, (ax1,ax3,ax5) = plt.subplots(1, 3,sharey=True,figsize=(3.5, 2), gridspec_kw={'width_ratios': width_radius})
  129. ###################### PLOT FIXTION TO DELAY PREVIOUS
  130. cut = range(borders[1], borders[5]-shorten_delay)
  131. # within
  132. ax1.plot(x,mean_R[cut], color=colors[0], label=labelR)#+', '+errorbars
  133. ax1.fill_between(x, mean_R[cut]-np.array(std_R[cut]), mean_R[cut]+np.array(std_R[cut]), color=colors[0], alpha=0.3)
  134. # across
  135. ax1.plot(x,mean_base[cut], color=colors[1], label=labelB)#+', '+errorbars
  136. ax1.fill_between(x, mean_base[cut]-np.array(std_base[cut]), mean_base[cut]+np.array(std_base[cut]),\
  137. color=colors[1], alpha=0.3)
  138. if currentTrial:
  139. ax1.set_xlabel('fix$_{n+1}$', labelpad=0)
  140. else:
  141. ax1.set_xlabel('fix$_{n}$', labelpad=0)
  142. ax1.spines['right'].set_visible(False)
  143. ax1.spines['top'].set_visible(False)
  144. ax1.xaxis.set_ticks_position('bottom')
  145. ax1.yaxis.set_ticks_position('left')
  146. ax1.axhline(refline, color='#333333', linestyle='--', alpha=0.7)
  147. ax1.set_ylabel(ylabel)
  148. if yticks != False:
  149. ax1.set_yticks(yticks)
  150. if (len(np.unique(base)) != 1):
  151. ax1.legend()
  153. ###################### PLOT DELAY TO SACCADE PREVIOUS
  154. cut = range(borders[5]+shorten_delay, borders[8]+1)
  155. if end_border < 8:
  156. cut = range(borders[5]+shorten_delay, borders[end_border])
  157. # within
  158. ax3.plot(x3,mean_R[cut], color=colors[0], label=labelR)
  159. ax3.fill_between(x3, mean_R[cut]-np.array(std_R[cut]), mean_R[cut]+np.array(std_R[cut]), color=colors[0], alpha=0.3)
  160. # across
  161. ax3.plot(x3,mean_base[cut], color=colors[1], label=labelB)
  162. ax3.fill_between(x3, mean_base[cut]-np.array(std_base[cut]), mean_base[cut]+np.array(std_base[cut]), color=colors[1], alpha=0.3)
  163. if currentTrial:
  164. ax3.set_xlabel('response$_{n+1}$', labelpad=0)
  165. else:
  166. ax3.set_xlabel('response$_{n}$', labelpad=0)
  167. ax3.spines['right'].set_visible(False)
  168. ax3.spines['left'].set_visible(False)
  169. ax3.spines['top'].set_visible(False)
  170. ax3.xaxis.set_ticks_position('bottom')
  171. ax3.axhline(refline, color='#333333', linestyle='--', alpha=0.7)
  172. #ax3.set_ylabel(ylabel)
  173. #ax3.set_yticks([0,0.1])
  175. ###################### PLOT ITI PREVIOUS TO DELAY CURRENT TRIAL
  176. if end_border >= 10:
  177. pc_base = mean_base[borders[9]:borders[10]+1]
  178. pcstd_base = std_base[borders[9]:borders[10]+1]
  179. pc_R = mean_R[borders[9]:borders[10]+1]
  180. pcstd_R = std_R[borders[9]:borders[10]+1]
  181. if end_border > 12:
  182. pc_base = np.append(mean_base[borders[9]:borders[10]+1],mean_base[borders[12]:borders[15]])
  183. pcstd_base = np.append(std_base[borders[9]:borders[10]+1],std_base[borders[12]:borders[15]])
  184. pc_R = np.append(mean_R[borders[9]:borders[10]+1],mean_R[borders[12]:borders[15]])
  185. pcstd_R = np.append(std_R[borders[9]:borders[10]+1],std_R[borders[12]:borders[15]])
  187. # plot reactivaton period
  188. # within
  189. ax5.plot(x5,pc_R, color=colors[0])
  190. ax5.fill_between(x5, pc_R-np.array(pcstd_R), pc_R+np.array(pcstd_R), color=colors[0], alpha=0.3)
  191. # across
  192. ax5.plot(x5,pc_base, color=colors[1])
  193. ax5.fill_between(x5, pc_base-np.array(pcstd_base), pc_base+np.array(pcstd_base), color=colors[1], alpha=0.3)
  194. # plot baseline
  195. ax5.set_xlabel('fix$_{n+1}$', labelpad=0)
  196. #ax5.set_xticks([])
  197. ax5.spines['right'].set_visible(False)
  198. ax5.spines['left'].set_visible(False)
  199. ax5.spines['top'].set_visible(False)
  200. ax5.xaxis.set_ticks_position('bottom')
  201. ax5.axhline(refline, color='#333333', linestyle='--', alpha=0.7)
  202. y0=ax5.get_ylim()[0]
  203. y1=ax5.get_ylim()[1]
  204. ax5.fill_between([(borders[14]-borders[12])*bins, (borders[15]-borders[12])*bins],\
  205. y0, y1, color='grey', alpha=0.2)
  207. ###################### MARK IMPORTANT TIME PERIODS
  208. f.text(0.55, 0.02, 'time (s) from', ha='center', fontsize=8)
  209. # # # plot significances
  210. y0=ax3.get_ylim()[0]
  211. y1=ax3.get_ylim()[1]
  212. off = (y1-y0)/10
  213. marker=250
  214. ax1.fill_between([(borders[3]-borders[2])*bins, (borders[4]-borders[2])*bins], y0, y1, color='grey', alpha=0.2)
  216. w = (y1-y0)/50
  217. off = [y1-w, y1]
  218. sigs = np.where(pvals[borders[1]:borders[5]-shorten_delay]<compare)[0]
  219. sig_bar(sigs,x,off,ax1,'#333333')
  220. sigs = np.where(pvals[borders[5]+shorten_delay:borders[8]+1]<compare)[0]
  221. sig_bar(sigs,x3,off,ax3,'#333333')
  222. if end_border >= 10:
  223. sigs = np.where(pvals[borders[9]:borders[10]+1]<compare)[0]
  224. sig_bar(sigs,x5,off,ax5,'#333333')
  225. f.text(0.23, 0.87, 'Stim$_{n}$', ha='center', fontsize=8)
  226. f.text(0.46, 0.87, 'Delay$_{n}$', ha='center', fontsize=8)
  227. f.text(0.85, 0.87, 'ITI$_{n}$', ha='center', fontsize=8)
  228. elif end_border > 12:
  229. sigs = np.where(np.append(pvals[borders[9]:borders[10]+1],pvals[borders[12]:borders[15]])<compare)[0]
  230. sig_bar(sigs,x5,off,ax5,'#333333')
  231. f.text(0.94, 0.87, 'Stim$_{n+1}$', ha='center', fontsize=8)
  232. else:
  233. f.text(0.29, 0.87, 'Stim$_{n}$', ha='center', fontsize=8)
  234. f.text(0.63, 0.87, 'Delay$_{n}$', ha='center', fontsize=8)
  236. if significances:
  237. # within area
  238. off = [y1-2.5*w, y1-1.5*w]
  239. sigs = np.where(pvals1[borders[1]:borders[5]-shorten_delay]<compare)[0]
  240. sig_bar(sigs,x,off,ax1,colors[0])
  241. sigs = np.where(pvals1[borders[5]+shorten_delay:borders[8]+1]<compare)[0]
  242. sig_bar(sigs,x3,off,ax3,colors[0])
  243. if end_border >= 10:
  244. sigs = np.where(pvals1[borders[9]:borders[10]+1]<compare)[0]
  245. sig_bar(sigs,x5,off,ax5,colors[0])
  246. elif end_border > 12:
  247. sigs = np.where(np.append(pvals1[borders[9]:borders[10]+1],pvals1[borders[12]:borders[15]])<compare)[0]
  248. sig_bar(sigs,x5,off,ax5,colors[0])
  250. # across areas
  251. off = [y1-4*w, y1-3*w]
  252. sigs = np.where(pvals2[borders[1]:borders[5]-shorten_delay]<compare)[0]
  253. sig_bar(sigs,x,off,ax1,colors[1])
  254. sigs = np.where(pvals2[borders[5]+shorten_delay:borders[8]+1]<compare)[0]
  255. sig_bar(sigs,x3,off,ax3,colors[1])
  256. if end_border >= 10:
  257. sigs = np.where(pvals2[borders[9]:borders[10]+1]<compare)[0]
  258. sig_bar(sigs,x5,off,ax5,colors[1])
  259. elif end_border > 12:
  260. sigs = np.where(np.append(pvals2[borders[9]:borders[10]+1],pvals2[borders[12]:borders[15]])<compare)[0]
  261. sig_bar(sigs,x5,off,ax5,colors[1])
  264. plt.suptitle(titel)
  266. return
  268. # %%
  269. def plot_threelines_full(R=[],base=[],base2=[], bins=200,labelR='',labelB='',labelB2='',\
  270. errorbars='SEM', MeanType = np.nanmean, shorten_delay = 0,\
  271. borders=[], ylabel='decoding', cutyaxis = False, titel='',\
  272. colors=['#333333', colors['Within'], colors['Across']]):
  274. offset_start = borders[1]-borders[0]
  275. x = (np.linspace(borders[1], borders[5]-shorten_delay, borders[5]-borders[1]+1-shorten_delay)-(borders[2]-borders[1])-offset_start)*bins
  276. x3 = (np.linspace(borders[5]+shorten_delay, borders[8], borders[8]-borders[5]+1-shorten_delay)-(borders[7]-1-borders[1])-offset_start)*bins
  277. # need end of one trial, start of next trial
  278. end_prev = (np.linspace(borders[9], borders[10], borders[10]-borders[9]+1)-(borders[13]-borders[1])-offset_start)*bins
  279. start_curr = (np.linspace(borders[12], borders[15], borders[15]-borders[12])-(borders[13]-borders[1])-offset_start)*bins
  280. x5 = np.append(end_prev,start_curr)
  281. len(x5)
  282. width_radius = [len(x), len(x3), len(x5)]# len(x3),
  284. compare = 0.05#/len(np.append(np.append(x,x2),x5))
  285. # TTESTS:
  286. pvals = ttest_rel(R,base, axis=0, nan_policy='omit')[1]
  287. pvals1 = ttest_1samp(R, 0, axis=0, nan_policy='omit')[1]
  288. pvals2 = ttest_1samp(base, 0, axis=0, nan_policy='omit')[1]
  290. mean_R = MeanType(R, axis=0)
  291. mean_base = MeanType(base, axis=0)
  292. mean_base2 = MeanType(base2, axis=0)
  293. if MeanType == circmean:
  294. mean_R = MeanType(R, axis=0, low=-np.pi, high=np.pi)
  295. mean_base = MeanType(base, axis=0, low=-np.pi, high=np.pi)
  296. mean_base2 = MeanType(base2, axis=0, low=-np.pi, high=np.pi)
  297. if errorbars=='CI':
  298. std_R = 2*sem(R, axis=0, nan_policy='omit')
  299. std_base = 2*sem(base, axis=0, nan_policy='omit')
  300. std_base2 = 2*sem(base2, axis=0, nan_policy='omit')
  301. else:
  302. std_R = sem(R, axis=0, nan_policy='omit')
  303. std_base = sem(base, axis=0, nan_policy='omit')
  304. std_base2 = sem(base2, axis=0, nan_policy='omit')
  307. f, (ax1,ax3,ax5) = plt.subplots(1, 3,sharey=True,figsize=(4.8, 2.1), gridspec_kw={'width_ratios': width_radius})
  309. ###################### PLOT FIXTION TO DELAY PREVIOUS
  310. cut = range(borders[1], borders[5]+1-shorten_delay)
  311. # within
  312. ax1.plot(x,mean_base[cut], color=colors[1], label=labelB)#+', '+errorbars
  313. ax1.fill_between(x, mean_base[cut]-np.array(std_base[cut]), mean_base[cut]+np.array(std_base[cut]), color=colors[1], alpha=0.2)
  314. # base2
  315. ax1.plot(x,mean_base2[cut], color=colors[2], label=labelB2)#+', '+errorbars
  316. ax1.fill_between(x, mean_base2[cut]-np.array(std_base2[cut]), mean_base2[cut]+np.array(std_base2[cut]), color=colors[2], alpha=0.2)
  317. # across
  318. ax1.plot(x,mean_R[cut], color=colors[0], label=labelR)#+', '+errorbars
  319. ax1.fill_between(x, mean_R[cut]-np.array(std_R[cut]), mean_R[cut]+np.array(std_R[cut]), color=colors[0], alpha=0.2)
  320. ax1.set_xlabel('fix$_{n}$', labelpad=0)
  321. ax1.spines['right'].set_visible(False)
  322. ax1.spines['top'].set_visible(False)
  323. ax1.xaxis.set_ticks_position('bottom')
  324. ax1.yaxis.set_ticks_position('left')
  325. ax1.axhline(0, color='#333333', linestyle='--', alpha=0.7)
  326. ax1.set_ylabel(ylabel)
  327. ax1.legend()
  329. # ###################### PLOT DELAY TO SACCADE PREVIOUS
  330. cut = range(borders[5]+shorten_delay, borders[8]+1)
  331. # within
  332. ax3.plot(x3,mean_base[cut], color=colors[1], label=labelB)
  333. ax3.fill_between(x3, mean_base[cut]-np.array(std_base[cut]), mean_base[cut]+np.array(std_base[cut]), color=colors[1], alpha=0.2)
  334. # base2
  335. ax3.plot(x3,mean_base2[cut], color=colors[2], label=labelB2)
  336. ax3.fill_between(x3, mean_base2[cut]-np.array(std_base2[cut]), mean_base2[cut]+np.array(std_base2[cut]), color=colors[2], alpha=0.2)
  337. # across
  338. ax3.plot(x3,mean_R[cut], color=colors[0], label=labelR)
  339. ax3.fill_between(x3, mean_R[cut]-np.array(std_R[cut]), mean_R[cut]+np.array(std_R[cut]), color=colors[0], alpha=0.2)
  340. ax3.set_xlabel('response$_{n}$', labelpad=0)
  341. ax3.spines['right'].set_visible(False)
  342. ax3.spines['left'].set_visible(False)
  343. ax3.spines['top'].set_visible(False)
  344. ax3.xaxis.set_ticks_position('bottom')
  345. ax3.axhline(0, color='#333333', linestyle='--', alpha=0.7)
  347. ###################### PLOT ITI PREVIOUS TO DELAY CURRENT TRIAL
  348. pc_base = np.append(mean_base[borders[9]:borders[10]+1],mean_base[borders[12]:borders[15]])
  349. std_base = np.append(std_base[borders[9]:borders[10]+1],std_base[borders[12]:borders[15]])
  350. pc_base2 = np.append(mean_base2[borders[9]:borders[10]+1],mean_base2[borders[12]:borders[15]])
  351. std_base2 = np.append(std_base2[borders[9]:borders[10]+1],std_base2[borders[12]:borders[15]])
  352. pc_R = np.append(mean_R[borders[9]:borders[10]+1],mean_R[borders[12]:borders[15]])
  353. std_R = np.append(std_R[borders[9]:borders[10]+1],std_R[borders[12]:borders[15]])
  354. # plot
  355. # within
  356. ax5.plot(x5,pc_base, color=colors[1])
  357. ax5.fill_between(x5, pc_base-np.array(std_base), pc_base+np.array(std_base), color=colors[1], alpha=0.2)
  358. # base2
  359. ax5.plot(x5,pc_base2, color=colors[2])
  360. ax5.fill_between(x5, pc_base2-np.array(std_base2), pc_base2+np.array(std_base2), color=colors[2], alpha=0.2)
  361. # across
  362. ax5.plot(x5,pc_R, color=colors[0])
  363. ax5.fill_between(x5, pc_R-np.array(std_R), pc_R+np.array(std_R), color=colors[0], alpha=0.2)
  364. # plot baseline
  365. ax5.set_xlabel('fix$_{n+1}$', labelpad=0)
  366. #ax5.set_xticks([])
  367. ax5.spines['right'].set_visible(False)
  368. ax5.spines['left'].set_visible(False)
  369. ax5.spines['top'].set_visible(False)
  370. ax5.xaxis.set_ticks_position('bottom')
  371. ax5.axhline(0, color='#333333', linestyle='--', alpha=0.7)
  372. if cutyaxis != False:
  373. ax5.set_yticks(cutyaxis)
  374. ax5.set_ylim(cutyaxis)
  376. ###################### MARK IMPORTANT TIME PERIODS
  377. y0=ax5.get_ylim()[0]
  378. y1=ax5.get_ylim()[1]
  379. off = (y1-y0)/10
  380. #marker=30
  381. marker=250
  382. # ax1.plot(0, color='midnightblue', alpha=0.5)
  383. ax1.fill_between([(borders[3]-borders[2])*bins, (borders[4]-borders[2])*bins], y0, y1, color='grey', alpha=0.2)
  384. ax3.plot(0, color='#A0B2A6', alpha=0.5)
  385. ax5.plot(0, color='midnightblue', alpha=0.5)
  386. ax5.fill_between([(borders[14]-borders[12])*bins, (borders[15]-borders[12])*bins], y0, y1, color='grey', alpha=0.2)
  388. f.text(0.55, -0., 'time (ms) from', ha='center', fontsize=12)
  389. f.text(0.23, 0.75, 'Stim$_{n}$', ha='center', fontsize=10)
  390. f.text(0.46, 0.75, 'Delay$_{n}$', ha='center', fontsize=10)
  391. f.text(0.85, 0.75, 'ITI$_{n}$', ha='center', fontsize=10)
  392. f.text(0.94, 0.75, 'Stim$_{n+1}$', ha='center', fontsize=10)
  394. plt.suptitle(titel)
  396. return
  398. # %% [markdown]
  399. # # LOAD DATA
  401. # %%
  402. import pandas as pd
  403. import matplotlib.pyplot as plt
  404. import seaborn as sns
  405. import scipy
  406. from scipy.io import loadmat
  407. from scipy.stats import *
  408. from scipy.optimize import curve_fit
  409. from cmath import phase
  410. from numpy import array
  411. from scipy.sparse import csr_matrix
  412. import urllib
  413. import pickle
  414. from scipy.io import loadmat
  415. import glob
  416. import sklearn
  417. from sklearn.linear_model import LinearRegression
  418. from sklearn.linear_model import LogisticRegression
  419. from sklearn.model_selection import train_test_split
  420. #from pymicro.view.vol_utils import compute_affine_transform
  421. from sklearn.model_selection import LeaveOneOut
  422. from sklearn.metrics import accuracy_score
  423. from sklearn import preprocessing
  424. import statsmodels.formula.api as sf
  425. from sklearn import metrics
  426. from random import randint
  427. from numpy.linalg import inv
  428. import math
  429. import io
  430. #import h5py
  431. from circ_stats import *
  432. from patsy import dmatrices
  433. import statsmodels.api as sm
  434. import helpers as hf
  435. import statsmodels.formula.api as smf
  436. import copy
  437. import matplotlib.ticker as ticker
  438. from pingouin import circ_corrcc
  441. #monkeys = ["Sa", "Pe", "Wa"]
  442. #for m in monkeys:
  443. # files = np.sort(glob.glob('../Data/new/%s*.mat' %m))
  444. # for f in files:
  445. # print(f)
  448. with open('./Results_Full/df_serial.pickle', 'rb') as handle:
  449. #with open('./Results/df_serial_Sa0.pickle', 'rb') as handle:
  450. df_sb = pickle.load(handle)
  451. df_sb = df_sb.reset_index()
  452. df_sb.rel_loc = np.round(df_sb.rel_loc,3)
  454. monkeys=df_sb.monkey.unique()
  457. #df_behav = df_behav.loc[(df_behav.monkey=='Sa') | (df_behav.monkey=='Wa')].reset_index(drop=True)
  458. df_sb['session_continuous'] = [df_sb.monkey[i]+str(df_sb.session[i]) for i in df_sb.index]
  460. same_id = 1
  461. opp_id = 0
  462. border_id = 2
  464. # DEFINE DOG FIT PARAMETERS FOR EACH ANIMAL (based on BIC, questionable for Pe, Wa (delta BIC <2)
  465. sigma={'Sa':0.9, 'Pe':2.15, 'Wa':0.45}
  466. neural_sigma = {'Sa': 1.35, 'Pe': 0.9, 'Wa': 1.75}#
  467. reactivation_sigma={'Sa':2.95, 'Pe':1.15, 'Wa':1.55}
  469. # %%
  470. def cut_task_timings(timecourse=[], timeperiods=np.array([]), max_time=np.nan):
  471. """
  472. Cuts the variable timecourse together based on the varying task timing in multiple sessions.
  473. Aligns to shortest session timing
  475. timecourse : np.array() of shape (sessions x trials x time)
  476. timeperiods: np.array() of shape (sessions x time)
  478. cut_array: np.array() of shape (sessions x trials x minimum_time)
  479. """
  481. # find minimum taskperiod onsets / durations across sessions
  483. min_duration = np.min(np.diff(timeperiods),axis=0)
  484. min_onset = np.append(0,np.cumsum(min_duration))#add one more 0 at start so borders and diff add up
  487. # if no end time is given, do for all timesteps
  488. if np.isnan(max_time):
  489. max_time = timeperiods.shape[-1]
  491. # cut sessions together
  492. # cut_array = np.array([np.concatenate([timecourse[sess][timing:timing+min_duration[idx]]\
  493. # if idx not in [5, 9, 16, 20]\
  494. # else timecourse[sess][timing-min_duration[idx]:timing]
  495. # for idx,timing in enumerate(timeperiods[sess][:max_time])])\
  496. # for sess in range(len(timecourse))])
  497. cut_array = np.array([np.concatenate([timecourse[sess][timing:timing+min_duration[idx]]\
  498. if idx not in [5, 9, 16, 20]\
  499. else timecourse[sess][timeperiods[sess][idx+1]-min_duration[idx]:timeperiods[sess][idx+1]]
  500. for idx,timing in enumerate(timeperiods[sess][:max_time])])\
  501. for sess in range(len(timecourse))])
  504. #cut_array = np.array([cut_array])
  506. return cut_array, min_onset
  508. # %% [markdown]
  509. # ----
  511. # %% [markdown]
  512. # # Single trial analyses
  514. # %%
  515. def abs_err(prediction, all_targets):
  516. return np.abs([circdist(prediction, targ) for targ in all_targets])
  519. df_behav = df_sb.copy()
  520. #df_behav = df_behav.loc[(df_behav.monkey=='Sa') | (df_behav.monkey=='Wa')].reset_index(drop=True)
  521. df_behav['session_continuous'] = df_behav['monkey']+df_behav['session'].astype(str)
  522. # drop session with shorter stimulus
  523. df_behav.drop(index = np.where(df_behav.session_continuous == 'Pe2')[0], inplace=True)
  524. df_behav.reset_index(drop=True, inplace=True)
  528. # LOAD NEURAL DECODER DATA
  529. # INSTEAD OF RESPONSE ORTHO USE DELAY DECODER
  530. with open('../../Desktop/PhD/2_Smith/Results/SingleTrialDecoding/SingletrialDecoder.pickle', 'rb') as handle:
  531. #with open('./Results/Figure3/SingletrialDecoder_Sa0.pickle', 'rb') as handle:
  532. df = pickle.load(handle)
  534. df['session_continuous'] = df['monkey']+df['session'].astype(str)
  536. bins=200
  538. # compute errors and shuffles
  539. for neuron_type in ['combined', 'left', 'right']:
  540. print('Computing neuron_type: '+neuron_type)
  541. pred_prev=[]
  542. pred_prevortho=[]
  543. pred_curr=[]
  544. basecorr_prev=[]
  545. basecorr_prevortho=[]
  546. basecorr_curr=[]
  547. for i in df.index:
  548. targ_prev = np.round(np.angle(df.loc[i, 'targ_prev_xy']),3)
  549. targ_curr = np.round(np.angle(df.loc[i, 'targ_curr_xy']), 3)
  550. pprev = np.angle(df.loc[i, 'pred_complex_prev_'+neuron_type])
  551. # TODO! Change for all that re not delay decoders
  552. #pprevortho = np.angle(df.loc[i,'pred_complex_prev_ortho_'+neuron_type])
  553. pprevortho = np.angle(df.loc[i,'pred_complex_delay_'+neuron_type])
  554. pcurr = np.angle(df.loc[i,'pred_complex_curr_'+neuron_type])
  556. # compute error of prediction to target
  557. err_prev = circdist(pprev,targ_prev)
  558. err_prevortho = circdist(pprevortho,targ_prev)
  559. err_curr = circdist(pcurr,targ_curr)
  560. pred_prev.append(err_prev)
  561. pred_prevortho.append(err_prevortho)
  562. pred_curr.append(err_curr)
  564. df['prederror_prev_'+neuron_type] = pred_prev
  565. df['prederror_ortho_'+neuron_type] = pred_prevortho
  566. df['prederror_curr_'+neuron_type] = pred_curr
  568. # define start/end decoer errors for bumpdrift
  569. for neuron_type in ['combined', 'left', 'right']:
  570. df[neuron_type+'DecoderErr'] = [df.loc[i,'prederror_prev_'+neuron_type]\
  571. for i in df.index]
  574. df['leftHemi'] = df['hemifield_prev_left']
  575. df['rightHemi'] = df['hemifield_prev_right']
  577. IPSI = 1
  578. BORDER = 0
  579. CONTRA = -1
  581. # drop session with shorter stimulus
  582. df.drop(index = np.where(df.session_continuous == 'Pe2')[0], inplace=True)
  583. df.reset_index(drop=True, inplace=True)
  585. assert list(df.index) == list(df_behav.index), 'For merging dataframes: must be of same length'
  587. df['behav_response_prev'] = df_behav['response_prev']
  588. df['behav_response_curr'] = df_behav['response_curr']
  589. df['response_prev_curr'] = np.round(circdist(df.behav_response_prev, np.angle(df.targ_curr_xy)),3)
  590. df['delay_curr'] = df_behav.delay_curr
  593. # drop border trials (for ipsi contra analysis)
  594. df_pred_noBorder = df.drop(index=np.where(df.leftHemi==BORDER)[0])
  596. df.head()
  598. # %% [markdown]
  599. # ---
  601. # %% [markdown]
  602. # # Serial dependence
  604. # %% [markdown]
  605. # ### Fig 5a: Neural SD at stim., delay
  607. # %%
  608. # decoder error of delay_n+1 to target_n+1
  609. folded=False
  611. for m,monkey in enumerate(monkeys):#enumerate(['Sa', 'Wa']):#
  612. borders_mono=[]
  613. df_mono = df.loc[(df.monkey==monkey)].copy().reset_index(drop=True)
  614. df_mono['sign_prevcurr'] = sign_rl(df_mono.prev_curr.values)
  615. if folded==True:
  616. df_mono['prev_curr'] = np.abs(df_mono.prev_curr)
  617. prevcurr = np.unique(df.loc[df.session_continuous=='Sa0'].prev_curr)#np.unique(df_mono.prev_curr)#np.unique(df_mono.loc[df_mono.session==0].prev_curr)#np.unique(df_mono.prev_curr) max_num_delaysteps = np.min((np.array([df_mono.borders_full.values[i][16]-df_mono.borders_full.values[i][15]\
  619. # shape: (sessions x delay groups x prev-curr differences)
  620. decodersd = np.empty((len(np.unique(df_mono.session)), 2, len(prevcurr)))*np.nan
  621. mono_prederr = {'early':np.empty((len(df_mono)))*np.nan, 'late':np.empty((len(df_mono)))*np.nan}
  622. # compute SD for sessions separately
  623. for session_id,session in enumerate(np.unique(df_mono.session)):#
  624. df_sess = df_mono.loc[(df_mono.session==session)].copy()
  625. borders_full = df_sess.borders_full.values[0]
  626. d_start = borders_full[14]+1
  627. d_end = borders_full[18]-2
  629. # start of delay
  630. prederr = [circdist(np.angle(df_sess.pred_complex_curr_combined[i][d_start]),\
  631. np.angle(df_sess.shufflepred_complex_curr_combined[i][d_start]))\
  632. for i in df_sess.index]
  634. mono_prederr['early'][df_sess.index] = np.squeeze(prederr)
  635. df_sess['prederr_currdelay_start'] = np.squeeze(prederr)
  637. # end of delay
  638. prederr = [circdist(np.angle(df_sess.pred_complex_curr_combined[i][d_end]),\
  639. np.angle(df_sess.shufflepred_complex_curr_combined[i][d_end]))\
  640. for i in df_sess.index]
  641. mono_prederr['late'][df_sess.index] = np.squeeze(prederr)
  642. df_sess['prederr_currdelay_end'] = np.squeeze(prederr)
  645. if folded==True: # flip error, rel_loc in case of folded
  646. df_sess['prederr_currdelay_start'] = df_sess['prederr_currdelay_start'].values*\
  647. df_sess.sign_prevcurr
  648. df_sess['prederr_currdelay_end'] = df_sess['prederr_currdelay_end'].values*\
  649. df_sess.sign_prevcurr
  651. # get mean PREDICTION error at each relative location for each session, delay split
  652. for loc_id,loc in enumerate(prevcurr):
  653. decodersd[session_id, 0, loc_id] = circmean(df_sess.loc[df_sess.prev_curr==loc]['prederr_currdelay_start'],\
  654. low=-np.pi, high=np.pi)
  655. decodersd[session_id, 1, loc_id] = circmean(df_sess.loc[df_sess.prev_curr==loc]['prederr_currdelay_end'],\
  656. low=-np.pi, high=np.pi)
  657. # plt.hist(df_sess['prederr_currdelay'+str(delay_id)].values, alpha=0.3)
  658. # plt.show()
  659. # get session-mean in each relative location (nans per session if the target didn't appear)
  660. mean = np.nanmean(decodersd, axis=0)
  661. # TODO change errorbars
  662. errorbars = 'SEM'
  663. if errorbars=='SEM':
  664. std = sem(decodersd, axis=0, nan_policy='omit')
  665. elif errorbars == "CI":
  666. std = 2*np.nanstd(decodersd, axis=0)
  667. # make consecutive colors for diff delays
  668. label = ['early', 'late']
  669. xx = np.linspace(-np.pi, np.pi, 1000)
  671. f, ax = plt.subplots(figsize=(2.1,1.95))
  672. plt.axhline(0, color='#333333', alpha=0.5)
  673. plt.axvline(0, color='#333333', alpha=0.5)
  674. plt.errorbar(np.rad2deg(prevcurr), np.rad2deg(mean[0]), yerr=np.rad2deg(std[0]),\
  675. color=colors['SerialBiasWeak'], label='stim.')
  676. plt.errorbar(np.rad2deg(prevcurr), np.rad2deg(mean[1]), yerr=np.rad2deg(std[1]),\
  677. color=colors['SerialBias'], label='delay')
  679. # fit DoG
  680. df_mono['neural_error'] = mono_prederr[label[0]]
  681. df_mono['neural_error_late'] = mono_prederr[label[1]]
  682. para = neural_sigma[monkey]
  683. df_mono['bias_estimate'] = -hf.dog1(para, df_mono.prev_curr)
  685. # fit stimulus
  686. model = smf.ols('neural_error ~ bias_estimate', data=df_mono).fit()
  687. plt.plot(np.rad2deg(xx), np.rad2deg(model.params['Intercept']+\
  688. model.params['bias_estimate']*-hf.dog1(para, xx)),\
  689. color=colors['SerialBiasWeak'], dashes=[1,1])
  691. # reporting summary EARLY
  692. ci_low, ci_high = model.conf_int().loc['bias_estimate']
  693. print('beta = '+str(np.round(np.rad2deg(model.params['bias_estimate']),2))+\
  694. ', 95% CI ['+str(np.round(np.rad2deg(ci_low),2))+', '+str(np.round(np.rad2deg(ci_high),2))+'], '+\
  695. 't('+str(int(model.df_resid))+') = '+str(np.round(model.tvalues['bias_estimate'],2))+\
  696. ', p = '+str("{:.2e}".format(model.pvalues['bias_estimate'])) +\
  697. ', # trials: '+str(len(df_mono))+', '+str(len(df_mono.session.unique()))+' sessions.')
  700. # fit delay
  701. model = smf.ols('neural_error_late ~ bias_estimate', data=df_mono).fit()
  702. plt.plot(np.rad2deg(xx), np.rad2deg(model.params['Intercept']+\
  703. model.params['bias_estimate']*-hf.dog1(para, xx)),\
  704. color=colors['SerialBias'], dashes=[1,1])
  706. # reporting summary LATE
  707. ci_low, ci_high = model.conf_int().loc['bias_estimate']
  708. print('beta = '+str(np.round(np.rad2deg(model.params['bias_estimate']),2))+\
  709. ', 95% CI ['+str(np.round(np.rad2deg(ci_low),2))+', '+str(np.round(np.rad2deg(ci_high),2))+'], '+\
  710. 't('+str(int(model.df_resid))+') = '+str(np.round(model.tvalues['bias_estimate'],2))+\
  711. ', p = '+str("{:.2e}".format(model.pvalues['bias_estimate'])) +\
  712. ', # trials: '+str(len(df_mono))+', '+str(len(df_mono.session.unique()))+' sessions.')
  714. plt.xlabel('rel. previous location (°)')
  715. plt.ylabel('decoder err. (°)')
  716. plt.legend()
  717. sns.despine()
  718. ax.xaxis.set_ticks_position('bottom')
  719. ax.yaxis.set_ticks_position('left')
  721. if monkey=='Sa':
  722. plt.ylim([-5, 5])
  723. plt.tight_layout()
  724. #plt.savefig(DATAPATH+'./Figures/Figure5/INLAYSerialBiasDrift_NeuralFitDoG_'+monkey+'.svg')
  725. sns.despine()
  726. plt.show()
  731. # %% [markdown]
  732. # ### Fig 5b: neural SD
  734. # %%
  735. # compute correlations of neuron and behavior
  736. # for each session get a correlation of neurons / behavior separately
  737. sd_all = {m: [] for m in monkeys}
  738. sd_joint = []
  739. SD_slopes = []
  741. borders_all = [sess_df.borders_full.iloc[0] for monkey in ['Sa', 'Wa']\
  742. for _, sess_df in df[df.monkey == monkey].groupby('session')]
  743. borders_mean = cut_task_timings(timeperiods = borders_all, max_time = 19)[1]
  744. min_stimDelay = borders_mean[16] - borders_mean[14]
  745. for m,monkey in enumerate(['Sa', 'Wa']):#enumerate(['Sa','Wa']):#, 'Wa']):#
  746. print('Computing monkey '+monkey+'...')
  747. df_mono = df.loc[df.monkey==monkey]
  748. # compute SD for sessions separately
  749. for session_id,session in enumerate(np.unique(df_mono.session)):#
  750. df_sess = df_mono.loc[(df_mono.session==session)].copy().reset_index(drop=True)
  752. df_sess['bias_estimate'] = -hf.dog1(neural_sigma[monkey], df_sess.prev_curr)
  753. # compute DECODER SD for different time points in the CURRENT delay
  754. sd_time=[]
  755. for delay_id,delay in enumerate(range(len(df_sess.pred_complex_curr_combined[1]))):
  756. # get average prediction error in defined timesteps
  757. df_sess['prederr_delaystep'] = circdist([np.angle(df_sess.pred_complex_curr_combined[i][delay])\
  758. for i in df_sess.index],\
  759. [np.angle(df_sess.shufflepred_complex_curr_combined[i][delay])\
  760. for i in df_sess.index])
  762. # fit model
  763. model = smf.ols('prederr_delaystep ~ bias_estimate', data=df_sess).fit()
  765. #save parameter of model
  766. sd_time.append(np.rad2deg(model.params['bias_estimate']))
  767. sd_joint.append(sd_time)
  768. sd_all[monkey].append(sd_time)
  770. # fit slopr through time
  771. # select timing (from stim-start of each session to minimum delay length)
  772. x_coords = np.array(range(df_sess.borders_full[0][14], df_sess.borders_full[0][14] + min_stimDelay))*bins/1000
  773. sd_stim2end = sd_time[df_sess.borders_full[0][14]:df_sess.borders_full[0][14] + min_stimDelay]
  774. slopel, intercept, r_value, p_value, std_err = stats.linregress(x_coords, sd_stim2end)
  775. SD_slopes.append(slopel)
  778. # cut sessions to same length
  780. sd_cut,borders_mean = cut_task_timings(sd_joint, borders_all, 19)
  781. # only look at 2nd trial's delay
  782. sd_trial2 = sd_cut[:, borders_mean[11]:]
  784. plot_twolines_full(R=sd_trial2,base=np.zeros((sd_trial2.shape)), bins=bins,\
  785. labelR='',labelB='', errorbars='SEM', end_border=6,shorten_delay = 2, currentTrial=True,\
  786. borders=borders_mean, ylabel='history drift (°)', colors=[colors['SerialBias'], '#333333'])
  787. plt.tight_layout()
  788. #plt.savefig('./Figures/Figure5/SerialBiasDelayDrift_SaWa_200ms.svg')
  789. plt.show()
  792. ## STATISTICAL TESTING
  794. print('Slope computed on '+str(min_stimDelay)+' independent '+str(bins)+'ms time bins.')
  795. degf = len(SD_slopes) - 1
  796. test, pval = ttest_1samp(SD_slopes, 0)
  797. mean_x = np.mean(SD_slopes)
  798. se = stats.sem(SD_slopes) # standard error of the mean
  799. ci_low, ci_high = t.interval(0.95, degf, loc=mean_x, scale=se)
  800. print('SD slope 1-sample: t('+str(degf)+') = '+str(np.round(test, 2))+', p = '+"{:.2e}".format(pval)+\
  801. ", mean = "+str(np.round(mean_x))+\
  802. ", 95% CI ["+str(np.round(ci_low,2))+","+str(np.round(ci_high, 2))+"]")
  804. # %% [markdown]
  805. # ### Supplement: Correlate behavioral SD vs neural SD
  807. # %%
  808. from sklearn.model_selection import KFold
  810. # decoder error of delay_n+1 to target_n+1
  811. folded=False
  813. # define delay timesteps
  814. delay_steps = int(600/bins)# size of each delay step
  816. colors_SD = ['#254441ff', '#43aa8bff', '#b2b09bff']
  818. all_behavfits=[]
  819. all_neuralfits=[]
  820. f, ax = plt.subplots(figsize=(2.,1.95))
  821. ax.axhline(0, color='#333333', alpha=0.5, lw=0.5)
  822. ax.axvline(0, color='#333333', alpha=0.5, lw=0.5)
  823. for m,monkey in enumerate(monkeys):#enumerate(['Pe','Sa', 'Wa']):#
  824. borders_mono=[]
  825. df_mono = df.loc[(df.monkey==monkey)].copy().reset_index(drop=True)
  826. df_mono['sign_prevcurr'] = sign_rl(df_mono.prev_curr.values)
  827. if folded==True:
  828. df_mono['prev_curr'] = np.abs(df_mono.prev_curr)
  829. prevcurr = np.unique(df_mono.loc[df_mono.session==0].prev_curr)#np.unique(df_mono.prev_curr) max_num_delaysteps = np.min((np.array([df_mono.borders_full.values[i][16]-df_mono.borders_full.values[i][15]\
  832. # SD FITS ARE OPTIMIZED DIFFERENTLY FOR NEURAL, BEHAV. SD
  833. df_mono['behav_error_curr_deg'] = np.rad2deg(df_mono.behav_error_curr)
  834. df_mono['behav_bias_estimate'] = -hf.dog1(sigma[monkey], df_mono.prev_curr)#sign_rl(df_sess.prev_curr.values)#
  835. df_mono['neural_bias_estimate'] = -hf.dog1(neural_sigma[monkey], df_mono.prev_curr)#sign_rl(df_sess.prev_curr.values)#
  838. mono_prederr = np.empty((len(df_mono)))*np.nan
  839. # compute SD for sessions separately
  840. for session_id,session in enumerate(np.unique(df_mono.session)):#
  841. df_sess = df_mono.loc[(df_mono.session==session)].copy()
  843. borders_full = df_sess.borders_full.values[0]
  844. d_end = range(borders_full[17]-2, borders_full[17])
  845. #d_end = borders_full[17]-1
  847. #mono_prederr[df_sess.index] = np.squeeze([circdist(np.angle(df_sess.pred_complex_curr_combined[i][d_end]),\
  848. # np.angle(df_sess.shufflepred_complex_curr_combined[i][d_end]))\
  849. # for i in df_sess.index])
  850. mono_prederr[df_sess.index] = np.squeeze([circmean(circdist(np.angle(df_sess.pred_complex_curr_combined[i][d_end]),\
  851. np.angle(df_sess.shufflepred_complex_curr_combined[i][d_end])),\
  852. low=-np.pi, high=np.pi)\
  853. for i in df_sess.index])
  854. df_mono['prederr'] = mono_prederr
  856. ################ FIT SERIAL DEPENDENCE ##########################
  857. kf = KFold(n_splits=5, shuffle=True)
  858. behavfits=[]
  859. neuralfits=[]
  860. for k, (train_idx, test_idx) in enumerate(kf.split(df_mono.index)):
  861. df_split = df_mono.loc[train_idx]
  863. # fit BEHAVIORAL model
  864. #behavmodel = smf.mixedlm('behav_error_curr_deg ~ behav_bias_estimate', data=df_mono, groups=df_mono['session']).fit()
  865. behavmodel = smf.ols('behav_error_curr_deg ~ behav_bias_estimate', data=df_split).fit()
  866. behavfits.append(behavmodel.params[1])
  868. # fit NEURAL model# fit BEHAVIORAL model
  869. df_split['prederr_end_deg'] = np.rad2deg(df_split['prederr'])
  870. neuralmodel = smf.ols('prederr_end_deg ~ neural_bias_estimate', data=df_split).fit()
  871. #if (behavmodel.pvalues[1]<0.05) & (neuralmodel.pvalues[1]<0.05):
  872. behavfits.append(behavmodel.params[1])
  873. neuralfits.append(neuralmodel.params[1])
  874. all_behavfits.append(behavfits)
  875. all_neuralfits.append(neuralfits)
  877. ax.errorbar(np.mean(all_behavfits[m]), np.mean(all_neuralfits[m]),\
  878. xerr= sem(all_behavfits[m]),yerr=2*sem(all_neuralfits[m]), color=colors_SD[m], marker='o',\
  879. label=monkey, markersize=4)
  880. ax.set_xlabel('behavioral bias (°)')
  881. ax.set_ylabel('history drift (°)')
  883. x0, x1 = ax.get_xlim()
  884. #ax.set_xlim([x0, x1])
  885. #ax.set_ylim([x0, x1])
  887. y0, y1 = ax.get_ylim()
  888. ax.plot([y0, y1], [y0, y1], color='#333333', alpha=0.3, dashes=[5,5], lw=0.5)
  890. ax.legend()
  891. ax.xaxis.set_ticks_position('bottom')
  892. ax.yaxis.set_ticks_position('left')
  893. sns.despine()
  894. plt.tight_layout()
  895. #plt.savefig('./NeuralBehavSerialBias_compareMonkeys.svg', dpi=300)
  896. plt.show()
  901. # %% [markdown]
  902. # ---
  904. # %% [markdown]
  905. # # Reactivation drift
  907. # %%
  908. # get minimum trial durations across all sessions
  909. borders_all = [sess_df.borders_full.iloc[0] for monkey in ['Sa', 'Wa']\
  910. for _, sess_df in df[df.monkey == monkey].groupby('session')]
  911. borders_mean = cut_task_timings(timeperiods = borders_all, max_time = 19)[1]
  913. # get reactivation location / strength for each hemisphere
  914. min_Reactlen = (borders_mean[14]+1) - borders_mean[13]
  915. print('Reactivation average is averaged across '+str(min_Reactlen)+' independent '+str(bins)+'ms time bins.')
  916. Rtime = [range(df.borders_full[i][13], df.borders_full[i][13]+min_Reactlen) for i in df.index]
  917. df['react'] = [np.angle(np.mean(df['pred_complex_delay_combined'][i][Rtime[i]]))\
  918. for i in df.index]
  919. df['react_left'] = [np.angle(np.mean(df['pred_complex_delay_left'][i][Rtime[i]]))\
  920. for i in df.index]
  921. df['react_right'] = [np.angle(np.mean(df['pred_complex_delay_right'][i][Rtime[i]]))\
  922. for i in df.index]
  923. df['react_strength'] = [np.abs(np.mean(df['pred_complex_delay_combined'][i][Rtime[i]]))\
  924. for i in df.index]
  925. df['react_strength_left'] = [np.abs(np.mean(df['pred_complex_delay_left'][i][Rtime[i]]))\
  926. for i in df.index]
  927. df['react_strength_right'] = [np.abs(np.mean(df['pred_complex_delay_right'][i][Rtime[i]]))\
  928. for i in df.index]
  931. # get delay angle, shuffle for each hemisphere
  932. min_Delaylen = borders_mean[17] - borders_mean[16]
  933. DELAYEND = [range(df.borders_full[trial][17] - min_Delaylen, df.borders_full[trial][17])\
  934. for trial in df.index]
  935. print('Delay average is averaged across '+str(min_Delaylen)+' independent '+str(bins)+'ms time bins.')
  936. df['delayAvg'] = [np.angle(np.mean(df.pred_complex_curr_combined[trial][DELAYEND[trial]]))\
  937. for trial in df.index]
  938. df['delayAvg_shuffle'] = [np.angle(np.mean(df.shufflepred_complex_curr_combined[trial][DELAYEND[trial]]))\
  939. for trial in df.index]
  940. df['delayAvg_left'] = [np.angle(np.mean(df.pred_complex_curr_left[trial][DELAYEND[trial]]))\
  941. for trial in df.index]
  942. df['delayAvg_shuffle_left'] = [np.angle(np.mean(df.shufflepred_complex_curr_left[trial][DELAYEND[trial]]))\
  943. for trial in df.index]
  944. df['delayAvg_right'] = [np.angle(np.mean(df.pred_complex_curr_right[trial][DELAYEND[trial]]))\
  945. for trial in df.index]
  946. df['delayAvg_shuffle_right'] = [np.angle(np.mean(df.shufflepred_complex_curr_right[trial][DELAYEND[trial]]))\
  947. for trial in df.index]
  950. # %% [markdown]
  951. # ### Fig 5c: Reactivation strength and precision
  953. # %%
  954. df_mono = df.loc[df.monkey=='Sa'].copy().reset_index(drop=True)
  956. # compute reactivation error (distance to previous target)
  957. df_mono['react_error'] = circdist(df_mono['react'].values, np.round(np.angle(df_mono.targ_prev_xy), 3))
  959. # determine cut-off for each hemisphere, session
  960. perc = 20
  961. df_mono[['cut']] = df_mono.groupby('session_continuous') \
  962. ['react_strength'].transform(lambda x: np.percentile(x, perc))
  963. df_high = df_mono.loc[df_mono.react_strength > df_mono.cut]
  964. df_low = df_mono.loc[df_mono.react_strength <= df_mono.cut]
  966. f,ax = plt.subplots(figsize=(2.1,1.95))
  967. ax.hist(np.rad2deg(df_high['react_error'].values), label='strong',\
  968. bins=17, density=True, histtype='step', color= colors['Reactivation'], linewidth= 1)
  969. ax.hist(np.rad2deg(df_low['react_error'].values), label='weak',\
  970. bins=17, density=True, histtype='step', color=colors['ReactivationWeak'], linewidth= 1)
  971. ax.legend()
  972. ax.set_xlim([-180, 180])
  973. ax.set_xlabel('reactivation error (°)')
  974. ax.set_ylabel('density')
  975. ax.xaxis.set_ticks_position('bottom')
  976. ax.yaxis.set_ticks_position('left')
  977. sns.despine()
  978. plt.tight_layout()
  979. #plt.savefig('./ReactivationAngleHist_Sa.svg', dpi=300)
  980. plt.show()
  982. # %% [markdown]
  983. # ### Fig 5d: Within vs across attraction to reactivation
  985. # %%
  986. # compute correlations of neuron and behavior
  987. # for each session get a correlation of neurons / behavior separately
  988. hemispheres = np.array(['left', 'right'])
  989. react_all = {'within':[], 'across':[]}
  990. react_term = {m: {side: [] for side in ['within', 'across']} for m in monkeys}
  991. borders_mono={m:[] for m in monkeys}
  992. borders_all=[]
  993. for m,monkey in enumerate(['Sa', 'Wa']):#enumerate(['Sa', 'Wa']):#
  994. print('Computing monkey '+monkey+'...')
  995. df_mono = df.loc[df.monkey==monkey]
  996. for session_id,session in enumerate(np.unique(df_mono.session)):#
  997. df_sess = df_mono.loc[(df_mono.session==session)].copy().reset_index(drop=True)
  998. borders_all.append(df_sess.borders_full[0])
  999. borders_mono[monkey].append(df_sess.borders_full[0])
  1001. # get relative distance: reactivation - current target
  1002. reactivationprevcurr_left = circdist(df_sess.react_left, np.angle(df_sess.targ_curr_xy))
  1003. df_sess['bias_left'] = -hf.dog1(reactivation_sigma[monkey],\
  1004. reactivationprevcurr_left)
  1005. reactivationprevcurr_right = circdist(df_sess.react_right, np.angle(df_sess.targ_curr_xy))
  1006. df_sess['bias_right'] = -hf.dog1(reactivation_sigma[monkey],\
  1007. reactivationprevcurr_right)
  1009. # subtract shuffle from delay error for all times
  1010. df_sess['pred_shufflesub_left'] = [circdist(np.angle(df_sess['pred_complex_curr_left'][i]),\
  1011. np.angle(df_sess.shufflepred_complex_curr_left[i]))\
  1012. for i in df_sess.index]
  1013. df_sess['pred_shufflesub_right'] = [circdist(np.angle(df_sess['pred_complex_curr_right'][i]),\
  1014. np.angle(df_sess.shufflepred_complex_curr_right[i]))\
  1015. for i in df_sess.index]
  1017. for s, same in enumerate(['within', 'across']):
  1018. sess_estimates=[]
  1019. # get estimates for each time step
  1020. for delay_id,delay in enumerate(range(len(df_sess.pred_complex_curr_combined[1]))):
  1021. bias_estimate = []
  1022. prederr_delaystep=[]
  1023. timeon=time.time()
  1024. for R,RSide in enumerate(hemispheres):
  1025. DSide = RSide if same == 'within' else hemispheres[hemispheres != RSide][0]
  1027. bias_estimate.append(df_sess['bias_'+RSide])
  1028. prederr_delaystep.append([df_sess['pred_shufflesub_'+DSide][i][delay]\
  1029. for i in df_sess.index])
  1030. #prederr_delaystep.append(np.squeeze([circdist(np.angle(df_sess['pred_complex_curr_'+DSide][i][delay]),\
  1031. # np.angle(df_sess.targ_curr_xy[i]))\
  1032. # for i in df_sess.index]))
  1034. # fit model on combined areas react/delay err
  1035. df_model = pd.DataFrame({'bias_estimate': np.concatenate(bias_estimate),
  1036. 'prederr_delaystep_deg': np.rad2deg(np.concatenate(prederr_delaystep))})
  1037. assert len(df_model) == 2*len(df_sess), "Using both left and right predictions needs to double trials (e.g. within = left-left, right-right)"
  1039. # fit model
  1040. model = smf.ols('prederr_delaystep_deg ~ bias_estimate', data=df_model).fit()
  1042. sess_estimates.append(model.params['bias_estimate'])
  1044. react_all[same].append(sess_estimates)
  1045. react_term[monkey][same].append(sess_estimates)
  1047. # cut sessions to same length
  1048. react_within_cut,borders_mean = cut_task_timings(react_all['within'], borders_all, 19)
  1049. react_across_cut,borders_mean = cut_task_timings(react_all['across'], borders_all, 19)
  1050. # only look at 2nd trial's delay
  1052. react_within_trial2 = react_within_cut[:, borders_mean[11]:]
  1053. react_across_trial2 = react_across_cut[:, borders_mean[11]:]
  1055. plot_twolines_full(R=react_within_trial2,base=react_across_trial2, bins=bins, currentTrial=True,\
  1056. labelR='within',labelB='across', errorbars='SEM', end_border=6,shorten_delay = 2,\
  1057. borders=borders_mean, ylabel='reactivation drift (°)',\
  1058. significances = True)
  1059. plt.tight_layout()
  1060. #plt.savefig('./Figures/Figure5/ReactivationBiasWithinAcross_SaWa_200ms.svg')
  1061. plt.show()
  1064. # %%
  1065. react_within_cut,borders_mean = cut_task_timings(react_all['within'], borders_all, 19)
  1066. react_across_cut,borders_mean = cut_task_timings(react_all['across'], borders_all, 19)
  1068. within_delayavg = np.mean(react_within_cut[:, borders_mean[16]:borders_mean[17]], axis=1)
  1069. across_delayavg = np.mean(react_across_cut[:, borders_mean[16]:borders_mean[17]], axis=1)
  1071. # plot
  1072. f, ax = plt.subplots(figsize=(1.1,2.1))
  1073. plt.axhline(0, color='#333333', alpha=0.5)
  1074. #within
  1075. ax.plot([np.zeros((len(within_delayavg))), np.ones((len(across_delayavg)))],\
  1076. [within_delayavg, across_delayavg], color='#333333', alpha=0.2)
  1077. ax.errorbar(0, np.mean(within_delayavg), yerr= 2*sem(within_delayavg), color=colors['Within'], marker='o')
  1079. # across
  1080. ax.errorbar(1, np.mean(across_delayavg), yerr= 2*sem(across_delayavg), color=colors['Across'], marker='o')
  1081. ax.set_xticks([0,1])
  1082. ax.set_xticklabels(['within', 'across'], rotation=25)
  1083. ax.set_ylabel('reactivation bias (°)')
  1084. ax.set_xlabel('hemisphere')
  1086. # pvalues
  1087. stats_w = ['***' if ttest_1samp(within_delayavg, 0)[1] < 0.005\
  1088. else '**' if ttest_1samp(within_delayavg, 0)[1] < 0.01\
  1089. else '*' if ttest_1samp(within_delayavg, 0)[1] < 0.05 else 'n.s.'][0]
  1090. stats_a = ['***' if ttest_1samp(across_delayavg, 0)[1] < 0.005\
  1091. else '**' if ttest_1samp(across_delayavg, 0)[1] < 0.01\
  1092. else '*' if ttest_1samp(across_delayavg, 0)[1] < 0.05 else 'n.s.'][0]
  1093. stats_diff = ['***' if ttest_rel(within_delayavg, across_delayavg)[1] < 0.005\
  1094. else '**' if ttest_rel(within_delayavg, across_delayavg)[1] < 0.01\
  1095. else '*' if ttest_rel(within_delayavg, across_delayavg)[1] < 0.05 else 'n.s.'][0]
  1097. y0, y1 = ax.get_ylim()
  1098. ax.annotate(stats_w,(0, y1-0.3), ha='center', va='bottom', fontsize=10, weight='bold', color='#333333')
  1099. ax.annotate(stats_a,(1, y1-0.3), ha='center', va='bottom', fontsize=10, weight='bold', color='#333333')
  1100. ax.plot([0,0,1,1],\
  1101. [y1+2*y1/10,y1+3*y1/10, y1+3*y1/10, y1+2*y1/10], linewidth=1, color='#333333')
  1102. ax.annotate(stats_diff,\
  1103. (0.5, y1+3*y1/10), ha='center', va='bottom', fontsize=10, weight='bold', color='#333333')
  1104. ax.set_xlim([-.2,1.2])
  1105. ax.set_ylim([-3.5,8])
  1106. ax.xaxis.set_ticks_position('bottom')
  1107. ax.yaxis.set_ticks_position('left')
  1108. sns.despine()
  1109. plt.tight_layout()
  1110. #plt.savefig('./Quantify_ReactivationBiasWithinAcross_SaWa_200ms.svg')
  1111. plt.show()
  1114. # STATISTICAL TESTING
  1115. print('Avg. of '+str(borders_mean[17] - borders_mean[16])+' independent '+str(bins)+' ms time bins.')
  1116. # WITHIN
  1117. degf = len(within_delayavg) - 1
  1118. test, pval = ttest_1samp(within_delayavg, 0)
  1119. mean_x = np.mean(within_delayavg)
  1120. se = stats.sem(within_delayavg) # standard error of the mean
  1121. ci_low, ci_high = t.interval(0.95, degf, loc=mean_x, scale=se)
  1122. print('WITHIN 1-sample: t('+str(degf)+') = '+str(np.round(test, 2))+', p = '+"{:.2e}".format(pval)+\
  1123. ', mean = '+str(np.round(mean_x,2))+\
  1124. ", 95% CI ["+str(np.round(ci_low,2))+","+str(np.round(ci_high, 2))+"]")
  1126. # ACROSS
  1127. test, pval = ttest_1samp(across_delayavg, 0)
  1128. mean_x = np.mean(across_delayavg)
  1129. se = stats.sem(across_delayavg) # standard error of the mean
  1130. ci_low, ci_high = t.interval(0.95, degf, loc=mean_x, scale=se)
  1131. print('ACROSS 1-sample: t('+str(degf)+') = '+str(np.round(test, 2))+', p = '+"{:.2e}".format(pval)+\
  1132. ', mean = '+str(np.round(mean_x,2))+\
  1133. ", 95% CI ["+str(np.round(ci_low,2))+","+str(np.round(ci_high, 2))+"]")
  1135. # WITHIN - ACROSS
  1136. test, pval = ttest_rel(within_delayavg, across_delayavg)
  1137. mean_x = np.mean(within_delayavg - across_delayavg)
  1138. se = stats.sem(within_delayavg - across_delayavg) # standard error of the mean
  1139. ci_low, ci_high = t.interval(0.95, degf, loc=mean_x, scale=se)
  1140. print('WITHIN - ACROSS paired: t('+str(degf)+') = '+str(np.round(test, 2))+', p = '+"{:.2e}".format(pval)+\
  1141. ', mean = '+str(np.round(mean_x,2))+\
  1142. ", 95% CI ["+str(np.round(ci_low,2))+","+str(np.round(ci_high, 2))+"]")
  1145. # %% [markdown]
  1146. # ### Supplement: Split by percentage of dual trials
  1148. # %%
  1149. df_helper = df.copy().reset_index(drop=True)
  1150. colors_wa = {'within': 'darkgreen', 'across':'darkorange'}
  1151. hemispheres = np.array(['left', 'right'])
  1153. percent = 30
  1154. print('Cut-off: '+str(percent)+'%...')
  1155. # determine cut-off for each session
  1156. df_helper[['cut_left', 'cut_right']] = df_helper.groupby('session_continuous') \
  1157. [['react_strength_left', 'react_strength_right']] \
  1158. .transform(lambda x: np.percentile(x, percent))
  1160. # mark reactivations as weak, strong
  1161. df_helper['R_strong_left'] = df_helper.react_strength_left > df_helper.cut_left
  1162. df_helper['R_strong_right'] = df_helper.react_strength_right > df_helper.cut_right
  1164. # assign reactivation types
  1165. chosen_trials = {'none': df_helper[(~df_helper.R_strong_left) & (~df_helper.R_strong_right)],\
  1166. 'single': df_helper[((df_helper.R_strong_left) & (~df_helper.R_strong_right)) |\
  1167. ((~df_helper.R_strong_left) & (df_helper.R_strong_right))],\
  1168. 'dual': df_helper[(df_helper.R_strong_left) & (df_helper.R_strong_right)]}
  1169. ['within', 'across']
  1170. sd_monkeys = {'within':[],'across':[]}
  1171. sd_term = {key:{'within': np.empty((len(np.unique(df_helper.session_continuous))))*np.nan,\
  1172. 'across': np.empty((len(np.unique(df_helper.session_continuous))))*np.nan}\
  1173. for key in chosen_trials.keys()}
  1174. pid=-1
  1175. f,ax = plt.subplots(figsize=(1.85,1.95), sharey=True)
  1176. for x, react_combi in enumerate(chosen_trials.keys()):
  1177. df_type = chosen_trials[react_combi]
  1178. sides=[]
  1179. for s, same in enumerate(['within', 'across']):
  1180. for session_id,session in enumerate(df_type.session_continuous.unique()):#
  1181. df_sess = df_type.loc[(df_type.session_continuous==session)].copy()
  1182. monkey = df_sess.monkey.unique()[0]
  1184. sess_estimates, reactivationprev_curr, delayerr = [], [], []
  1185. # for each hemisphere
  1186. for R,RSide in enumerate(hemispheres):
  1187. # use either same hemisphere during delay (within) or opposite (across)
  1188. DSide = RSide if same == 'within' else hemispheres[hemispheres != RSide][0]
  1190. reactivations = df_sess['react_'+RSide].values
  1191. delay = df_sess['delayAvg_'+DSide].values
  1192. shuffledelay = df_sess['delayAvg_shuffle_'+DSide].values# np.angle(df_sess.targ_curr_xy)
  1193. # RELATIVE LOCATION OF PREVIOUS REACTIVATION TO CURRENT TARGET
  1194. if (react_combi == 'single'): # if single, only view attraction to reactivated side
  1195. df_react_side = df_sess.loc[df_sess['R_strong_'+RSide]]
  1196. reactivationprev_curr.append(circdist(df_react_side['react_'+RSide].values,\
  1197. np.angle(df_react_side.targ_curr_xy)))
  1199. delay_react_side = df_react_side['delayAvg_'+DSide].values
  1200. shuffle_react_side = df_react_side['delayAvg_shuffle_'+DSide].values
  1201. delayerr.append(circdist(delay_react_side, shuffle_react_side))
  1202. else:
  1203. reactivationprev_curr.append(circdist(reactivations, np.angle(df_sess.targ_curr_xy)))
  1204. delayerr.append(circdist(delay, shuffledelay))
  1206. # combine left right in each condition
  1207. reactivationprev_curr = np.concatenate(reactivationprev_curr)
  1208. delayerr = np.concatenate(delayerr)
  1210. # MODEL HOW THE DELAY ACTIVITY IS ATTRACTED TO THE PREVIOUS REACTIVATION LOCATION
  1211. bias_estimate = -hf.dog1(reactivation_sigma[monkey], reactivationprev_curr)#sign_rl(reactivationprev_curr)#
  1213. # Prepare DataFrame for model fitting
  1214. df_model = pd.DataFrame({'bias_estimate': bias_estimate,
  1215. 'prederr_delaystep_deg': np.rad2deg(delayerr)})
  1216. # fit model
  1217. model = smf.ols('prederr_delaystep_deg ~ bias_estimate', data=df_model).fit()
  1219. #save parameter of model
  1220. sd_term[react_combi][same][session_id] = model.params['bias_estimate']
  1221. sess_estimates.append(model.params['bias_estimate'])
  1222. sd_monkeys[same].append(sess_estimates)
  1225. # get session-mean in each relative location (nans per session if the target didn't appear)
  1227. sd_timing = sd_term[react_combi][same]
  1228. mean, CI = np.nanmean(sd_timing, axis=0), sem(sd_timing, axis=0, nan_policy='omit')
  1230. ax.errorbar(x+s*0.2, mean, yerr=CI, color=colors[same.title()], marker='o', label=same)
  1231. #if (react_combi=='none'):
  1232. # plt.legend()
  1233. ax.axhline(0, color='#333333', alpha=0.5)
  1235. # pvalues of each condition against 0
  1236. y0, y1 = ax.get_ylim()
  1237. for x, react_combi in enumerate(chosen_trials.keys()):
  1238. test_0 = ttest_rel(sd_term[react_combi]['within'], sd_term[react_combi]['across'])[1]
  1239. if test_0 < 0.05:
  1240. ax.plot([x,x,x+0.2,x+0.2],\
  1241. [y1-0.2,y1-0.1, y1-0.1, y1-0.2], linewidth=1, color='#333333')
  1242. ax.annotate('*',\
  1243. (x+0.1, y1-0.2), ha='center', va='bottom', fontsize=10, color='#333333')
  1244. # pvalues of single vs dual across condition (n.s.)
  1245. ttest_across = ttest_rel(sd_term['single']['across'], sd_term['dual']['across'])[1]
  1246. stats_a = ['***' if ttest_across < 0.005\
  1247. else '**' if ttest_across < 0.01\
  1248. else '*' if ttest_across < 0.05 else 'n.s.'][0]
  1249. ax.plot([1.2,1.2,2.2,2.2],\
  1250. [1.7,1.8, 1.8, 1.7], linewidth=1, color='#333333')
  1251. ax.annotate(stats_a,\
  1252. (1.7, 1.8), ha='center', va='bottom', fontsize=10, color='#333333')
  1254. plt.xticks([0,1,2], chosen_trials.keys())#, rotation=25
  1255. plt.xlabel('reactivation type')
  1256. plt.ylabel('react. drift (°)')
  1257. sns.despine()
  1258. ax.xaxis.set_ticks_position('bottom')
  1259. ax.yaxis.set_ticks_position('left')
  1260. plt.tight_layout()
  1261. #plt.savefig('./ReactivationBiasSplitByType_'+str(percent)+'%.svg')
  1262. plt.show()
  1265. # STATISTICAL TESTING
  1266. # WITHIN - ACROSS
  1267. print('##### Across-0 statistics: #####')
  1268. for x, react_combi in enumerate(chosen_trials.keys()):
  1269. test, pval = ttest_1samp(sd_term[react_combi]['across'], 0)
  1270. mean_x = np.mean(sd_term[react_combi]['across'])
  1271. se = stats.sem(sd_term[react_combi]['across']) # standard error of the mean
  1272. ci_low, ci_high = t.interval(0.95, degf, loc=mean_x, scale=se)
  1273. print(react_combi+': t('+str(degf)+') = '+str(np.round(test, 2))+', p = '+"{:.2e}".format(pval)+\
  1274. ', mean = '+str(np.round(mean_x,2))+\
  1275. ", 95% CI ["+str(np.round(ci_low,2))+", "+str(np.round(ci_high, 2))+"]")
  1277. print('##### Within-Across statistics: #####')
  1278. for x, react_combi in enumerate(chosen_trials.keys()):
  1279. test, pval = ttest_rel(sd_term[react_combi]['within'], sd_term[react_combi]['across'])
  1280. mean_x = np.mean(sd_term[react_combi]['within'] - sd_term[react_combi]['across'])
  1281. se = stats.sem(sd_term[react_combi]['within'] - sd_term[react_combi]['across']) # standard error of the mean
  1282. ci_low, ci_high = t.interval(0.95, degf, loc=mean_x, scale=se)
  1283. print(react_combi+': t('+str(degf)+') = '+str(np.round(test, 2))+', p = '+"{:.2e}".format(pval)+\
  1284. ', mean = '+str(np.round(mean_x,2))+\
  1285. ", 95% CI ["+str(np.round(ci_low,2))+", "+str(np.round(ci_high, 2))+"]")
  1287. print('##### Single vs Dual Across #####')
  1288. test, pval = ttest_rel(sd_term['dual']['across'], sd_term['single']['across'])
  1289. mean_x = np.mean(sd_term['dual']['within'] - sd_term['single']['across'])
  1290. se = stats.sem(sd_term['dual']['within'] - sd_term['single']['across']) # standard error of the mean
  1291. ci_low, ci_high = t.interval(0.95, degf, loc=mean_x, scale=se)
  1292. print('Single vs Dual Across: t('+str(degf)+') = '+str(np.round(test, 2))+', p = '+"{:.2e}".format(pval)+\
  1293. ', mean = '+str(np.round(mean_x,2))+\
  1294. ", 95% CI ["+str(np.round(ci_low,2))+", "+str(np.round(ci_high, 2))+"]")
  1296. print('##### Single vs Dual Within-Across Difference #####')
  1297. single_WA = sd_term['single']['within'] - sd_term['single']['across']
  1298. dual_WA = sd_term['dual']['within'] - sd_term['dual']['across']
  1299. test, pval = ttest_rel(single_WA, dual_WA)
  1300. mean_x = np.mean(dual_WA - single_WA)
  1301. se = stats.sem(dual_WA - single_WA) # standard error of the mean
  1302. ci_low, ci_high = t.interval(0.95, degf, loc=mean_x, scale=se)
  1303. print('Single vs Dual within-across: t('+str(degf)+') = '+str(np.round(test, 2))+', p = '+"{:.2e}".format(pval)+\
  1304. ', mean = '+str(np.round(mean_x,2))+\
  1305. ", 95% CI ["+str(np.round(ci_low,2))+", "+str(np.round(ci_high, 2))+"]")
  1307. # %% [markdown]
  1308. # ### Fig 5e: GMM - Single or Dual Reactivations
  1310. # %%
  1311. from sklearn.mixture import GaussianMixture
  1312. from sklearn.model_selection import train_test_split
  1313. from sklearn.model_selection import GridSearchCV
  1315. # get reactivation strengths
  1316. # get minimum trial durations across all sessions
  1317. borders_all = [sess_df.borders_full.iloc[0] for monkey in ['Sa', 'Wa']\
  1318. for _, sess_df in df[df.monkey == monkey].groupby('session')]
  1319. borders_mean = cut_task_timings(timeperiods = borders_all, max_time = 19)[1]
  1321. # get reactivation location / strength for each hemisphere
  1322. min_Reactlen = (borders_mean[14]+1) - borders_mean[13]
  1323. print('Reactivation average is averaged across '+str(min_Reactlen)+' independent '+str(bins)+'ms time bins.')
  1324. Rtimes = [range(df.borders_full[i][13], df.borders_full[i][13]+min_Reactlen) for i in df.index]
  1326. # Rtimes = [range(df.borders_full.values[trial][13], df.borders_full.values[trial][14]+1)\
  1327. # for trial in df.index]
  1328. # REACTIVATION
  1329. # get reactivations: predict PREV target during current fixation time (Rtimes)
  1330. df['react_strength_left'] = [np.abs(np.mean(df.pred_complex_delay_left[trial][Rtimes[trial]]))\
  1331. for trial in df.index]
  1332. df['react_strength_right'] = [np.abs(np.mean(df.pred_complex_delay_right[trial][Rtimes[trial]]))\
  1333. for trial in df.index]
  1334. df['react_angle_left'] = [np.angle(np.mean(df.pred_complex_delay_left[trial][Rtimes[trial]]))\
  1335. for trial in df.index]
  1336. df['react_angle_right'] = [np.angle(np.mean(df.pred_complex_delay_right[trial][Rtimes[trial]]))\
  1337. for trial in df.index]
  1339. # create vector of left, right react strength per trial
  1340. X = np.array([df.react_strength_left.values, df.react_strength_right.values]).T
  1342. # %%
  1343. # use four components, expect: None, single-left, single-right, Dual
  1344. chosen_components = 4
  1346. # get left, right reactivation strengths
  1347. df_mono = df.copy().reset_index(drop=True)
  1350. # only plot random subsample of datapoints to get smaller svg file
  1351. randomIDs = np.random.choice(range(len(X)), 2000)
  1353. # figure
  1354. f,ax=plt.subplots(figsize=(2.1,1.95))
  1355. X_subsample = X[randomIDs]
  1356. ax.scatter(X_subsample[:,0], X_subsample[:,1], alpha=0.2, color=colors['Reactivation'], marker='.')
  1357. # fit 250 GMMs on full data with differnt start points to check stability
  1358. gm = GaussianMixture(n_components=chosen_components)
  1359. for i in range(100): # repeat for random starting points
  1360. gm = GaussianMixture(n_components=chosen_components, covariance_type='full', max_iter=500).fit(X)#, covariance_type='spherical'
  1361. # highlight the found classes
  1362. # highlight the found classes
  1363. dual_class = np.argmax(np.sum(gm.means_, axis=1))
  1364. for m_id,m in enumerate(gm.means_):
  1365. #if m_id == dual_class :
  1366. # ax.scatter(m[0], m[1], color='r', marker='x', s=50, alpha=0.1)
  1367. #else:
  1368. ax.scatter(m[0], m[1], color='#333333', marker='x', s=50, alpha=0.1)
  1369. # plot diagonal for showing dual reactivation line
  1370. x0,x1 = ax.get_xlim()
  1371. y0,y1 = ax.get_ylim()
  1372. ax.plot(np.linspace(0,np.max([x1, y1]),10), np.linspace(0,np.max([x1, y1]),10),color='#333333', dashes=[1,1])
  1373. ax.set_xlabel('left reactivation strength')
  1374. ax.set_ylabel('right reactivation strength')
  1375. ax.xaxis.set_ticks_position('bottom')
  1376. ax.yaxis.set_ticks_position('left')
  1377. sns.despine()
  1378. plt.tight_layout()
  1379. #plt.savefig('./GMM4FitsSingleVSDualReact.svg')
  1380. plt.show()
  1382. # %% [markdown]
  1383. # ### Supplement: Individual monkeys with three classes
  1385. # %%
  1386. # use four components, expect: None, single-left, single-right, Dual
  1387. chosen_components = 4
  1389. f,ax=plt.subplots(1,3,figsize=(5.8,1.95), sharex=True, sharey=True)
  1390. for m, monkey in enumerate(monkeys):
  1391. # get left, right reactivation strengths
  1392. df_mono = df.loc[df.monkey==monkey].copy().reset_index(drop=True)
  1393. X_mono = np.array([df_mono.react_strength_left.values, df_mono.react_strength_right.values]).T
  1394. # only plot random subsample of datapoints to get smaller svg file
  1395. randomIDs = np.random.choice(range(len(X_mono)), 2000)
  1397. # figure
  1398. X_subsample = X_mono[randomIDs]
  1399. ax[m].scatter(X_subsample[:,0], X_subsample[:,1], alpha=0.2, color=colors['Reactivation'], marker='.')
  1400. # fit 250 GMMs on full data with differnt start points to check stability
  1401. gm = GaussianMixture(n_components=chosen_components)
  1402. for i in range(100): # repeat for random starting points
  1403. gm = GaussianMixture(n_components=chosen_components, covariance_type='full', max_iter=500).fit(X_mono)#, covariance_type='spherical'
  1404. # highlight the found classes
  1405. # highlight the found classes
  1406. dual_class = np.argmax(np.sum(gm.means_, axis=1))
  1407. for m_id,means in enumerate(gm.means_):
  1408. ax[m].scatter(means[0], means[1], color='#333333', marker='x', s=50, alpha=0.1)
  1409. # plot diagonal for showing dual reactivation line
  1410. x0,x1 = ax[m].get_xlim()
  1411. y0,y1 = ax[m].get_ylim()
  1412. ax[m].plot(np.linspace(0,np.max([x1, y1]),10), np.linspace(0,np.max([x1, y1]),10),color='#333333', dashes=[1,1])
  1413. ax[m].set_xlabel('left reactivation strength')
  1414. ax[m].set_ylabel('right reactivation strength')
  1415. ax[m].xaxis.set_ticks_position('bottom')
  1416. ax[m].yaxis.set_ticks_position('left')
  1417. ax[m].set_title(monkey)
  1418. sns.despine()
  1419. plt.tight_layout()
  1420. #plt.savefig('./SupplementGMM4FitsSingleVSDualReact.svg')
  1421. plt.show()
  1423. # %%
  1424. # use four components, expect: None, single-left, single-right, Dual
  1425. chosen_components = 3
  1427. f,ax=plt.subplots(1,3,figsize=(5.8,1.95), sharex=True, sharey=True)
  1428. for m, monkey in enumerate(monkeys):
  1429. # get left, right reactivation strengths
  1430. df_mono = df.loc[df.monkey==monkey].copy().reset_index(drop=True)
  1431. X_mono = np.array([df_mono.react_strength_left.values, df_mono.react_strength_right.values]).T
  1432. # only plot random subsample of datapoints to get smaller svg file
  1433. randomIDs = np.random.choice(range(len(X_mono)), 2000)
  1435. # figure
  1436. X_subsample = X_mono[randomIDs]
  1437. ax[m].scatter(X_subsample[:,0], X_subsample[:,1], alpha=0.2, color=colors['Reactivation'], marker='.')
  1438. # fit 250 GMMs on full data with differnt start points to check stability
  1439. gm = GaussianMixture(n_components=chosen_components)
  1440. for i in range(100): # repeat for random starting points
  1441. gm = GaussianMixture(n_components=chosen_components, covariance_type='full', max_iter=500).fit(X_mono)#, covariance_type='spherical'
  1442. # highlight the found classes
  1443. # highlight the found classes
  1444. dual_class = np.argmax(np.sum(gm.means_, axis=1))
  1445. for m_id,means in enumerate(gm.means_):
  1446. ax[m].scatter(means[0], means[1], color='#333333', marker='x', s=50, alpha=0.1)
  1447. # plot diagonal for showing dual reactivation line
  1448. x0,x1 = ax[m].get_xlim()
  1449. y0,y1 = ax[m].get_ylim()
  1450. ax[m].plot(np.linspace(0,np.max([x1, y1]),10), np.linspace(0,np.max([x1, y1]),10),color='#333333', dashes=[1,1])
  1451. ax[m].set_xlabel('left reactivation strength')
  1452. ax[m].set_ylabel('right reactivation strength')
  1453. ax[m].xaxis.set_ticks_position('bottom')
  1454. ax[m].yaxis.set_ticks_position('left')
  1455. ax[m].set_title(monkey)
  1456. sns.despine()
  1457. plt.tight_layout()
  1458. #plt.savefig('./SupplementGMM3FitsSingleVSDualReact.svg')
  1459. plt.show()
  1461. # %% [markdown]
  1462. # ### Supplement: Illustration of cut-off
  1464. # %%
  1465. percentile = 30
  1466. cut_off = np.percentile(X, percentile)
  1468. f,ax=plt.subplots(figsize=(2.,2))
  1469. plt.scatter(X_subsample[:,0], X_subsample[:,1], color='#333333', alpha=0.2, marker='.')
  1470. plt.axhline(cut_off, lw=1.25, dashes=[5,1], color='lightcoral')
  1471. plt.axvline(cut_off, lw=1.25, dashes=[5,1], color='lightcoral')
  1472. ax.set_xlabel('left reactivation strength')
  1473. ax.set_ylabel('right reactivation strength')
  1474. plt.xlim([0, 300])
  1475. plt.ylim([0, 300])
  1476. plt.title('All monkeys')
  1477. ax.xaxis.set_ticks_position('bottom')
  1478. ax.yaxis.set_ticks_position('left')
  1479. sns.despine()
  1480. plt.tight_layout()
  1481. #plt.savefig('./Figures/Supplement/Illustration_cutOff.svg')
  1483. # %% [markdown]
  1484. # ### Figure 5f: Delay GMM
  1486. # %%
  1487. borders_all = [sess_df.borders_full.iloc[0] for monkey in ['Sa', 'Wa']\
  1488. for _, sess_df in df[df.monkey == monkey].groupby('session')]
  1489. borders_mean = cut_task_timings(timeperiods = borders_all, max_time = 19)[1]
  1491. # get reactivation location / strength for each hemisphere
  1492. min_Delaylen = borders_mean[6] - borders_mean[5]
  1493. print('Delay average is averaged across '+str(min_Delaylen)+' independent '+str(bins)+'ms time bins.')
  1494. # get delay strengths
  1495. Dtimes = [range(df.borders_full.values[trial][6]-min_Delaylen, df.borders_full.values[trial][6])\
  1496. for trial in df.index]
  1497. # REACTIVATION
  1498. # get reactivations: predict PREV target during current fixation time (Rtimes)
  1499. df['delay_strength_left'] = [np.abs(np.mean(df.pred_complex_delay_left[trial][Dtimes[trial]]))\
  1500. for trial in df.index]
  1501. df['delay_strength_right'] = [np.abs(np.mean(df.pred_complex_delay_right[trial][Dtimes[trial]]))\
  1502. for trial in df.index]
  1503. df['delay_strength_combined'] = [np.abs(np.mean(df.pred_complex_delay_combined[trial][Dtimes[trial]]))\
  1504. for trial in df.index]
  1506. X_delay = np.array([df.delay_strength_left.values, df.delay_strength_right.values]).T
  1508. # %%
  1509. # use four components, expect: None, single-left, single-right, Dual
  1510. chosen_components = 4
  1512. # get left, right reactivation strengths
  1513. df_mono = df.copy().reset_index(drop=True)
  1515. # only plot random subsample of datapoints to get smaller svg file
  1516. randomIDs = np.random.choice(range(len(X_delay)), 2000)
  1518. # figure
  1519. f,ax=plt.subplots(figsize=(2.1,1.95))
  1521. X_subsample = X_delay[randomIDs]
  1522. ax.scatter(X_subsample[:,0], X_subsample[:,1], alpha=0.2, color='#0B506F', marker='.')
  1523. # fit 250 GMMs on full data with differnt start points to check stability
  1524. gm = GaussianMixture(n_components=chosen_components)
  1525. for i in range(250): # repeat for random starting points
  1526. gm = GaussianMixture(n_components=chosen_components, covariance_type='full', max_iter=500).fit(X_delay)#, covariance_type='spherical'
  1527. # highlight the found classes
  1528. dual_class = np.argmax(np.sum(gm.means_, axis=1))
  1529. for m_id,m in enumerate(gm.means_):
  1530. #if m_id == dual_class :
  1531. # ax.scatter(m[0], m[1], color='r', marker='x', s=50, alpha=0.1)
  1532. #else:
  1533. ax.scatter(m[0], m[1], color='#333333', marker='x', s=50, alpha=0.1)
  1534. # plot diagonal for showing dual reactivation line
  1535. x0,x1 = ax.get_xlim()
  1536. y0,y1 = ax.get_ylim()
  1537. ax.plot(np.linspace(0,np.max([x1, y1]),10), np.linspace(0,np.max([x1, y1]),10),color='#333333', dashes=[1,1])
  1538. ax.set_xlabel('left delay strength')
  1539. ax.set_ylabel('right delay strength')
  1540. sns.despine()
  1541. plt.tight_layout()
  1542. #plt.savefig('./GMM4FitsSingleVSDualReact_DELAY.svg')
  1543. plt.show()
  1545. # %% [markdown]
  1546. # # Fig. 6d: Behavioral precision based on memory in either or both hemispheres
  1548. # %%
  1549. #colors_ic={'ipsi':'#EF8354', 'contra':'#2E9FDC'}
  1551. # PREDICTONS
  1552. contraPred = [df_pred_noBorder.pred_complex_delay_left[idx] if df_pred_noBorder.leftHemi[idx]==CONTRA\
  1553. else df_pred_noBorder.pred_complex_delay_right[idx] if df_pred_noBorder.rightHemi[idx]==CONTRA\
  1554. else [np.nan for i in df_pred_noBorder.pred_complex_delay_left[idx]] for idx in df_pred_noBorder.index]
  1555. ipsiPred = [df_pred_noBorder.pred_complex_delay_left[idx] if df_pred_noBorder.leftHemi[idx]==IPSI\
  1556. else df_pred_noBorder.pred_complex_delay_right[idx] if df_pred_noBorder.rightHemi[idx]==IPSI\
  1557. else [np.nan for i in df_pred_noBorder.pred_complex_delay_left[idx]] for idx in df_pred_noBorder.index]
  1560. df_pred_noBorder['contraDecoderPred'] = contraPred
  1561. df_pred_noBorder['ipsiDecoderPred'] = ipsiPred
  1562. df_pred_noBorder['combinedDecoderPred'] = [np.mean([contraPred[i], ipsiPred[i]], axis=0)\
  1563. for i in range(len(contraPred))]
  1565. # ERRORS
  1566. contraErr = [df_pred_noBorder.prederror_ortho_left[idx] if df_pred_noBorder.leftHemi[idx]==CONTRA\
  1567. else df_pred_noBorder.prederror_ortho_right[idx] if df_pred_noBorder.rightHemi[idx]==CONTRA\
  1568. else [np.nan for i in df_pred_noBorder.prederror_ortho_left[idx]] for idx in df_pred_noBorder.index]
  1569. ipsiErr = [df_pred_noBorder.prederror_ortho_left[idx] if df_pred_noBorder.leftHemi[idx]==IPSI\
  1570. else df_pred_noBorder.prederror_ortho_right[idx] if df_pred_noBorder.rightHemi[idx]==IPSI\
  1571. else [np.nan for i in df_pred_noBorder.prederror_ortho_left[idx]] for idx in df_pred_noBorder.index]
  1574. df_pred_noBorder['contraDecoderErr'] = contraErr
  1575. df_pred_noBorder['ipsiDecoderErr'] = ipsiErr
  1576. df_pred_noBorder['combinedDecoderErr'] = [circdist(np.angle(df_pred_noBorder['combinedDecoderPred'][i]),\
  1577. np.angle(df_pred_noBorder.targ_prev_xy[i]))\
  1578. for i in df_pred_noBorder.index]
  1580. # %%
  1581. precision = {'none':[], 'ipsi':[], 'contra':[], 'both':[]}
  1582. precisionerrors = {'none':[], 'ipsi':[], 'contra':[], 'both':[]}
  1583. binspace = np.linspace(-np.pi/10, np.pi/10, 50)
  1584. precision_mono = {m: {'none':[], 'ipsi':[], 'contra':[], 'both':[]} for m in monkeys}
  1585. precision_orig = []
  1586. labels, label_name, label_errors = [], [], []
  1587. cc=[]
  1588. h=0
  1589. for monkey in monkeys:
  1590. df_mono = df_pred_noBorder.loc[df_pred_noBorder.monkey==monkey]
  1591. for s, session in enumerate(np.unique(df_mono.session_continuous)): # for each session
  1592. df_sess = df_mono.loc[(df_mono.session_continuous==session)].copy().reset_index(drop=True)
  1593. borders = df_sess.borders_full[0]
  1595. # define end of delay
  1596. baseline = range(df_sess.borders_full[0][0], df_sess.borders_full[0][3])
  1597. latedelay = range(df_sess.borders_full[0][6]-3, df_sess.borders_full[0][6])
  1599. # get decoder prediction error at the end of the delay
  1600. contra_strength = [np.mean(cP[latedelay])\
  1601. for cP in df_sess['contraDecoderErr']]#- np.mean(np.abs(cP[baseline]))
  1602. ipsi_strength = [np.mean(iP[latedelay])\
  1603. for iP in df_sess['ipsiDecoderErr']]# - np.mean(np.abs(iP[baseline]))
  1604. errors = np.rad2deg(df_sess.behav_error_prev.values)
  1605. err_std = np.std(errors)
  1606. precision_orig.append(err_std)
  1608. # split into high vs low memory
  1609. cut_perc = 20
  1610. contra_cut, ipsi_cut = np.percentile(np.abs(contra_strength), cut_perc), np.percentile(np.abs(ipsi_strength), cut_perc)
  1611. cc.append(contra_cut)
  1612. # good memory means precision errors < cut_perc
  1613. none_trials = np.where((np.abs(contra_strength) >= contra_cut) & (np.abs(ipsi_strength) >= ipsi_cut))[0]
  1614. ipsi_trials = np.where((np.abs(contra_strength) >= contra_cut) & (np.abs(ipsi_strength) < ipsi_cut))[0]
  1615. contra_trials = np.where((np.abs(contra_strength) < contra_cut) & (np.abs(ipsi_strength) >= ipsi_cut))[0]
  1616. both_trials = np.where((np.abs(contra_strength) < contra_cut) & (np.abs(ipsi_strength) < ipsi_cut))[0]
  1619. # save
  1620. precision['none'].append(np.std(errors[none_trials]))
  1621. precision['ipsi'].append(np.std(errors[ipsi_trials]))
  1622. precision['contra'].append(np.std(errors[contra_trials]))
  1623. precision['both'].append(np.std(errors[both_trials]))
  1625. precisionerrors['none'].append(errors[none_trials])
  1626. precisionerrors['ipsi'].append(errors[ipsi_trials])
  1627. precisionerrors['contra'].append(errors[contra_trials])
  1628. precisionerrors['both'].append(errors[both_trials])
  1630. precision_mono[monkey]['none'].append(np.std(errors[none_trials]))
  1631. precision_mono[monkey]['ipsi'].append(np.std(errors[ipsi_trials]))
  1632. precision_mono[monkey]['contra'].append(np.std(errors[contra_trials]))
  1633. precision_mono[monkey]['both'].append(np.std(errors[both_trials]))
  1635. # labels for each trial
  1636. sess_labels = np.empty(len(errors), dtype=object)
  1637. sess_labels[none_trials] = 'none'
  1638. sess_labels[ipsi_trials] = 'ipsi'
  1639. sess_labels[contra_trials] = 'contra'
  1640. sess_labels[both_trials] = 'both'
  1642. assert len(errors) == len(sess_labels)
  1643. labels.append(sess_labels)
  1644. label_name.append([session for _ in sess_labels])
  1645. label_errors.append(errors)
  1646. h+=1
  1648. # %%
  1649. from scipy import stats
  1651. def print_stats(a=[], b = [], ttest=ttest_rel, text=''):
  1652. test, pval = ttest(a, b)
  1653. degf = len(a) - 1
  1654. mean_x = np.mean(a - np.array(b))
  1655. se = stats.sem(a - np.array(b)) # standard error of the mean
  1656. ci_low, ci_high = t.interval(0.95, degf, loc=mean_x, scale=se)
  1657. print(text+', '+str(ttest)+': t('+str(degf)+') = '+str(np.round(test, 2))+', p = '+"{:.2e}".format(pval)+\
  1658. ", mean = "+str(np.round(mean_x, 2))+\
  1659. ", 95% CI ["+str(np.round(ci_low,2))+","+str(np.round(ci_high, 2))+"]")
  1660. return
  1662. f, ax = plt.subplots(figsize=(1.7,1.3))
  1663. ax.yaxis.set_major_formatter(ticker.FormatStrFormatter('%0.1f'))
  1664. colors_mono = {'Sa': 'darkred', 'Pe': 'lightsalmon', 'Wa':'k'}
  1665. for monkey in monkeys:
  1666. sesss = len(precision_mono[monkey]['none'])
  1667. precisions = [precision_mono[monkey]['none'], precision_mono[monkey]['ipsi'],\
  1668. precision_mono[monkey]['contra'], precision_mono[monkey]['both']]
  1670. #plt.plot([np.zeros((sesss)), np.ones((sesss)), 2*np.ones((sesss)), 3*np.ones((sesss))],\
  1671. # precisions, color=colors_mono[monkey], alpha=0.3)
  1672. plt.axhline(np.mean(precision_orig), color='k', lw=0.5)
  1673. plt.errorbar(0,np.mean(precision['none']), yerr=sem(precision['none']), color='k', marker='o')
  1674. plt.errorbar(1,np.mean(precision['ipsi']), yerr=sem(precision['ipsi']), color=colors['Across'], marker='o')
  1675. plt.errorbar(2,np.mean(precision['contra']), yerr=sem(precision['contra']), color=colors['Across'], marker='o')
  1676. plt.errorbar(3,np.mean(precision['both']), yerr=sem(precision['both']), color=colors['Reactivation'], marker='o')
  1677. plt.ylabel('error std. (°)')
  1678. plt.xticks([0,1,2,3], ['none', 'ipsi', 'contra', 'dual'])
  1679. plt.xlabel('hemispheric memories')
  1681. # p-values
  1682. y0, y1 = ax.get_ylim()
  1683. for l, label in enumerate(['none', 'ipsi', 'contra', 'both']):
  1684. stat = ['***' if ttest_rel(precision[label], precision_orig)[1] < 0.005\
  1685. else '**' if ttest_rel(precision[label], precision_orig)[1] < 0.01\
  1686. else '*' if ttest_rel(precision[label], precision_orig)[1] < 0.05 else 'n.s.'][0]
  1687. ax.annotate(stat,(l, y1), ha='center', va='bottom', fontsize=10, color='#333333')
  1688. ax.plot([2,2,3,3],\
  1689. [y1+ 1*y1/8,y1 + 1*y1/5, y1 + 1*y1/5, y1+ 1*y1/8], linewidth=1, color='#333333')
  1690. ax.annotate('p = '+str(np.round(ttest_rel(precision['both'], precision['contra'])[1],3)),\
  1691. (2.5, y1+y1/5), ha='center', va='bottom', fontsize=9, color='#333333')
  1693. # PRINT STATS
  1695. # BOTH
  1696. print_stats(a=precision['both'], b = precision_orig, ttest=ttest_rel, text='Both-Orig')
  1697. # IPSI
  1698. print_stats(a=precision['ipsi'], b = precision_orig, ttest=ttest_rel, text='Ipsi-Orig')
  1699. # CONTRA
  1700. print_stats(a=precision['contra'], b = precision_orig, ttest=ttest_rel, text='Contra-Orig')
  1701. # BOTH - CONTRA
  1702. print_stats(a=precision['both'], b = precision['contra'], ttest=ttest_rel, text='Both-Contra')
  1703. # NONE
  1704. print_stats(a=precision['none'], b = precision_orig, ttest=ttest_rel, text='None-Orig')
  1705. # BOTH - NONE
  1706. print_stats(a=precision['both'], b = precision['none'], ttest=ttest_rel, text='Both-None')
  1709. plt.tight_layout()
  1710. #plt.savefig(DATAPATH+'Figures/Paper/Neurons/BumpDrift/DualMemoryImprovement'+'.svg')
  1712. # %% [markdown]
  1713. # ### Supplement 9: Different cut-offs
  1715. # %%
  1716. percentages = np.arange(20, 80, 10)
  1717. precision = {'none':[[] for p in percentages], 'ipsi':[[] for p in percentages],\
  1718. 'contra':[[] for p in percentages], 'both':[[] for p in percentages]}
  1719. precisionerrors = {'none':[[] for p in percentages], 'ipsi':[[] for p in percentages],\
  1720. 'contra':[[] for p in percentages], 'both':[[] for p in percentages]}
  1721. binspace = np.linspace(-np.pi/10, np.pi/10, 50)
  1722. precision_orig = []
  1723. labels, label_name, label_errors = [], [], []
  1724. cc=[]
  1725. h=0
  1726. for monkey in monkeys:
  1727. df_mono = df_pred_noBorder.loc[df_pred_noBorder.monkey==monkey]
  1728. for s, session in enumerate(np.unique(df_mono.session_continuous)): # for each session
  1729. df_sess = df_mono.loc[(df_mono.session_continuous==session)].copy().reset_index(drop=True)
  1730. borders = df_sess.borders_full[0]
  1732. # define end of delay
  1733. baseline = range(df_sess.borders_full[0][0], df_sess.borders_full[0][3])
  1734. latedelay = range(df_sess.borders_full[0][6]-3, df_sess.borders_full[0][6])
  1736. # get decoder prediction error at the end of the delay
  1737. contra_strength = [np.mean(cP[latedelay])\
  1738. for cP in df_sess['contraDecoderErr']]#- np.mean(np.abs(cP[baseline]))
  1739. ipsi_strength = [np.mean(iP[latedelay])\
  1740. for iP in df_sess['ipsiDecoderErr']]# - np.mean(np.abs(iP[baseline]))
  1741. errors = np.rad2deg(df_sess.behav_error_prev.values)
  1742. err_std = np.std(errors)
  1743. precision_orig.append(err_std)
  1745. # split into high vs low memory
  1746. for p, cut_perc in enumerate(percentages):
  1747. contra_cut, ipsi_cut = np.percentile(np.abs(contra_strength), cut_perc), np.percentile(np.abs(ipsi_strength), cut_perc)
  1748. cc.append(contra_cut)
  1749. # good memory means precision errors < cut_perc
  1750. none_trials = np.where((np.abs(contra_strength) >= contra_cut) & (np.abs(ipsi_strength) >= ipsi_cut))[0]
  1751. ipsi_trials = np.where((np.abs(contra_strength) >= contra_cut) & (np.abs(ipsi_strength) < ipsi_cut))[0]
  1752. contra_trials = np.where((np.abs(contra_strength) < contra_cut) & (np.abs(ipsi_strength) >= ipsi_cut))[0]
  1753. both_trials = np.where((np.abs(contra_strength) < contra_cut) & (np.abs(ipsi_strength) < ipsi_cut))[0]
  1755. # save
  1756. precision['none'][p].append(np.std(errors[none_trials])-err_std)
  1757. precision['ipsi'][p].append(np.std(errors[ipsi_trials])-err_std)
  1758. precision['contra'][p].append(np.std(errors[contra_trials])-err_std)
  1759. precision['both'][p].append(np.std(errors[both_trials])-err_std)
  1761. precisionerrors['none'][p].append(errors[none_trials])
  1762. precisionerrors['ipsi'][p].append(errors[ipsi_trials])
  1763. precisionerrors['contra'][p].append(errors[contra_trials])
  1764. precisionerrors['both'][p].append(errors[both_trials])
  1766. # labels for each trial
  1767. sess_labels = np.empty(len(errors), dtype=object)
  1768. sess_labels[none_trials] = 'none'
  1769. sess_labels[ipsi_trials] = 'ipsi'
  1770. sess_labels[contra_trials] = 'contra'
  1771. sess_labels[both_trials] = 'both'
  1773. assert len(errors) == len(sess_labels)
  1774. labels.append(sess_labels)
  1775. label_name.append([session for _ in sess_labels])
  1776. label_errors.append(errors)
  1777. h+=1
  1779. cmap_ACC = matplotlib.cm.get_cmap('Greys')
  1780. colors_perc = [cmap_ACC(0.3+i/(len(percentages)+1)) for i in range(len(percentages))]
  1782. f, ax = plt.subplots(figsize=(1.7,2))
  1783. plt.axhline(0, color='k', lw=0.5)
  1784. for p, cut_perc in enumerate(percentages):
  1785. precisions = [precision['none'][p], precision['ipsi'][p],\
  1786. precision['contra'][p], precision['both'][p]]
  1787. plt.plot([0,1,2,3],np.mean(precisions, axis=1), color=colors_perc[p],\
  1788. marker='o', linestyle='None', label=str(cut_perc)+'%')
  1789. plt.ylabel('$\Delta$error std. (°)')
  1790. plt.xticks([0,1,2,3], ['none', 'ipsi', 'contra', 'same'])
  1791. plt.xlabel('hemispheric memories')
  1792. plt.legend(fontsize=5, frameon=True)
  1793. plt.tight_layout()
  1794. #plt.savefig(DATAPATH+'Figures/Paper/Neurons/BumpDrift/SUPPLEMENT_DualMemoryImprovement'+'.svg')
  1796. # %%

Figure5.ipynb at commit bd7d532, no license · at the source

Overview

Authors: Melanie Tschiersch1,2, Akash Umakantha3,4, Ryan C Williamson3,4, Matthew A Smith3,4,5, Joao Barbosa6,7,8, Albert Compte1,9
ORCID iDs: Melanie Tschiersch, Matthew A Smith, Joao Barbosa, Albert Compte
  1. Institut d’Investigacions Biomèdiques August Pi i Sunyer (IDIBAPS), Barcelona, Spain
  2. Programa de doctorat en Biomedicina, Universitat de Barcelona (UB), Barcelona, Spain
  3. Center for the Neural Basis of Cognition, Carnegie Mellon University & University of Pittsburgh, Pittsburgh, PA USA
  4. Carnegie Mellon University Neuroscience Institute, Pittsburgh, PA USA
  5. Carnegie Mellon University Department of Biomedical Engineering, Pittsburgh, PA USA
  6. Laboratoire de Neurosciences Cognitives et Computationnelles, INSERM U960, Ecole Normale Superieure - PSL Research University, Paris, France
  7. Cognitive Neuroimaging Unit, INSERM, CEA, CNRS, Université Paris-Saclay, NeuroSpin center, Gif/Yvette, France
  8. Institut de neuromodulation, GHU Paris, psychiatrie et neurosciences, centre hospitalier Sainte-Anne, pôle hospitalo-universitaire 15, Université Paris Cité, Paris, France
  9. Institut d’Investigacions Biomèdiques de Barcelona (IIBB), CSIC, Barcelona, Spain
Institutions: Consorci Institut D'Investigacions Biomediques August Pi I Sunyer (Spain); Universitat de Barcelona (Spain); University of Pittsburgh (United States); Center for the Neural Basis of Cognition (United States); Carnegie Mellon University (United States); Inserm (France); Université Paris Cité (France); Commissariat à l'Énergie Atomique et aux Énergies Alternatives (France); Université Paris Sciences et Lettres (France); Université Paris-Saclay (France); École Normale Supérieure - PSL (France); Centre Hospitalier Sainte-Anne (France); CEA Paris-Saclay (France); Laboratoire de Neurosciences Cognitives et Computationnelles; Consejo Superior de Investigaciones Científicas (Spain); Institut d'Investigacions Biomèdiques de Barcelona (Spain)
Journal: Nature communications, volume 17, issue 1, article 8858
Dates: received 4 August 2025; accepted 7 July 2026; published online 20 July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41467-026-75705-2 · PMID 42477323 · PMCID PMC13500473 · OpenAlex W7169787509
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: computational modeling (no new data) (modality), human (organism), non-human primate (organism), cognitive (subfield)
Methods: Connectivity, Statistics, Machine learning, Preprocessing, Single-unit activity, calcium imaging, Smoothing, state filtering, decompositions, Physiology & signal measures
Keywords: Network models, Short-term memory, Working memory, Cortex, Neural circuits
MeSH: Functional Laterality*, Memory, Short-Term*, Prefrontal Cortex*, Animals, Humans, Macaca mulatta, Male (* major topic)
Topic: Neural and Behavioral Psychology Studies (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: Generalitat de Catalunya (2021SGR01522, 2017SGR01565); Ministry of Economy and Competitiveness | Agencia Estatal de Investigación (Spanish Agencia Estatal de Investigación) (PID2021-125453OB-I00, PCI2025-167127-2, FPI Program); U.S. Department of Health &amp; Human Services | National Institutes of Health (R01 NS147766, R01 EB026953, R01 EY022928, R01 MH118929); Ministry of Economy and Competitiveness | Instituto de Salud Carlos III (AC20/00071); National Science Foundation (NCS BCS 1734916/1954107); NEI NIH HHS (R01 EY022928); NIMH NIH HHS (R01 MH118929); U.S. Department of Health & Human Services | National Institutes of Health (NIH) (R01 MH118929, R01 EY022928, R01 EB026953, R01 NS147766); NINDS NIH HHS (R01 NS147766); Simons Foundation; NIBIB NIH HHS (R01 EB026953)
Citations: cited by 1 paper (Europe PMC); 111 references in the paper

Abstract

The prefrontal hemispheres must coordinate dynamically to maintain a unified representation of visual space. Recently, two opposing theories using distinct storage strategies have been proposed: A high-capacity specialized architecture, where each hemisphere governs contralateral behavior, and a fail-safe redundant one, where both hemispheres jointly guide behavior across the visual space. Here, we analyzed simultaneous bilateral prefrontal cortex recordings from three male macaque monkeys performing a visuo-spatial working memory task. Both hemispheres equally predicted behavioral imprecision, decoding errors were weakly correlated between hemispheres, and serial dependence remained local within hemispheres, suggesting a redundant, weakly coupled organization. Attractor network simulations showed that redundancy improved precision when task demands were below memory capacity, while weak interhemispheric coupling increased capacity in more demanding tasks by allowing hemispheric specialization. These predicted patterns were validated in human and monkey data, reconciling previous findings and revealing a versatile interhemispheric architecture that adapts to varying cognitive demands.

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

Zenodo 20513995

License: CC-BY-4.0
State: the link answers, verified on 27 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 (33 files), SciPy (32 files), Matplotlib (29 files), pandas (28 files), seaborn (19 files), scikit-learn (17 files), statsmodels (16 files), Pingouin (5 files), Brian 2 (4 files), h5py (4 files)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
35 files
At the source:

melanietschiersch/redundantprefrontalhemispheres

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: bd7d5329176b24441eb457aecef734d02967ee04, 2 June 2026
Languages: Python (25), Jupyter (9)
Size: 40 files, 34 scripts
Software Heritage: not archived
Found in: the Zenodo archive record
Holds: README, environment (environment.yml), 9 notebooks
Not found: license file, CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (33 files), SciPy (32 files), Matplotlib (29 files), pandas (28 files), seaborn (19 files), scikit-learn (17 files), statsmodels (16 files), Pingouin (5 files), Brian 2 (4 files), h5py (4 files)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
35 files
At the source: github.com/MelanieTschiersch/RedundantPrefrontalHemispheres/tree/v1.0.0

Code availability

The code accompanying the paper, specifying analyses and models is available on Github111.

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;
  • 68 scripts, each with its path and the digest of its content;
  • 16 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 non-human primate data54 analyzed in this study are available at the KiltHub database from Carnegie Mellon University55. The human data90,91 analyzed in this study have been retrieved from the OSF database under accession code krv7g and 67tn3 (retrieved as OSF data sets osf.io/krv7g, and osf.io/67tn3).

The code accompanying the paper, specifying analyses and models is available on Github111.

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

Versions

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

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 6 authors, 5 keywords, 7 MeSH terms, 11 funders, 107 references.

Cite

This paper

Tschiersch, M., Umakantha, A., Williamson, R. C., Smith, M. A., Barbosa, J., & Compte, A. (2026). Redundant prefrontal hemispheres adapt storage strategy to working memory demands. Nature communications, 17(1), 8858. https://doi.org/10.1038/s41467-026-75705-2

BibTeX

@article{tschiersch2026redundant,
author = {Tschiersch, Melanie and Umakantha, Akash and Williamson, Ryan C and Smith, Matthew A and Barbosa, Joao and Compte, Albert},
title = {{Redundant prefrontal hemispheres adapt storage strategy to working memory demands}},
journal = {Nature communications},
year = {2026},
month = jul,
volume = {17},
number = {1},
pages = {8858},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-75705-2},
url = {https://doi.org/10.1038/s41467-026-75705-2},
pmid = {42477323},
pmcid = {PMC13500473}
}

RIS

TY - JOUR
AU - Tschiersch, Melanie
AU - Umakantha, Akash
AU - Williamson, Ryan C
AU - Smith, Matthew A
AU - Barbosa, Joao
AU - Compte, Albert
TI - Redundant prefrontal hemispheres adapt storage strategy to working memory demands
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/07/20
VL - 17
IS - 1
SP - 8858
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-75705-2
UR - https://doi.org/10.1038/s41467-026-75705-2
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-75705-2",
"type": "article-journal",
"title": "Redundant prefrontal hemispheres adapt storage strategy to working memory demands",
"container-title": "Nature communications",
"author": [
{
"family": "Tschiersch",
"given": "Melanie"
},
{
"family": "Umakantha",
"given": "Akash"
},
{
"family": "Williamson",
"given": "Ryan C"
},
{
"family": "Smith",
"given": "Matthew A"
},
{
"family": "Barbosa",
"given": "Joao"
},
{
"family": "Compte",
"given": "Albert"
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "8858",
"DOI": "10.1038/s41467-026-75705-2",
"PMID": "42477323",
"PMCID": "PMC13500473",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-75705-2",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
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.1038/s42003-026-10071-9 [code]
Alpha phase coding supports feature binding during working memory maintenance.
Journal: Communications biology
In common: Brian 2, scikit-learn, SciPy, 2 other tools, cognitive, 6 references, 2 authors
[2] doi:10.1371/journal.pbio.3003333
Attractive serial dependence arises during decision-making
Journal: n/a
In common: cognitive, 14 references
[3] 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, pandas, SciPy, 2 other tools, non-human primate, cognitive, 6 references, author Matthew A Smith
[4] doi:10.1186/s12915-026-02515-9 [code]
Different modality-specific mechanisms mediate serial dependence effects in visual and auditory perception.
Journal: BMC biology
In common: cognitive, 9 references
[5] doi:10.1038/s41467-026-75704-3 [code]
A minimal model of working memory in neural systems and neuromorphic circuits.
Journal: Nature communications
In common: Brian 2, SciPy, Matplotlib, 1 other tool, computational modeling (no new data), 4 references
[6] doi:10.1016/j.isci.2026.117492 [code]
Neural subspace reorganization reflects value-based decision-making.
Journal: iScience
In common: statsmodels, scikit-learn, pandas, 3 other tools, cognitive, 4 references
[7] doi:10.1038/s41562-026-02414-7 [code]
Optimized feature gains explain and predict successes and failures of human selective listening.
Journal: Nature human behaviour
In common: Pingouin, h5py, statsmodels, 6 other tools, cognitive, 1 reference
[8] doi:10.7554/elife.99278 [code]
Brain-wide arousal signals are segregated from movement planning in the superior colliculus of the macaque.
Journal: eLife
In common: h5py, NumPy, non-human primate, 2 references, author Matthew A Smith
[9] doi:10.1038/s41593-026-02314-z [code]
Low-dimensional population dynamics in the brainstem gate REM sleep.
Journal: Nature neuroscience
In common: Pingouin, h5py, statsmodels, 6 other tools, 1 reference
[10] doi:10.1038/s41467-026-74823-1 [code]
Cerebellar activity is triggered by reach endpoint during learning of a complex locomotor task.
Journal: Nature communications
In common: Pingouin, h5py, statsmodels, 6 other tools, 1 reference

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.