OSCR

Sensorimotor remapping drives task specialization in prefrontal cortex.

Code ↔ Paper

21 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 21 matches
  1. [1] § Results › mPFC neurons are most involved in the representation of variables in their preferred task ↔ Response-switch/Pseudo_pop/regression.ipynb, lines 762–793 · score 0.70 · post stimulus activity, pre stimulus, task variables, task selective, motor, trained
  2. [2] § Methods › Ridge regression › cvR2 scores and statistics ↔ Response-switch/Individuals/regression.ipynb, lines 68–212 · score 0.68 · cvR2, cross validation, variance explained, permutation, pipeline, splitting
  3. [3] § Methods › Ridge regression › cvR2 scores and statistics ↔ Response-switch/Pseudo_pop/regression.ipynb, lines 102–243 · score 0.68 · cvR2, cross validation, variance explained, permutation, pipeline, splitting
  4. [4] § Results › Population modularity can be explained by an interaction between task signal and global inhibition ↔ Simulation/model_response_switch.ipynb, lines 207–285 · score 0.65 · activation function, projection weight, standard deviations, task signal, networks, modulation
  5. [5] § Methods › Subpopulations and Gaussian mixture › Gaussian mixture ↔ Response-switch/Pseudo_pop/gaussian_mixture.ipynb, lines 128–148 · score 0.65 · BayesianGaussianMixture, cross validated, split, GM, components, scores
  6. [6] § Methods › Subpopulations and Gaussian mixture › Gaussian mixture ↔ Rule-switch/Pseudo_pop/gaussian_mixture_rs.ipynb, lines 134–154 · score 0.65 · BayesianGaussianMixture, cross validated, split, GM, components, scores
  7. [7] § Results › Population modularity can be explained by an interaction between task signal and global inhibition ↔ Simulation/model_rule_switch.ipynb, lines 192–277 · score 0.65 · activation function, projection weight, standard deviations, task signal, networks, Rule switch
  8. [8] § Results › Task selectivity delineates three neural subpopulations with distinct functional roles ↔ Response-switch/Pseudo_pop/gaussian_mixture.ipynb, lines 654–689 · score 0.63 · population LR, Gaussian Mixture, Population GNG, NoGo, GM
  9. [9] § Methods › Analyses of neural selectivity ↔ Response-switch/Pseudo_pop/regression.ipynb, lines 762–793 · score 0.62 · pre stimulus, post stimulus, task selectivity, interaction, Lick, trained
  10. [10] § Methods › Dimension reduction ↔ Response-switch/Pseudo_pop/PCA.ipynb, lines 226–279 · score 0.62 · cross validated, variance explained, favored, cvPCA, covariance, trained
  11. [11] § Methods › Recurrent neural network modeling › Model performance on each dual-task paradigm ↔ Simulation/model_response_switch.ipynb, lines 1188–1208 · score 0.62 · direction readout, Lick readout, perceptron, spout, simulated, response switch
  12. [12] § Results › Population modularity can be explained by an interaction between task signal and global inhibition ↔ Simulation/model_response_switch.ipynb, lines 207–285 · score 0.60 · activation functions, projections weights, task signal, networks, modulated, model
  13. [13] § Results › Task selectivity delineates three neural subpopulations with distinct functional roles ↔ Response-switch/Pseudo_pop/gaussian_mixture.ipynb, lines 809–921 · score 0.60 · Support Vector Classifiers, classifiers trained, Gaussian mixture, SVCs, decoding, GNG
  14. [14] § Results › Population modularity can be explained by an interaction between task signal and global inhibition ↔ Simulation/model_rule_switch.ipynb, lines 192–277 · score 0.59 · activation functions, projections weights, task signal, networks, modulated, model
  15. [15] § Methods › Ridge regression › Regressor matrix construction ↔ Response-switch/Individuals/regression.ipynb, lines 68–212 · score 0.58 · standard scaler, cross validation, fold, predict, R2, regressor
  16. [16] § Methods › Ridge regression › Regressor matrix construction ↔ Response-switch/Pseudo_pop/regression.ipynb, lines 102–243 · score 0.58 · standard scaler, cross validation, fold, predict, R2, regressor
  17. [17] § Methods › Ridge regression › Regressor matrix construction ↔ Response-switch/Pseudo_pop/regression.ipynb, lines 433–467 · score 0.55 · Akaike Information Criterion, AIC, split, regressors, scores, model
  18. [18] § Methods › Pseudo-population ↔ Response-switch/Pseudo_pop/gaussian_mixture.ipynb, lines 217–245 · score 0.53 · conservative criterion, Gaussian mixture, pseudo, Response switch, neurons, population
  19. [19] § Methods › Analyses of neural selectivity ↔ Response-switch/Pseudo_pop/gaussian_mixture.ipynb, lines 548–569 · score 0.52 · Lick LR, absolute selectivity, variables, regression
  20. [20] § Results › Modularity of mPFC population decreases in a paradigm with less sensorimotor remapping ↔ Rule-switch/Pseudo_pop/gaussian_mixture_rs.ipynb, lines 579–668 · score 0.51 · Support vector classifiers, Gaussian mixture, SVCs, decoding, splits, trained
  21. [21] § Results › Modularity of mPFC population decreases in a paradigm with less sensorimotor remapping ↔ Response-switch/Pseudo_pop/gaussian_mixture.ipynb, lines 809–921 · score 0.51 · Support vector classifiers, Gaussian mixture, SVCs, decoding, splits, 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,366 lines · 45 KB · no license · 6 matches

  1. # %%
  2. import numpy as np
  3. import pandas as pd
  4. import pickle
  5. import matplotlib.pyplot as plt
  6. from tqdm.auto import tqdm
  7. import matplotlib as mpl
  8. from matplotlib.patches import Ellipse
  9. from matplotlib.lines import Line2D
  10. import matplotlib.transforms as transforms
  11. import seaborn as sns
  12. from scipy.interpolate import make_interp_spline, BSpline, splrep, splev
  13. from scipy.stats import mannwhitneyu, wilcoxon
  15. plt.style.use('seaborn-v0_8-ticks')
  16. mpl_config = pd.read_csv('../../mpl_config.csv').to_dict(orient='records')[0]
  17. mpl.rcParams.update(mpl_config)
  19. from pathlib import Path
  20. import sys
  22. ROOT = Path.cwd().parents[1] # go up 2 levels
  23. sys.path.insert(0, str(ROOT))
  25. from export_source import export_plot_data
  27. # %%
  28. def standardize(M) :
  29. """ Standardization by column for 2D matrices """
  31. M_copy = np.copy(M)
  32. M_copy = M_copy.astype(float)
  33. for j in range(np.shape(M)[1]) :
  34. if np.std(M[:,j]) != 0 :
  35. M_copy[:,j] = (M[:,j].astype(float) - float(np.mean(M[:,j])))/float(np.std(M[:,j]))
  36. else :
  37. M_copy[:,j] = M[:,j].astype(float) - float(np.mean(M[:,j]))
  39. return M_copy
  42. def confidence_ellipse(x, y, ax, n_std=2.0, edgecolor = 'black',facecolor='none', **kwargs):
  43. """
  44. Create a plot of the covariance confidence ellipse of `x` and `y`
  45. """
  47. if x.size != y.size:
  48. raise ValueError("x and y must be the same size")
  50. cov = np.cov(x, y)
  51. pearson = cov[0, 1]/np.sqrt(cov[0, 0] * cov[1, 1])
  52. # Using a special case to obtain the eigenvalues of this
  53. # two-dimensionl dataset.
  54. ell_radius_x = np.sqrt(1 + pearson)
  55. ell_radius_y = np.sqrt(1 - pearson)
  56. ellipse = Ellipse((0, 0),
  57. width=ell_radius_x * 2,
  58. height=ell_radius_y * 2,
  59. facecolor=facecolor,
  60. edgecolor=edgecolor,
  61. **kwargs)
  63. # Calculating the stdandard deviation of x from
  64. # the squareroot of the variance and multiplying
  65. # with the given number of standard deviations.
  66. scale_x = np.sqrt(cov[0, 0]) * n_std
  67. mean_x = np.mean(x)
  69. # calculating the stdandard deviation of y ...
  70. scale_y = np.sqrt(cov[1, 1]) * n_std
  71. mean_y = np.mean(y)
  73. transf = transforms.Affine2D() \
  74. .rotate_deg(45) \
  75. .scale(scale_x, scale_y) \
  76. .translate(mean_x, mean_y)
  78. ellipse.set_transform(transf + ax.transData)
  79. return ax.add_patch(ellipse)
  81. def smooth(x,y,nb_point) :
  82. """Smooth a (x,y) serie by inserting nb_point with spline interpolation"""
  83. x_smooth = np.linspace(np.min(x), np.max(x), nb_point)
  84. y_smooth = make_interp_spline(x, y, k=3)(x_smooth)
  86. return x_smooth, y_smooth
  88. # %%
  89. ## We load here the selected ridge and PCA models
  91. with open('Models/ridge.pickle', 'rb') as f:
  92. ridge_output = pickle.load(f)
  94. with open('Models/pca_notime.pickle', 'rb') as f:
  95. pca1 = pickle.load(f)
  97. with open('../DATA/Dataframes/df_pseudo.pickle', 'rb') as f:
  98. merged_data = pickle.load(f)
  100. res = ridge_output['Models']
  101. reg_scores = ridge_output['Scores']
  103. # %% [markdown]
  104. # ### Selectivity matrix construction
  106. # %%
  107. plt.figure(figsize=(7,7))
  109. m = res[np.argmax(reg_scores)]
  110. S = m['Ridge'].coef_.T
  112. ax = plt.gca()
  114. plt.scatter(S[0,:],S[5,:],color='grey', s = 4)
  116. ax.spines['left'].set_position('center')
  117. ax.spines['bottom'].set_position('center')
  119. ax.spines['left'].set_linewidth(1.4)
  120. ax.spines['bottom'].set_linewidth(1.4)
  122. ax.spines['right'].set_color('none')
  123. ax.spines['top'].set_color('none')
  125. plt.xticks([])
  126. plt.yticks([])
  128. # %% [markdown]
  129. # ### Components estimation
  131. # %%
  132. from sklearn.model_selection import KFold, GridSearchCV, RepeatedKFold, cross_validate
  133. from sklearn.mixture import GaussianMixture, BayesianGaussianMixture
  134. from tqdm.auto import tqdm
  136. def gm_components_estimation(S_mat, n_repeat=50) :
  137. """Determine the number of components to keep in Gaussian mixture through a 2-fold cross-validation procedure (50 reps)"""
  138. gm_scores = []
  140. rkf = RepeatedKFold(n_splits = 2, n_repeats = n_repeat)
  141. for n_comp in tqdm(range(1,6)) :
  142. gm = BayesianGaussianMixture(n_components=n_comp ,n_init = 5, max_iter=100)
  144. cv_scores = cross_validate(gm,S_mat.T,cv=rkf)
  146. gm_scores.append(cv_scores['test_score'])
  148. return np.array(gm_scores)
  151. # %%
  152. import warnings
  153. warnings.filterwarnings("ignore")
  155. gm_scores_full = gm_components_estimation(S)
  157. diff_gm_scores = np.diff(gm_scores_full,axis=0) ## Computing relative increase in model performance after the addition of each component
  159. # %%
  160. import matplotlib as mpl
  162. plt.figure(figsize=(8,3))
  164. yerr = np.abs(np.mean(diff_gm_scores,axis=1) - np.percentile(diff_gm_scores,(5,95),axis=1))
  166. _,caps,_ = plt.errorbar([2,3,4,5],np.mean(diff_gm_scores,axis=1),yerr=yerr,linewidth=3,capsize=5,elinewidth=2.5,color='#323232')
  169. for cap in caps:
  170. cap.set_color('#323232')
  171. cap.set_markeredgewidth(2)
  173. plt.scatter([2,3,4,5],np.mean(diff_gm_scores,axis=1),marker = 's',s=150, color='#323232')
  174. #plt.scatter([2,3,4,5],np.percentile(np.diff(gm_scores_full,axis=0),5,axis=1),marker='s', s=120,color='#276690')
  176. plt.fill_between([1.9,3,4,5.1],-0.2,0, color='black',alpha=0.05)
  177. plt.axhline(0,color="black",linestyle='--',alpha=0.2)
  179. plt.xlabel('# Subpopulations',fontsize=30)
  180. plt.xlim((1.9,5.1))
  181. plt.xticks([2,3,4,5],fontsize=28)
  183. plt.ylabel('LL increase \n in model fitting',fontsize=30)
  185. #plt.ylim((-0.15,0.4))
  186. plt.ylim((-0.1,0.6))
  188. #plt.yticks([-0.05,0,0.05,0.1])
  189. plt.yticks([0,0.2,0.4,0.6],fontsize=28)
  192. #plt.savefig('Plots/SVG/gm_comp.SVG', dpi = 300,bbox_inches='tight')
  193. #plt.savefig('Plots/PNG/gm_comp.PNG', dpi = 300,bbox_inches='tight')
  195. # %%
  196. # Store data (serialize)
  197. with open('Models/gm_components.pickle', 'wb') as handle:
  198. pickle.dump(np.diff(gm_scores_full,axis=0) , handle, protocol=pickle.HIGHEST_PROTOCOL)
  200. # %%
  201. with open('Models/gm_components.pickle', 'rb') as f:
  202. diff_gm_scores = pickle.load(f)
  204. # %%
  205. export_plot_data('../../Source_Data.xlsx','Figure4c',LL_increase=diff_gm_scores)
  207. # %% [markdown]
  208. # ### GM model
  210. # %%
  211. from sklearn.mixture import GaussianMixture, BayesianGaussianMixture
  213. gm = BayesianGaussianMixture(n_components = 3,n_init = 20, max_iter=100) ## Keeping 3 components (see above)
  214. gm.fit(S.T)
  215. prob = gm.predict_proba(S.T)
  217. # %% [markdown]
  218. # ### Define population
  220. # %%
  221. pop1 = prob[:,0] >= 0.9 #Conservative criterion to associate each neuron with a GM component.
  222. pop2 = prob[:,1] >= 0.9 #A neuron is associated to a population if the weight associated to it is more than 90%
  223. pop3 = prob[:,2] >= 0.9
  225. pops = [pop1,pop2,pop3]
  227. S_pop = [S[:,pop1],S[:,pop2],S[:,pop3]]
  230. mean_ctx_s = [np.mean(S_pop[0][5,:]),np.mean(S_pop[1][5,:]),np.mean(S_pop[2][5,:])]
  231. print(mean_ctx_s)
  233. S_B = S_pop[np.argsort(mean_ctx_s)[2]]
  234. popB = pops[np.argsort(mean_ctx_s)[2]]
  236. S_A = S_pop[np.argsort(mean_ctx_s)[0]]
  237. popA = pops[np.argsort(mean_ctx_s)[0]]
  239. S_0 = S_pop[np.argsort(mean_ctx_s)[1]]
  240. pop0 = pops[np.argsort(mean_ctx_s)[1]]
  242. print(len(S.T))
  243. print(np.sum(popA),np.sum(popA)/len(S.T))
  244. print(np.sum(popB),np.sum(popB)/len(S.T))
  245. print(np.sum(pop0),np.sum(pop0)/len(S.T))
  248. # %%
  249. with open('../DATA/Dataframes/df_full.pickle', 'rb') as f:
  250. data_full = pickle.load(f)
  252. files = ['M12','M13','M14','M15','M16','M17','M18','M19','M20']
  254. neuron_per_animal = []
  255. for file in files :
  256. neuron_per_animal.append(len(data_full[file]['Spike rate'][0]))
  257. animal_slice = np.insert(np.cumsum(neuron_per_animal), 0, 0, axis=0)
  258. popA_per_animal = [np.sum(popA[animal_slice[i]:animal_slice[i+1]]) for i in range(len(animal_slice)-1)]
  259. popB_per_animal = [np.sum(popB[animal_slice[i]:animal_slice[i+1]]) for i in range(len(animal_slice)-1)]
  260. pop0_per_animal = [np.sum(pop0[animal_slice[i]:animal_slice[i+1]]) for i in range(len(animal_slice)-1)]
  262. def plot_composition_bars(list1, list2, list3, labels=("Population GNG", "Population LR", "Population 0")):
  263. list1 = np.array(list1)
  264. list2 = np.array(list2)
  265. list3 = np.array(list3)
  267. n = len(list1)
  268. x = np.arange(n)
  270. plt.figure(figsize=(10, 7))
  272. # Bottom segment
  273. plt.bar(x, list1, label=labels[0],color='#2E548A')
  275. # Middle segment
  276. plt.bar(x, list2, bottom=list1, label=labels[1],color='#E63946')
  278. # Top segment
  279. plt.bar(x, list3, bottom=list1 + list2, label=labels[2],color='#5F8162')
  281. plt.xticks(range(len(files)),files)
  282. plt.ylabel("# Neurons")
  283. plt.ylim(0,400)
  284. plt.legend()
  285. plt.tight_layout()
  287. plot_composition_bars(popA_per_animal,popB_per_animal,pop0_per_animal)
  289. #plt.savefig('Plots/SVG/pop_repartition_resp.SVG', dpi = 300,bbox_inches='tight')
  290. #plt.savefig('Plots/PNG/pop_repartition_resp.PNG', dpi = 300,bbox_inches='tight')
  292. # %%
  293. from scipy.stats import binomtest
  294. from statsmodels.stats.multitest import multipletests
  296. pop_table = np.stack((popA_per_animal,popB_per_animal,pop0_per_animal)).T
  298. def binomial_enrichment_pmatrix(contingency, method = "bonferroni"):
  299. """
  300. Binomial enrichment test for all (sample, category) pairs.
  302. Parameters
  303. ----------
  304. contingency : np.ndarray
  305. 2D array with shape (n_samples, n_categories)
  306. method : str
  307. Multiple testing correction method (default: FDR Benjamini–Hochberg)
  309. Returns
  310. -------
  311. np.ndarray
  312. 2D array of corrected p-values with same shape as contingency
  313. """
  314. contingency = np.asarray(contingency, dtype=int)
  316. T = contingency.sum()
  317. row_totals = contingency.sum(axis=1) # samples
  318. col_totals = contingency.sum(axis=0) # categories
  320. n_samples, n_categories = contingency.shape
  321. pvals = np.ones((n_samples, n_categories))
  323. for i in range(n_samples):
  324. p0 = row_totals[i] / T if T > 0 else 0.0
  326. for j in range(n_categories):
  327. k = contingency[i, j]
  328. n = col_totals[j]
  330. if n > 0 and p0 > 0:
  331. pvals[i, j] = binomtest(
  332. k, int(n), p=0.5, alternative="greater"
  333. ).pvalue
  334. else:
  335. pvals[i, j] = 1.0
  337. # multiple testing correction
  338. pvals_corr = multipletests(
  339. pvals.ravel(), method=method
  340. )[1].reshape(pvals.shape)
  342. return pvals_corr
  344. binomial_enrichment_pmatrix(pop_table)
  346. # %%
  347. def permutation_global_driver_test(
  348. counts,
  349. threshold = 0.5,
  350. n_permutations = 10000,
  351. random_state = None):
  352. """
  353. Global permutation test for whether at least one category
  354. is mainly driven by one sample.
  356. Parameters
  357. ----------
  358. counts : np.ndarray
  359. Contingency table of shape (n_samples, n_categories)
  360. threshold : float
  361. Dominance threshold (default: 0.5)
  362. n_permutations : int
  363. Number of permutations
  364. random_state : int or None
  365. Random seed
  367. Returns
  368. -------
  369. float
  370. Global p-value
  371. """
  372. rng = np.random.default_rng(random_state)
  373. counts = np.asarray(counts, dtype=int)
  375. n_samples, n_categories = counts.shape
  376. col_totals = counts.sum(axis=0)
  377. grand_total = counts.sum()
  379. # global sample proportions (null model)
  380. sample_probs = counts.sum(axis=1) / grand_total
  382. # observed global statistic
  383. with np.errstate(divide="ignore", invalid="ignore"):
  384. observed_shares = counts / col_totals
  385. observed_shares[:, col_totals == 0] = 0.0
  387. T_obs = observed_shares.max()
  389. # early exit: nothing exceeds the threshold
  390. if T_obs <= threshold:
  391. return 1.0
  393. exceed_count = 0
  395. for _ in range(n_permutations):
  396. T_perm = 0.0
  398. for j in range(n_categories):
  399. n_j = col_totals[j]
  400. if n_j == 0:
  401. continue
  403. # multinomial draw under null
  404. perm_counts = rng.multinomial(n_j, sample_probs)
  405. T_perm = max(T_perm, perm_counts.max() / n_j)
  407. # early stopping for speed
  408. if T_perm >= T_obs:
  409. break
  411. if T_perm >= T_obs:
  412. exceed_count += 1
  414. # +1 correction for unbiased p-value
  415. p_value = (exceed_count + 1) / (n_permutations + 1)
  416. return p_value
  419. permutation_global_driver_test(pop_table)
  421. # %% [markdown]
  422. # ### Plots of the cloud of dots
  424. # %%
  425. plt.figure(figsize=(6,6))
  427. plt.xlim((-1.2,1.2))
  428. plt.ylim((-1.2,1.2))
  430. ax = plt.gca()
  432. mode = 'center'
  434. if mode == 'center' :
  436. ax.spines['left'].set_position('center')
  437. ax.spines['bottom'].set_position('center')
  439. ax.spines['left'].set_linewidth(1.4)
  440. ax.spines['bottom'].set_linewidth(1.4)
  442. ax.spines['right'].set_color('none')
  443. ax.spines['top'].set_color('none')
  445. plt.xticks([])
  446. plt.yticks([])
  448. elif mode == 'normal' :
  449. ax.spines[['top','right']].set_visible(False)
  451. plt.yticks([-1,0,1],fontsize=25)
  452. plt.xticks([-1,0,1],fontsize=25)
  453. ax.set_xlabel(r'$\beta_{Category G/NG}$')
  454. ax.set_ylabel(r'$\beta_{Task}$')
  458. plt.scatter(S_0[0,:],S_0[5,:], color = 'green', s = 8, alpha = 0.4,label='Population 0'\
  459. ,edgecolors='#EEEEEE',lw=0.7,antialiased=True, zorder=3)
  460. plt.scatter(S_B[0,:],S_B[5,:], color = 'red', s = 20, alpha = 1,label='Population LR'\
  461. ,edgecolors='#EEEEEE',lw=0.7,antialiased=True, zorder=3)
  462. plt.scatter(S_A[0,:],S_A[5,:], color = 'blue', s = 20, alpha = 1,label='Population GNG'\
  463. ,edgecolors='#EEEEEE',lw=0.7,antialiased=True, zorder=3)
  465. #plt.legend(frameon=False)
  467. confidence_ellipse(S_B[0,:],S_B[5,:],ax, edgecolor='red', n_std=2, alpha = 0.5,linewidth=2)
  468. confidence_ellipse(S_A[0,:],S_A[5,:],ax, edgecolor='blue', n_std=2, alpha = 0.5,linewidth=2)
  471. plt.savefig('Plots/SVG/gm_catGNG_task_' + mode + '.SVG', dpi = 300,bbox_inches='tight')
  472. plt.savefig('Plots/PNG/gm_catGNG_task_' + mode + '.PNG', dpi = 300,bbox_inches='tight')
  474. # %% [markdown]
  475. # #### 3D
  477. # %%
  478. pv.set_jupyter_backend('client')
  479. p = pv.Plotter(window_size = (800,600))
  480. p.set_background('white')
  482. p.subplot(0,0)
  484. #p.set_position(np.mean(S[GNG_trials,:3],axis=0))
  486. cord = [2,3,5]
  488. #labels = dict(ztitle=r'$\beta_{Task}$', xtitle=r'$\beta_{Category GNG}$', ytitle=r'$\beta_{Choice GNG}$')
  490. #----- This block generates the sphere for each trial label and render them in the 3D plot----
  491. radius = 0.01
  493. n_perpop = 400
  494. popGNG_spheres = [pv.Sphere(center = S[:,popA][cord,i], radius = radius) for i in range(n_perpop)]
  495. popLR_spheres = [pv.Sphere(center = S[:,popB][cord,i], radius = radius) for i in range(n_perpop)]
  496. pop0_spheres = [pv.Sphere(center = S[:,pop0][cord,i], radius = radius) for i in range(n_perpop)]
  498. for sphere in popGNG_spheres :
  499. p.add_mesh(sphere, color='blue', opacity = 0.4, lighting=False)
  501. for sphere in popLR_spheres :
  502. p.add_mesh(sphere, color='red', opacity = 0.4, lighting=False)
  504. for sphere in pop0_spheres :
  505. p.add_mesh(sphere, color='green', opacity = 0.4, lighting=False)
  508. torus_GNG = pv.ParametricTorus(ringradius=2*np.std(S[:,popA][cord[:2],:n_perpop]), crosssectionradius=5*10**(-3),center=np.mean(S[:,popA][cord,:n_perpop],axis=1))
  509. p.add_mesh(torus_GNG, color='blue', lighting=False)
  511. torus_LR = pv.ParametricTorus(ringradius=2*np.std(S[:,popB][cord[:2],:n_perpop]), crosssectionradius=5*10**(-3),center=np.mean(S[:,popB][cord,:n_perpop],axis=1))
  512. p.add_mesh(torus_LR, color='red', lighting=False)
  514. torus_0 = pv.ParametricTorus(ringradius=2*np.std(S[:,pop0][cord[:2],:n_perpop]), crosssectionradius=5*10**(-3),center=np.mean(S[:,pop0][cord,:n_perpop],axis=1))
  515. p.add_mesh(torus_0, color='green', lighting=False)
  516. #p.set_position(np.array([40,0,0]))
  518. #p.show_grid(**labels)
  519. p.show()
  522. # %% [markdown]
  523. # ### Violin
  525. # %%
  526. ## Building a dataframe with data to use seaborn module for prettier plots
  528. plot_labels=['Category GNG', 'Choice GNG','Category LR', 'Choice LR', 'Lick LR', 'Task']
  530. s_todf = []
  531. reg_todf = []
  532. cluster_todf = []
  533. for k in range(np.shape(S)[1]) :
  534. s_todf.append(S[:6,k])
  535. reg_todf.append(plot_labels)
  536. if popA[k] :
  537. cluster_todf.append(['A']*(np.shape(S)[0] - 1))
  538. elif popB[k]:
  539. cluster_todf.append(['B']*(np.shape(S)[0] - 1))
  540. else :
  541. cluster_todf.append(['else']*(np.shape(S)[0] - 1))
  543. df_const = {'Selectivity' : np.array(s_todf).flatten(), 'Abs selectivity' : np.array(np.abs(s_todf)).flatten() ,'Regressor' : np.array(reg_todf).flatten(), 'Cluster' : np.array(cluster_todf).flatten()}
  544. plot_df = pd.DataFrame(df_const)
  548. # %%
  549. ## Absolute selectivity across task variables for different populations
  551. plt.figure(figsize=(12,9))
  553. plot_dfAB = plot_df[(plot_df['Cluster'] != 'else')&(plot_df['Regressor'] != 'Task')&(plot_df['Regressor'] != 'Lick LR')]
  555. p = sns.violinplot(data = plot_dfAB, x='Regressor',y='Abs selectivity',hue='Cluster',\
  556. hue_order=['A','B'],alpha=0.7,palette=['blue','red'],cut=0,common_norm=True,\
  557. inner_kws=dict(box_width=7, whis_width=2))
  560. plt.ylabel(r'|$\beta$|',fontsize=40)
  561. plt.xlabel('')
  562. plt.legend([],[], frameon=False)
  563. plt.axvline(x = 1.5, color = 'black', linewidth=1.6)
  564. plt.ylim(0,1)
  565. plt.yticks([0,0.2,0.4,0.6,0.8,1])
  568. plt.savefig('Plots/SVG/violin_absolute.SVG', dpi = 300,bbox_inches='tight')
  569. plt.savefig('Plots/PNG/violin_absolute.PNG', dpi = 300,bbox_inches='tight')
  571. # %%
  572. export_plot_data('../../Source_Data.xlsx','Figure4f',CategoryGNG_popA=np.abs(S[0,popA]),CategoryGNG_popB=np.abs(S[0,popB]),
  573. ChoiceGNG_popA=np.abs(S[1,popA]),ChoiceGNG_popB=np.abs(S[1,popB]),
  574. CategoryLR_popA=np.abs(S[2,popA]),CategoryLR_popB=np.abs(S[2,popB]),
  575. ChoiceLR_popA=np.abs(S[3,popA]),ChoiceLR_popB=np.abs(S[3,popB]))
  577. # %%
  578. ## Absolute selectivity across task variables for different populations
  580. plt.figure(figsize=(12,8))
  582. plot_dfAB = plot_df[(plot_df['Regressor'] != 'Task')&(plot_df['Regressor'] != 'Lick LR')]
  584. p = sns.violinplot(data = plot_dfAB, x='Regressor',y='Abs selectivity',hue='Cluster',\
  585. hue_order=['A','B','else'],alpha=0.7,palette=['blue','red','green'],cut=0,common_norm=False,\
  586. inner_kws=dict(box_width=7, whis_width=2))
  589. plt.ylabel(r'|$\beta$|',fontsize=40)
  590. plt.xlabel('')
  591. plt.legend([],[], frameon=False)
  592. plt.axvline(x = 1.5, color = 'black', linewidth=1.6)
  593. plt.ylim(0,1)
  594. plt.yticks([0,0.2,0.4,0.6,0.8,1])
  596. plt.savefig('Plots/SVG/violin_absolute_all.SVG', dpi = 300,bbox_inches='tight')
  597. plt.savefig('Plots/PNG/violin_absolute_all.PNG', dpi = 300,bbox_inches='tight')
  599. # %%
  600. from scipy.stats import mannwhitneyu
  602. for reg_label in plot_dfAB['Regressor'].unique() :
  603. one_reg = plot_dfAB[(plot_dfAB['Regressor'] == reg_label)]
  604. print(reg)
  605. print(mannwhitneyu(one_reg[one_reg['Cluster'] == 'A']['Abs selectivity'],\
  606. one_reg[one_reg['Cluster'] == 'B']['Abs selectivity']))
  609. # %%
  610. plt.figure(figsize=(6,10))
  612. plot_dfall = plot_df
  614. task_only = True
  615. if task_only :
  616. plot_dfall = plot_dfall[plot_dfall['Regressor'] == 'Task']
  618. plt.figure(figsize=(2.5,8))
  620. p = sns.violinplot(data = plot_dfall, x='Regressor',y='Selectivity',hue='Cluster',hue_order=['A','B','else'],\
  621. alpha=0.7,palette=['blue','red','green'],inner_kws=dict(box_width=7, whis_width=2))
  623. plt.ylabel(r'$\beta_i$', fontsize=35)
  624. plt.xlabel('')
  625. plt.ylim((-1.2,1.2))
  626. plt.legend([],[], frameon=False)
  628. if task_only :
  630. plt.axhline(0,color='black',linestyle='--',linewidth=1.5)
  631. plt.ylabel(r'$\beta_{Task}$', fontsize=32)
  632. plt.xticks([])
  633. ax = plt.gca()
  635. ax.spines['bottom'].set_visible(False)
  636. ax.get_xaxis().set_visible(False)
  637. plt.savefig('Plots/SVG/violin_task.SVG', dpi = 300,bbox_inches='tight')
  638. plt.savefig('Plots/PNG/violin_task.PNG', dpi = 300,bbox_inches='tight')
  641. else :
  642. plt.axvline(x = 1.5, color = 'black', linewidth=1.5)
  643. plt.axvline(x = 4.5, color = 'black', linewidth=1.5)
  645. plt.savefig('Plots/SVG/violin.SVG', dpi = 300,bbox_inches='tight')
  646. plt.savefig('Plots/PNG/violin.PNG', dpi = 300,bbox_inches='tight')
  648. # %%
  649. print(mannwhitneyu(S[5,popA],S[5,popB],alternative = 'two-sided'))
  651. # %%
  652. export_plot_data('../../Source_Data.xlsx','Figure4b',popA = S[5,popA], popB = S[5,popB], pop0 = S[5,pop0])
  654. # %% [markdown]
  655. # ### Variance across time
  657. # %%
  658. fig, ax = plt.subplots(figsize=(7,7))
  659. M = np.stack(merged_data['Trajectory'])
  660. t_GNG = merged_data['Context'] == 0
  661. t_LR = merged_data['Context'] == 1
  663. var_GNG = np.std(M[t_GNG,:,:],axis=(0,1))
  664. var_LR = np.std(M[t_LR,:,:],axis=(0,1))
  666. nb_neurons = 200
  668. ax.scatter(var_GNG[pop0][:nb_neurons],var_LR[pop0][:nb_neurons],color='green',s=35, alpha = 0.65,\
  669. edgecolors='white',lw=1,antialiased=True,label='Population 0')
  670. ax.scatter(var_GNG[popA][:nb_neurons],var_LR[popA][:nb_neurons],color='blue',s=35, alpha = 0.65,\
  671. edgecolors='white',lw=1,antialiased=True,label='Population GNG')
  672. ax.scatter(var_GNG[popB][:nb_neurons],var_LR[popB][:nb_neurons],color='red',s=35, alpha = 0.65,\
  673. edgecolors='white',lw=1,antialiased=True,label='Population LR')
  675. confidence_ellipse(var_GNG[popA][:nb_neurons],var_LR[popA][:nb_neurons],ax, edgecolor='blue', n_std=1, alpha = 0.6,linewidth=3, zorder=3)
  676. confidence_ellipse(var_GNG[popB][:nb_neurons],var_LR[popB][:nb_neurons],ax, edgecolor='red', n_std=1, alpha = 0.6,linewidth=3, zorder=3)
  677. confidence_ellipse(var_GNG[pop0][:nb_neurons],var_LR[pop0][:nb_neurons],ax, edgecolor='green', n_std=1, alpha = 0.6,linewidth=3, zorder=3)
  678. plt.xlim((0,0.035))
  679. plt.ylim((0,0.035))
  680. plt.plot(np.linspace(0,0.035,100),np.linspace(0,0.035,100),alpha=0.5,color='black',linestyle='--')
  681. plt.xlabel('Time variance in Go/NoGo')
  682. plt.ylabel('Time variance in Left/Right')
  683. plt.xticks([0,0.035],[0,0.035])
  684. plt.yticks([0,0.035],[0,0.035])
  686. lgnd = plt.legend(frameon=False,fontsize=18)
  688. #plt.savefig('Plots/PNG/variance_time_gm.PNG', dpi = 300, bbox_inches='tight')
  689. #plt.savefig('Plots/SVG/variance_time_gm.SVG', dpi = 300, bbox_inches='tight')
  691. # %%
  692. plt.figure(figsize=(6,5))
  694. u = np.array([1/np.sqrt(2),-1/np.sqrt(2)])
  696. M_A = np.array([var_GNG[popA][:nb_neurons],var_LR[popA][:nb_neurons]]).T
  697. M_B = np.array([var_GNG[popB][:nb_neurons],var_LR[popB][:nb_neurons]]).T
  698. M_0 = np.array([var_GNG[pop0][:nb_neurons],var_LR[pop0][:nb_neurons]]).T
  700. ymin,ymax = -0.01, 0.01
  702. nb_bin=20
  704. hA = plt.hist(M_A@u,color='blue',density=True,range=(ymin,ymax), bins=nb_bin,alpha = 0.3,edgecolor = "black");
  705. hB = plt.hist(M_B@u,color='red',density=True, range=(ymin,ymax), bins=nb_bin,alpha = 0.3,edgecolor = "black");
  706. h0 = plt.hist(M_0@u,color='green',density=True,range=(ymin,ymax), bins=nb_bin,alpha = 0.3,edgecolor = "black");
  709. plt.clf()
  711. spl_A = splrep(hA[1], list(hA[0]) + [0], s=0.01, per=False)
  712. spl_B = splrep(hB[1], list(hB[0]) + [0], s=0.01, per=False)
  713. spl_0 = splrep(h0[1], list(h0[0]) + [0], s=0.01, per=False)
  714. x = np.linspace(ymin, ymax, 200)
  715. yA = splev(x, spl_A)
  716. yB = splev(x, spl_B)
  717. y0 = splev(x, spl_0)
  719. plt.fill_between(x,y1=0,y2=yA,color='blue',alpha=0.5)
  720. plt.fill_between(x,y1=0,y2=yB,color='red',alpha=0.5)
  721. plt.fill_between(x,y1=0,y2=y0,color='green',alpha=0.5)
  723. plt.ylim((0,300))
  725. plt.yticks([])
  727. plt.ylabel("Density")
  728. plt.xlim((ymin,ymax))
  729. plt.xticks([0],[''])
  731. #plt.savefig('Plots/PNG/variance_time_hist.PNG', dpi = 300, bbox_inches='tight')
  732. #plt.savefig('Plots/SVG/variance_time_hist.SVG', dpi = 300, bbox_inches='tight')
  734. # %%
  735. nb_bin=25
  736. bin_v = np.array([(k+0.5)*np.pi/(2*nb_bin) for k in range(0,nb_bin)])
  738. thetas = np.arctan(var_GNG/var_LR)
  740. hA = plt.hist(thetas[popA],color='blue',density=True,range=(0,np.pi/2), bins=nb_bin,alpha = 0.4);
  741. hB = plt.hist(thetas[popB],color='red',density=True, range=(0,np.pi/2), bins=nb_bin,alpha = 0.4);
  742. h0 = plt.hist(thetas[pop0],color='green',density=True,range=(0,np.pi/2), bins=nb_bin,alpha = 0.4);
  744. plt.clf()
  746. # %%
  747. from scipy.optimize import curve_fit
  749. plt.figure(figsize=(7,7))
  751. plt.hist(thetas[popA],color='red',density=True,range=(0,np.pi/2), bins=nb_bin, alpha = 0.3,edgecolor = "black")
  752. #plt.plot(bin_v,fit_y_popA,color='red',linewidth=3)
  754. plt.hist(thetas[popB],color='blue',density=True,range=(0,np.pi/2), bins=nb_bin, alpha = 0.3,edgecolor = "black")
  755. #plt.plot(bin_v,fit_y_popB,color='blue',linewidth=3)
  757. plt.hist(thetas[pop0],color='green',density=True,range=(0,np.pi/2), bins=nb_bin, alpha = 0.3,edgecolor = "black")
  758. #plt.plot(bin_v,fit_y_pop0,color='green',linewidth=3)
  760. plt.xlabel(r'$\alpha$',fontsize=30)
  761. plt.xticks([0,np.pi/4,np.pi/2],['0',r'$\frac{\pi}{4}$',r'$\frac{\pi}{2}$'])
  762. plt.ylabel('Density',fontsize=25)
  763. plt.xlim(0,np.pi/2)
  764. plt.yticks([0,1,2,3])
  766. """
  767. plt.savefig('Plots/PNG/variance_time_angle_gm.PNG', dpi = 300, bbox_inches='tight')
  768. plt.savefig('Plots/SVG/variance_time_angle_gm.SVG', dpi = 300, bbox_inches='tight')
  769. """
  771. # %%
  772. from scipy.stats import ttest_ind
  774. print(ttest_ind(thetas[popA],thetas[pop0]))
  775. print(ttest_ind(thetas[popB],thetas[pop0]))
  776. print(ttest_ind(thetas[popA],thetas[popB]))
  778. # %% [markdown]
  779. # ### Population decoding
  781. # %%
  782. M = np.stack(merged_data['Spike rate'])
  783. M = pca1.transform(M)
  784. M = pca1.inverse_transform(M)
  785. M = standardize(M)
  787. hit_GNG = merged_data['Label'] == 1
  788. cr_GNG = merged_data['Label'] == 3
  790. left_LR = merged_data['Label'] == 5
  791. right_LR = merged_data['Label'] == 8
  793. GNG_task = merged_data['Context'] == 0
  794. LR_task = 1-GNG_task
  796. M_GNG_hit = M[GNG_task&hit_GNG,:]
  797. M_GNG_cr = M[GNG_task&cr_GNG,:]
  799. M_LR_left = M[~GNG_task&left_LR,:]
  800. M_LR_right = M[~GNG_task&right_LR,:]
  802. X_GNG = np.concatenate((M_GNG_hit,M_GNG_cr),axis=0)
  803. X_LR = np.concatenate((M_LR_right,M_LR_left),axis=0)
  806. Y_GNG = np.array([1 for k in range(len(M_GNG_hit))] + [-1 for k in range(len(M_GNG_cr))])
  807. Y_LR = np.array([1 for k in range(len(M_LR_right))] + [-1 for k in range(len(M_LR_left))])
  809. # %%
  810. from tqdm.auto import tqdm
  811. from sklearn.svm import SVC
  814. nb_decoding = 100
  815. nb_neuron = 100
  816. def pop_decoding(nb_neuron = nb_neuron, nb_decoding = nb_decoding, return_inter_task=False) :
  817. """Method performing population decoding with Support Vector Classifier trained on nb_neuron different cells .
  818. We train and test models on every combinations of epochs"""
  819. scores_A_GNG = []
  820. scores_A_LR = []
  822. scores_B_GNG = []
  823. scores_B_LR = []
  825. scores_0_GNG = []
  826. scores_0_LR = []
  828. scores_A_GNG_to_LR = []
  829. scores_A_LR_to_GNG = []
  831. scores_B_GNG_to_LR = []
  832. scores_B_LR_to_GNG = []
  834. scores_0_GNG_to_LR = []
  835. scores_0_LR_to_GNG = []
  837. masks = []
  838. for n in tqdm(range(nb_decoding)) :
  839. svc_GNG_popA = SVC(kernel='linear',gamma='auto')
  840. svc_GNG_popB = SVC(kernel='linear',gamma='auto')
  841. svc_GNG_pop0 = SVC(kernel='linear',gamma='auto')
  843. svc_LR_popA = SVC(kernel='linear',gamma='auto')
  844. svc_LR_popB = SVC(kernel='linear',gamma='auto')
  845. svc_LR_pop0 = SVC(kernel='linear',gamma='auto')
  847. mask_A = np.random.choice(range(np.sum(popA)),nb_neuron,replace = False) # We do not replace to ensure independent decoding scores
  848. mask_B = np.random.choice(range(np.sum(popB)),nb_neuron,replace = False)
  849. mask_0 = np.random.choice(range(np.sum(pop0)),nb_neuron,replace = False)
  851. split_GNG = np.random.choice(range(len(X_GNG)),len(X_GNG)//2,replace=False)
  852. split_LR = np.random.choice(range(len(X_LR)),len(X_LR)//2,replace=False)
  854. X_GNG_train = X_GNG[split_GNG,:]
  855. Y_GNG_train = Y_GNG[split_GNG]
  857. X_GNG_test = np.delete(X_GNG,split_GNG,axis=0)
  858. Y_GNG_test = np.delete(Y_GNG,split_GNG,axis=0)
  860. X_LR_train = X_LR[split_LR,:]
  861. Y_LR_train = Y_LR[split_LR]
  863. X_LR_test = np.delete(X_LR,split_LR,axis=0)
  864. Y_LR_test = np.delete(Y_LR,split_LR,axis=0)
  867. svc_GNG_popA.fit(X_GNG_train[:,popA][:,mask_A],Y_GNG_train)
  868. svc_GNG_popB.fit(X_GNG_train[:,popB][:,mask_B],Y_GNG_train)
  869. svc_GNG_pop0.fit(X_GNG_train[:,pop0][:,mask_0],Y_GNG_train)
  871. svc_LR_popA.fit(X_LR_train[:,popA][:,mask_A],Y_LR_train)
  872. svc_LR_popB.fit(X_LR_train[:,popB][:,mask_B],Y_LR_train)
  873. svc_LR_pop0.fit(X_LR_train[:,pop0][:,mask_0],Y_LR_train)
  875. scores_A_GNG.append(svc_GNG_popA.score(X_GNG_test[:,popA][:,mask_A],Y_GNG_test))
  876. scores_B_GNG.append(svc_GNG_popB.score(X_GNG_test[:,popB][:,mask_B],Y_GNG_test))
  877. scores_0_GNG.append(svc_GNG_pop0.score(X_GNG_test[:,pop0][:,mask_0],Y_GNG_test))
  879. scores_A_LR.append(svc_LR_popA.score(X_LR_test[:,popA][:,mask_A],Y_LR_test))
  880. scores_B_LR.append(svc_LR_popB.score(X_LR_test[:,popB][:,mask_B],Y_LR_test))
  881. scores_0_LR.append(svc_LR_pop0.score(X_LR_test[:,pop0][:,mask_0],Y_LR_test))
  883. masks.append([mask_A, mask_B, mask_0])
  885. if return_inter_task :
  887. scores_A_GNG_to_LR.append(svc_GNG_popA.score(X_LR_test[:,popA][:,mask_A],Y_LR_test))
  888. scores_B_GNG_to_LR.append(svc_GNG_popB.score(X_LR_test[:,popB][:,mask_B],Y_LR_test))
  889. scores_0_GNG_to_LR.append(svc_GNG_pop0.score(X_LR_test[:,pop0][:,mask_0],Y_LR_test))
  891. scores_A_LR_to_GNG.append(svc_LR_popA.score(X_GNG_test[:,popA][:,mask_A],Y_GNG_test))
  892. scores_B_LR_to_GNG.append(svc_LR_popB.score(X_GNG_test[:,popB][:,mask_B],Y_GNG_test))
  893. scores_0_LR_to_GNG.append(svc_LR_pop0.score(X_GNG_test[:,pop0][:,mask_0],Y_GNG_test))
  895. scores_A_GNG = np.array(scores_A_GNG)
  896. scores_A_LR = np.array(scores_A_LR)
  898. scores_B_GNG = np.array(scores_B_GNG)
  899. scores_B_LR = np.array(scores_B_LR)
  901. scores_0_GNG = np.array(scores_0_GNG)
  902. scores_0_LR = np.array(scores_0_LR)
  904. scores_A_GNG_to_LR = np.array(scores_A_GNG_to_LR)
  905. scores_A_LR_to_GNG = np.array(scores_A_LR_to_GNG)
  907. scores_B_GNG_to_LR = np.array(scores_B_GNG_to_LR)
  908. scores_B_LR_to_GNG = np.array(scores_B_LR_to_GNG)
  910. scores_0_GNG_to_LR = np.array(scores_0_GNG_to_LR)
  911. scores_0_LR_to_GNG = np.array(scores_0_LR_to_GNG)
  913. masks = np.array(masks)
  915. if return_inter_task :
  916. return [scores_A_GNG, scores_A_LR, scores_A_GNG_to_LR, scores_A_LR_to_GNG],\
  917. [scores_B_GNG, scores_B_LR, scores_B_GNG_to_LR, scores_B_LR_to_GNG],\
  918. [scores_0_GNG, scores_0_LR, scores_0_GNG_to_LR, scores_0_LR_to_GNG], masks
  920. else :
  921. return [scores_A_GNG, scores_A_LR], [scores_B_GNG, scores_B_LR], [scores_0_GNG, scores_0_LR]
  923. # %%
  924. scores_A, scores_B, scores_0, _ = pop_decoding(2,100,return_inter_task=True)
  925. scores_A_many, scores_B_many, scores_0_many, _ = pop_decoding(20,return_inter_task=True)
  927. # %% [markdown]
  928. # ### Go/NoGo - Go/NoGo
  930. # %%
  931. from scipy.interpolate import splev, splrep
  933. fig,axs = plt.subplots(1,2,figsize=(6,8),gridspec_kw={'width_ratios': [3, 1]})
  935. #fig.subplots_adjust(hspace=0)
  937. df_dict = {'Population\nGNG' : scores_A[0], 'Population\nLR' : scores_B[0]}
  938. test = pd.DataFrame(df_dict)
  939. sns.stripplot(test,alpha=0.7,palette=['blue','red'],jitter=0.2,size=12,ax=axs[0])
  940. axs[0].errorbar([0,1],[np.mean(scores_A[0]),np.mean(scores_B[0])],yerr=[np.std(scores_A[0])/np.sqrt(nb_decoding),np.std(scores_B[0])/np.sqrt(nb_decoding)],\
  941. color='black',marker='_',markersize=85,elinewidth=3.5,capsize=7,zorder=10,linestyle='',markeredgewidth=3.5)
  943. ymin,ymax = 0.3,1
  945. axs[0].set_xlim((-0.5,1.5))
  946. axs[0].set_ylim((ymin,ymax+0.05))
  947. axs[0].set_yticks([0.5,0.75,1],['50%','75%','100%'])
  948. axs[0].set_ylabel('Population decoding accuracy')
  949. axs[0].fill_between([-0.5,1.5],ymin,0.5, color='black',alpha=0.05)
  950. axs[0].axhline(0.5,color="black",linestyle='--',alpha=0.3,linewidth=1.5)
  955. nb_bins=25
  956. histA = axs[1].hist(scores_A[0],range=(ymin,ymax),bins=nb_bins,orientation='horizontal',color='blue',alpha=0.5,density=True)
  957. histB = axs[1].hist(scores_B[0],range=(ymin,ymax),bins=nb_bins,orientation='horizontal',color='red',alpha=0.5,density=True)
  959. axs[1].clear()
  961. axs[1].set_yticks([])
  962. axs[1].set_xticks([])
  963. axs[1].spines['bottom'].set_visible(False)
  964. axs[1].set_ylim((ymin,ymax+0.02))
  965. axs[1].set_xlim((0,70))
  966. # Interpolation
  967. spl_A = splrep(histA[1], list(histA[0]) + [0], s=0.01, per=False)
  968. spl_B = splrep(histB[1], list(histB[0]) + [0], s=0.01, per=False)
  969. x = np.linspace(ymin, ymax, 200)
  970. yA = splev(x, spl_A)
  971. yB = splev(x, spl_B)
  973. axs[1].fill_betweenx(x,x1=0,x2=yA,color='blue',alpha=0.5)
  974. axs[1].fill_betweenx(x,x1=0,x2=yB,color='red',alpha=0.5)
  977. #plt.savefig('Plots/SVG/decoding_2neuron_GNG.SVG', dpi = 300,bbox_inches='tight')
  978. #plt.savefig('Plots/PNG/decoding_2neuron_GNG.PNG', dpi = 300,bbox_inches='tight')
  980. print(mannwhitneyu(scores_A[0],scores_B[0]))
  982. # %%
  983. from scipy.interpolate import splev, splrep
  985. fig,ax = plt.subplots(1,1,figsize=(6.2,8))
  987. #fig.subplots_adjust(hspace=0)
  989. df_dict = {'Population\nGNG' : scores_A[0], 'Population\nLR' : scores_B[0], 'Population\n0' : scores_0[0]}
  990. test = pd.DataFrame(df_dict)
  991. sns.stripplot(test,alpha=0.7,palette=['blue','red','green'],jitter=0.2,size=12,ax=ax)
  992. ax.errorbar([0,1,2],[np.mean(scores_A[0]),np.mean(scores_B[0]),np.mean(scores_0[0])],yerr=[np.std(scores_A[0])/np.sqrt(nb_decoding),np.std(scores_B[0])/np.sqrt(nb_decoding),np.std(scores_0[0])/np.sqrt(nb_decoding)],\
  993. color='black',marker='_',markersize=85,elinewidth=3.5,capsize=7,zorder=10,linestyle='',markeredgewidth=3.5)
  995. ymin,ymax = 0.3,1
  997. ax.set_xlim((-0.5,2.5))
  998. ax.set_ylim((ymin,ymax+0.05))
  999. ax.set_yticks([0.5,0.75,1],['50%','75%','100%'])
  1000. ax.set_ylabel('Population decoding accuracy')
  1001. ax.fill_between([-0.5,2.5],ymin,0.5, color='black',alpha=0.05)
  1002. ax.axhline(0.5,color="black",linestyle='--',alpha=0.3,linewidth=1.5)
  1005. plt.savefig('Plots/SVG/decoding_2neuron_GNG.SVG', dpi = 300,bbox_inches='tight')
  1006. plt.savefig('Plots/PNG/decoding_2neuron_GNG.PNG', dpi = 300,bbox_inches='tight')
  1008. print(mannwhitneyu(scores_A[0],scores_B[0]))
  1009. print(mannwhitneyu(scores_A[0],scores_0[0]))
  1010. print(mannwhitneyu(scores_B[0],scores_0[0]))
  1012. # %%
  1013. export_plot_data('../../Source_Data.xlsx','Figure4g',scores_A = scores_A[0], scores_B = scores_B[0], scores_0 = scores_0[0])
  1015. # %% [markdown]
  1016. # ### LR - LR
  1018. # %%
  1019. from scipy.interpolate import splev, splrep
  1021. fig,axs = plt.subplots(1,2,figsize=(6,8),gridspec_kw={'width_ratios': [3, 1]})
  1022. #fig.subplots_adjust(hspace=0)
  1024. df_dict = {'Population\nGNG' : scores_A[1], 'Population\nLR' : scores_B[1]}
  1025. test = pd.DataFrame(df_dict)
  1026. sns.stripplot(test,alpha=0.7,palette=['blue','red'],jitter=0.2,size=12,ax=axs[0])
  1027. axs[0].errorbar([0,1],[np.mean(scores_A[1]),np.mean(scores_B[1])],yerr=[np.std(scores_A[1])/np.sqrt(nb_decoding),np.std(scores_B[1])/np.sqrt(nb_decoding)],\
  1028. color='black',marker='_',markersize=85,elinewidth=3.5,capsize=7,zorder=10,linestyle='',markeredgewidth=3.5)
  1030. ymin,ymax = 0.3,1
  1032. axs[0].set_xlim((-0.5,1.5))
  1033. axs[0].set_ylim((ymin,ymax+0.05))
  1034. axs[0].set_yticks([0.5,0.75,1],['50%','75%','100%'])
  1035. axs[0].set_ylabel('Population decoding accuracy')
  1036. axs[0].fill_between([-0.5,1.5],ymin,0.5, color='black',alpha=0.05)
  1037. axs[0].axhline(0.5,color="black",linestyle='--',alpha=0.3,linewidth=1.5)
  1040. nb_bins=25
  1041. histA = axs[1].hist(scores_A[1],range=(ymin,ymax),bins=nb_bins,orientation='horizontal',color='blue',alpha=0.5,density=True)
  1042. histB = axs[1].hist(scores_B[1],range=(ymin,ymax),bins=nb_bins,orientation='horizontal',color='red',alpha=0.5,density=True)
  1044. axs[1].clear()
  1047. axs[1].set_yticks([])
  1048. axs[1].set_xticks([])
  1049. axs[1].spines['bottom'].set_visible(False)
  1050. axs[1].set_ylim((ymin,ymax+0.02))
  1051. axs[1].set_xlim((0,30))
  1052. # Interpolation
  1053. spl_A = splrep(histA[1], list(histA[0]) + [0], s=0.01, per=False)
  1054. spl_B = splrep(histB[1], list(histB[0]) + [0], s=0.01, per=False)
  1055. x = np.linspace(ymin, ymax, 200)
  1056. yA = splev(x, spl_A)
  1057. yB = splev(x, spl_B)
  1059. axs[1].fill_betweenx(x,x1=0,x2=yA,color='blue',alpha=0.5)
  1060. axs[1].fill_betweenx(x,x1=0,x2=yB,color='red',alpha=0.5)
  1063. #plt.savefig('Plots/SVG/decoding_2neuron_LR.SVG', dpi = 300,bbox_inches='tight')
  1064. #plt.savefig('Plots/PNG/decoding_2neuron_LR.PNG', dpi = 300,bbox_inches='tight')
  1066. print(mannwhitneyu(scores_A[1],scores_B[1]))
  1068. # %%
  1069. from scipy.interpolate import splev, splrep
  1071. fig,ax = plt.subplots(1,1,figsize=(6.2,8))
  1073. df_dict = {'Population\nGNG' : scores_A[1], 'Population\nLR' : scores_B[1], 'Population\n0' : scores_0[1]}
  1074. test = pd.DataFrame(df_dict)
  1075. sns.stripplot(test,alpha=0.7,palette=['blue','red','green'],jitter=0.2,size=12,ax=ax)
  1076. ax.errorbar([0,1,2],[np.mean(scores_A[1]),np.mean(scores_B[1]),np.mean(scores_0[1])],yerr=[np.std(scores_A[1])/np.sqrt(nb_decoding),np.std(scores_B[1])/np.sqrt(nb_decoding),np.std(scores_0[1])/np.sqrt(nb_decoding)],\
  1077. color='black',marker='_',markersize=85,elinewidth=3.5,capsize=7,zorder=10,linestyle='',markeredgewidth=3.5)
  1079. ymin,ymax = 0.3,1
  1081. ax.set_xlim((-0.5,2.5))
  1082. ax.set_ylim((ymin,ymax+0.05))
  1083. ax.set_yticks([0.5,0.75,1],['50%','75%','100%'])
  1084. ax.set_ylabel('Population decoding accuracy')
  1085. ax.fill_between([-0.5,2.5],ymin,0.5, color='black',alpha=0.05)
  1086. ax.axhline(0.5,color="black",linestyle='--',alpha=0.3,linewidth=1.5)
  1088. plt.savefig('Plots/SVG/decoding_2neuron_LR.SVG', dpi = 300,bbox_inches='tight')
  1089. plt.savefig('Plots/PNG/decoding_2neuron_LR.PNG', dpi = 300,bbox_inches='tight')
  1091. print(mannwhitneyu(scores_A[1],scores_B[1]))
  1092. print(mannwhitneyu(scores_A[1],scores_0[1]))
  1093. print(mannwhitneyu(scores_B[1],scores_0[1]))
  1095. # %%
  1096. export_plot_data('../../Source_Data.xlsx','Figure4h',scores_A = scores_A[1], scores_B = scores_B[1], scores_0 = scores_0[1])
  1098. # %% [markdown]
  1099. # ### Go/NoGo - LR
  1101. # %%
  1102. from scipy.interpolate import splev, splrep
  1104. fig,axs = plt.subplots(1,2,figsize=(10,8),gridspec_kw={'width_ratios': [3, 1]})
  1105. #fig.subplots_adjust(hspace=0)
  1107. df_dict = {'Population\nGNG' : scores_A_many[2], 'Population\nLR' : scores_B_many[2],'Population\n0' : scores_0_many[2]}
  1108. df_plot = pd.DataFrame(df_dict)
  1109. sns.stripplot(df_plot,alpha=0.7,palette=['blue','red','green'],jitter=0.2,size=12,ax=axs[0])
  1110. axs[0].errorbar([0,1,2],[np.mean(scores_A_many[2]),np.mean(scores_B_many[2]),np.mean(scores_0_many[2])],\
  1111. yerr=[np.std(scores_A_many[2])/np.sqrt(nb_decoding),np.std(scores_B_many[2])/np.sqrt(nb_decoding),np.std(scores_0_many[2])/np.sqrt(nb_decoding)],\
  1112. color='black',marker='_',markersize=80,elinewidth=3.5,capsize=7,zorder=10,linestyle='',markeredgewidth=3.5)
  1114. ymin,ymax = 0.3,1
  1116. axs[0].set_xlim((-0.5,2.5))
  1117. axs[0].set_ylim((ymin,ymax+0.05))
  1118. axs[0].set_yticks([0.5,0.75,1],['50%','75%','100%'],fontsize=25)
  1119. axs[0].set_ylabel('Population decoding accuracy')
  1120. axs[0].fill_between([-0.5,2.5],ymin,0.5, color='black',alpha=0.05)
  1121. axs[0].axhline(0.5,color="black",linestyle='--',alpha=0.3,linewidth=1.5)
  1127. nb_bins=25
  1128. histA = axs[1].hist(scores_A_many[2],range=(ymin,ymax),bins=nb_bins,orientation='horizontal',color='blue',alpha=0.5,density=True)
  1129. histB = axs[1].hist(scores_B_many[2],range=(ymin,ymax),bins=nb_bins,orientation='horizontal',color='red',alpha=0.5,density=True)
  1130. hist0 = axs[1].hist(scores_0_many[2],range=(ymin,ymax),bins=nb_bins,orientation='horizontal',color='green',alpha=0.5,density=True)
  1132. axs[1].clear()
  1135. axs[1].set_yticks([])
  1136. axs[1].set_xticks([])
  1137. axs[1].spines['bottom'].set_visible(False)
  1138. axs[1].set_ylim((ymin,ymax+0.02))
  1139. axs[1].set_xlim((0,30))
  1140. # Interpolation
  1141. spl_A = splrep(histA[1], list(histA[0]) + [0], s=0.01, per=False)
  1142. spl_B = splrep(histB[1], list(histB[0]) + [0], s=0.01, per=False)
  1143. spl_0 = splrep(hist0[1], list(hist0[0]) + [0], s=0.01, per=False)
  1144. x = np.linspace(ymin, ymax, 200)
  1145. yA = splev(x, spl_A)
  1146. yB = splev(x, spl_B)
  1147. y0 = splev(x, spl_0)
  1148. axs[1].fill_betweenx(x,x1=0,x2=yA,color='blue',alpha=0.5)
  1149. axs[1].fill_betweenx(x,x1=0,x2=yB,color='red',alpha=0.5)
  1150. axs[1].fill_betweenx(x,x1=0,x2=y0,color='green',alpha=0.5)
  1152. #plt.savefig('Plots/SVG/decoding_20neuron_GNG_to_LR.SVG', dpi = 300,bbox_inches='tight')
  1153. #plt.savefig('Plots/PNG/decoding_20neuron_GNG_to_LR.PNG', dpi = 300,bbox_inches='tight')
  1155. print(mannwhitneyu(scores_0_many[2],scores_A_many[2]))
  1156. print(mannwhitneyu(scores_0_many[2],scores_B_many[2]))
  1158. print(wilcoxon(scores_0_many[2]-0.5,alternative='greater'))
  1160. # %%
  1161. export_plot_data('../../Source_Data.xlsx','Figure4i',scores_A = scores_A_many[2], scores_B = scores_B_many[2], scores_0 = scores_0_many[2])
  1163. # %% [markdown]
  1164. # ### Information
  1166. # %%
  1167. ## Exploration on the nature of information each population provides
  1169. nb_bins = 20
  1170. def compute_prob(data) :
  1171. hist = np.histogram(data,nb_bins,range=(-2,2))[0]
  1172. return hist/np.sum(hist)
  1174. eps = 0
  1175. def MI(XY_mat,pX,pY) :
  1176. # XY matrix must be the matrix of X conditionned on Y
  1177. mi = 0
  1178. for i in range(np.shape(XY_mat)[0]) :
  1179. for j in range(np.shape(XY_mat)[1]) :
  1180. if (pY[i] <= eps) or (pX[j] <= eps) or (XY_mat[i,j] <= eps) :
  1181. mi+=0
  1182. else :
  1183. #print(pY[i])
  1184. #print(pX[j])
  1185. mi+=pY[i]*XY_mat[i,j]*np.log2(pY[i]*XY_mat[i,j]/(pY[i]*pX[j]))
  1186. return mi
  1188. def cond_MI(XYZ_mat,XZ_mat,pX,pY,pZ) :
  1189. # We compute the MI between X and Y conditionned on Z
  1190. # Here we suppose the independence between Y and Z. XYZ matrix is the matrix of X conditionned on Y and Z.
  1191. mi = 0
  1192. for i in range(np.shape(XYZ_mat)[0]) :
  1193. for j in range(np.shape(XYZ_mat)[1]) :
  1194. for k in range(np.shape(XYZ_mat)[2]) :
  1195. if (XZ_mat[j,k] <= eps) or (XYZ_mat[i,j,k] <= eps):
  1196. mi += 0
  1197. else :
  1198. mi += pY[i]*pZ[j]*XYZ_mat[i,j,k]*np.log2(XYZ_mat[i,j,k]/XZ_mat[j,k])
  1200. return mi
  1202. # %%
  1203. from sklearn.mixture import GaussianMixture, BayesianGaussianMixture
  1204. from tqdm.auto import tqdm
  1206. def get_pop(prob,S) :
  1207. pop1 = prob[:,0] >= 0.9
  1208. pop2 = prob[:,1] >= 0.9
  1209. pop3 = prob[:,2] >= 0.9
  1211. pops = [pop1,pop2,pop3]
  1213. S_pop = [S[:,pop1],S[:,pop2],S[:,pop3]]
  1215. mean_ctx_s = [np.mean(S_pop[0][0,:]),np.mean(S_pop[1][0,:]),np.mean(S_pop[2][0,:])]
  1217. popA = pops[np.argsort(mean_ctx_s)[0]]
  1218. popB = pops[np.argsort(mean_ctx_s)[2]]
  1219. pop0 = pops[np.argsort(mean_ctx_s)[1]]
  1221. return popA, popB, pop0
  1223. def compute_UI_SYN_pop(S_grid,A_data,s_data,c_data) :
  1224. pC = np.array([0.5,0.5])
  1225. pS = np.array([0.5,0.5])
  1227. MI_AC_popA = []
  1228. MI_AS_popA = []
  1229. MI_AC_S_popA = []
  1231. MI_AC_popB = []
  1232. MI_AS_popB = []
  1233. MI_AC_S_popB = []
  1235. MI_AC_pop0 = []
  1236. MI_AS_pop0 = []
  1237. MI_AC_S_pop0 = []
  1239. A_matrices = np.array([compute_prob(A_data[:,n]) for n in range(np.shape(A_data)[1])])
  1240. AS_matrices = np.array([[compute_prob(A_data[s_data==-1,n]),compute_prob(A_data[s_data==1,n])] for n in range(np.shape(A_data)[1])])
  1241. AC_matrices = np.array([[compute_prob(A_data[c_data==-1,n]),compute_prob(A_data[c_data==1,n])] for n in range(np.shape(A_data)[1])])
  1244. ACS_matrices = np.array([[[compute_prob(A_data[(s_data==-1)&(c_data==-1),n]),compute_prob(A_data[(s_data==1)&(c_data==-1),n])],\
  1245. [compute_prob(A_data[(s_data==-1)&(c_data==1),n]),compute_prob(A_data[(s_data==1)&(c_data==1),n])]] for n in range(np.shape(A_data)[1])])
  1247. MI_AS_s = []
  1248. MI_AC_s = []
  1249. MI_AC_S_s = []
  1250. for n in range(np.shape(A_matrices)[0]) :
  1251. pA = A_matrices[n]
  1252. AS_mat = AS_matrices[n]
  1253. AC_mat = AC_matrices[n]
  1254. ACS_mat = ACS_matrices[n]
  1256. MI_AS_s.append(MI(AS_mat,pA,pS))
  1257. MI_AC_s.append(MI(AC_mat,pA,pC))
  1258. MI_AC_S_s.append(cond_MI(ACS_mat,AS_mat,pA,pC,pS))
  1261. MI_AS_s = np.array(MI_AS_s)
  1262. MI_AC_s = np.array(MI_AC_s)
  1263. MI_AC_S_s = np.array(MI_AC_S_s)
  1265. MI_AC_popA.append(MI_AC_s[popA])
  1266. MI_AC_popB.append(MI_AC_s[popB])
  1267. MI_AC_pop0.append(MI_AC_s[pop0])
  1269. MI_AS_popA.append(MI_AS_s[popA])
  1270. MI_AS_popB.append(MI_AS_s[popB])
  1271. MI_AS_pop0.append(MI_AS_s[pop0])
  1273. MI_AC_S_popA.append(MI_AC_S_s[popA])
  1274. MI_AC_S_popB.append(MI_AC_S_s[popB])
  1275. MI_AC_S_pop0.append(MI_AC_S_s[pop0])
  1277. MI_AC_popA = np.array(MI_AC_popA)
  1278. MI_AC_popB = np.array(MI_AC_popB)
  1279. MI_AC_pop0 = np.array(MI_AC_pop0)
  1281. MI_AS_popA = np.array(MI_AS_popA)
  1282. MI_AS_popB = np.array(MI_AS_popB)
  1283. MI_AS_pop0 = np.array(MI_AS_pop0)
  1285. MI_AC_S_popA = np.array(MI_AC_S_popA)
  1286. MI_AC_S_popB = np.array(MI_AC_S_popB)
  1287. MI_AC_S_pop0 = np.array(MI_AC_S_pop0)
  1289. return np.squeeze(np.array([MI_AC_popA, MI_AS_popA, MI_AC_S_popA])), \
  1290. np.squeeze(np.array([MI_AC_popB, MI_AS_popB, MI_AC_S_popB])), \
  1291. np.squeeze(np.array([MI_AC_pop0, MI_AS_pop0, MI_AC_S_pop0]))
  1294. # %%
  1295. M = np.stack(merged_data['Spike rate'])
  1296. M = standardize(M)
  1297. M = pca1.inverse_transform(pca1.transform(M))
  1300. c_d = 2*merged_data['Context']-1
  1301. s_d = 2*merged_data['Category ID']-3
  1303. MI_popA, MI_popB, MI_pop0 = compute_UI_SYN_pop(S, M, s_d, c_d)
  1305. # %%
  1306. plt.figure(figsize=(10,7))
  1309. colors=['blue','red','green']
  1310. def set_color(bplot) :
  1311. for patch, color in zip(bplot['boxes'], colors):
  1312. patch.set_facecolor(color)
  1313. patch.set_alpha(0.3)
  1315. for median, color in zip(bplot['medians'], colors):
  1316. median.set_color(color)
  1317. median.set_linewidth(2)
  1319. for mean, color in zip(bplot['means'], colors):
  1320. mean.set_color(color)
  1321. mean.set_alpha(0.5)
  1322. mean.set_linewidth(1.5)
  1325. plt.ylim(-0.05,1.05)
  1329. bplot = plt.boxplot([MI_popA[0,:],MI_popB[0,:],MI_pop0[0,:]], positions=[0,0.3,0.6],\
  1330. patch_artist=True,notch=True,bootstrap=1000,sym='',showmeans=True,meanline=True)
  1331. set_color(bplot)
  1334. bplot = plt.boxplot([MI_popA[1,:],MI_popB[1,:],MI_pop0[1,:]], positions=[1.5,1.8,2.1],\
  1335. patch_artist=True,notch=True,bootstrap=1000,sym='',showmeans=True,meanline=True)
  1336. set_color(bplot)
  1339. bplot = plt.boxplot([MI_popA[2,:]-MI_popA[0,:],MI_popB[2,:]-MI_popB[0,:],MI_pop0[2,:]-MI_pop0[0,:]], positions=[3,3.3,3.6],\
  1340. patch_artist=True,notch=True,bootstrap=1000,sym='',showmeans=True,meanline=True)
  1341. set_color(bplot)
  1343. plt.xticks([0.3,1.8,3.3],[r'$UI_{Category}$',r'$UI_{Task}$',r'$Syn_{\{Category,Task\}}$'])
  1344. plt.ylabel('Information (bits)')
  1345. plt.yticks([0,0.2,0.4,0.6,0.8,1])
  1347. plt.savefig('Plots/SVG/MI_pop.SVG', dpi = 300,bbox_inches='tight')
  1348. plt.savefig('Plots/PNG/MI_pop.PNG', dpi = 300,bbox_inches='tight')
  1350. # %%
  1351. from scipy.stats import mannwhitneyu, permutation_test
  1353. print('------- UI category -------')
  1354. print('A-B : ' + str(mannwhitneyu(MI_popA[0,:],MI_popB[0,:])))
  1355. print('A-0 : ' + str(mannwhitneyu(MI_popA[0,:],MI_pop0[0,:])))
  1356. print('B-0 : ' + str(mannwhitneyu(MI_popB[0,:],MI_pop0[0,:])))
  1358. print('------- UI task -------')
  1359. print('A-B : ' + str(mannwhitneyu(MI_popA[1,:],MI_popB[1,:])))
  1360. print('A-0 : ' + str(mannwhitneyu(MI_popA[1,:],MI_pop0[1,:])))
  1361. print('B-0 : ' + str(mannwhitneyu(MI_popB[1,:],MI_pop0[1,:])))
  1363. print('------- Synergy -------')
  1364. print('A-B : ' + str(mannwhitneyu(MI_popA[2,:]-MI_popA[0,:],MI_popB[2,:]-MI_popB[0,:])))
  1365. print('A-0 : ' + str(mannwhitneyu(MI_popA[2,:]-MI_popA[0,:],MI_pop0[2,:]-MI_pop0[0,:])))
  1366. print('B-0 : ' + str(mannwhitneyu(MI_popB[2,:]-MI_popB[0,:],MI_pop0[2,:]-MI_pop0[0,:])))

gaussian_mixture.ipynb at commit 852d401, no license · at the source

Overview

Authors: Hugo Tissot1, Jeff Boucher1,2, Sandra Reinert3,4, Pieter M Goltstein3, Yves Boubenec1
ORCID iDs: Hugo Tissot, Jeff Boucher, Sandra Reinert, Pieter M Goltstein, Yves Boubenec
  1. Laboratoire des Systemes Perceptifs, Ecole Normale Superieure, PSL University, CNRS, Paris, France
  2. Ear Institute, University College London, London, UK
  3. Max Planck Institute for Biological Intelligence, Martinsried, Germany
  4. Present Address: Sainsbury Wellcome Centre, University College London, London, UK
Institutions: Centre National de la Recherche Scientifique (France); Université Paris Sciences et Lettres (France); École Normale Supérieure - PSL (France); Laboratoire des Systèmes Perceptifs (France); University College London (United Kingdom); Max Planck Institute for Biological Intelligence; Sainsbury Wellcome Centre (United Kingdom)
Journal: Nature communications, volume 17, issue 1, article 9904
Dates: received 21 August 2025; accepted 20 July 2026; published online 19 August 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1038/s41467-026-76104-3 · PMID 42749711 · PMCID PMC13582938 · OpenAlex W4414383722
Open access: gold, a free copy (OpenAlex)
Preprint: doi:10.1101/2025.09.22.677699
Status: code verified
Categories: computational modeling (no new data) (modality), non-human primate (organism), systems (subfield)
Methods: Statistics, Smoothing, state filtering, decompositions, Machine learning, Preprocessing, Single-unit activity, calcium imaging, Connectivity
Keywords: Neural encoding, Neural circuits, Network models, Decision
MeSH: Decision Making*, Prefrontal Cortex*, Animals, Macaca mulatta, Male, Models, Neurological, Neurons (* major topic)
Topic: Neural dynamics and brain function (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: Agence Nationale de la Recherche (French National Research Agency) (ANR-10-IDEX- 0001-02, ANR-17-EURE-0017, ANR-24-CE45-7453); Université de Recherche Paris Sciences et Lettres (PSL Research University) (PSL-NEURO, qLife)
Citations: not cited yet (Europe PMC); 51 references in the paper

Abstract

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

Repository

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

HTissot42/Sensorimotor-remapping-drives-task-specialization

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 852d401f2c575a49377ae0e7b91e2b309db210a5, 23 June 2026
Languages: Jupyter (66), Python (2)
Size: 345 files, 68 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, 34 notebooks
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (58 files), pandas (58 files), Matplotlib (45 files), SciPy (38 files), scikit-learn (37 files), seaborn (24 files), statsmodels (4 files), NetworkX (2 files)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
59 files
At the source: github.com/HTissot42/Sensorimotor-remapping-drives-task-specialization

Code availability statement

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

Read it in the paper: doi.org/10.1038/s41467-026-76104-3.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 58 scripts, each with its path and the digest of its content;
  • 21 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Data availability statement

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

Read it in the paper: doi.org/10.1038/s41467-026-76104-3.

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

Cite

This paper

Tissot, H., Boucher, J., Reinert, S., Goltstein, P. M., & Boubenec, Y. (2026). Sensorimotor remapping drives task specialization in prefrontal cortex. Nature communications, 17(1), 9904. https://doi.org/10.1038/s41467-026-76104-3

BibTeX

@article{tissot2026sensorimotor,
author = {Tissot, Hugo and Boucher, Jeff and Reinert, Sandra and Goltstein, Pieter M and Boubenec, Yves},
title = {{Sensorimotor remapping drives task specialization in prefrontal cortex}},
journal = {Nature communications},
year = {2026},
month = aug,
volume = {17},
number = {1},
pages = {9904},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-76104-3},
url = {https://doi.org/10.1038/s41467-026-76104-3},
pmid = {42749711},
pmcid = {PMC13582938}
}

RIS

TY - JOUR
AU - Tissot, Hugo
AU - Boucher, Jeff
AU - Reinert, Sandra
AU - Goltstein, Pieter M
AU - Boubenec, Yves
TI - Sensorimotor remapping drives task specialization in prefrontal cortex
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/08/19
VL - 17
IS - 1
SP - 9904
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-76104-3
UR - https://doi.org/10.1038/s41467-026-76104-3
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-76104-3",
"type": "article-journal",
"title": "Sensorimotor remapping drives task specialization in prefrontal cortex",
"container-title": "Nature communications",
"author": [
{
"family": "Tissot",
"given": "Hugo"
},
{
"family": "Boucher",
"given": "Jeff"
},
{
"family": "Reinert",
"given": "Sandra"
},
{
"family": "Goltstein",
"given": "Pieter M"
},
{
"family": "Boubenec",
"given": "Yves"
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "9904",
"DOI": "10.1038/s41467-026-76104-3",
"PMID": "42749711",
"PMCID": "PMC13582938",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-76104-3",
"language": "en",
"issued": {
"date-parts": [
[
2026,
8,
19
]
]
}
}

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/s41593-026-02333-w [code]
Learning shapes neural geometry in the primate prefrontal cortex.
Journal: Nature neuroscience
In common: seaborn, scikit-learn, pandas, 3 other tools, non-human primate, 9 references
[2] doi:10.1371/journal.pbio.3003831 [code]
Disinhibitory signaling enables flexible coding of top-down information in cortical networks.
Journal: PLoS biology
In common: NetworkX, scikit-learn, pandas, 3 other tools, systems, 6 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, 6 references
[4] 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, 5 references
[5] doi:10.1038/s41467-026-74347-8 [code]
Compositionality of social gaze in the prefrontal-amygdala circuits.
Journal: Nature communications
In common: statsmodels, seaborn, scikit-learn, 4 other tools, non-human primate, systems, 3 references
[6] doi:10.1371/journal.pcbi.1014162 [code]
Exploring neural manifolds across a wide range of intrinsic dimensions.
Journal: PLoS computational biology
In common: scikit-learn, pandas, SciPy, 2 other tools, 5 references
[7] doi:10.1016/j.celrep.2026.117420 [code]
Neural population dynamics of direct electrical stimulation of neocortex.
Journal: Cell reports
In common: NetworkX, statsmodels, seaborn, 5 other tools, systems, 1 reference
[8] doi:10.1038/s41540-026-00727-x [code]
Association-sensory spatiotemporal hierarchy and functional gradient-regularised recurrent neural network with implications for schizophrenia.
Journal: NPJ systems biology and applications
In common: statsmodels, seaborn, pandas, 3 other tools, 3 references
[9] doi:10.7554/elife.109717 [code]
Retrosplenial cortex enables context-dependent goal-directed sensorimotor transformation.
Journal: eLife
In common: statsmodels, seaborn, scikit-learn, 4 other tools, systems, 2 references
[10] doi:10.1038/s41467-026-72057-9 [code]
Sex-specific behavioral feedback modulates sensorimotor processing and drives flexible social behavior.
Journal: Nature communications
In common: statsmodels, seaborn, scikit-learn, 4 other tools, systems, 2 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.