OSCR

Spatially structured heterogeneity shapes large-scale cortical dynamics in a model of the human cortex.

Code ↔ Paper

4 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 4 matches
  1. [1] § Materials and Methods › Cholinergic Muscarinic Regional Heterogeneity. ↔ HeterogeneousWBrainSimulations.ipynb, lines 19–62 · score 0.85 · probe selection, robust sigmoid, ABAGEN toolbox, mirrored, Gene, subtypes
  2. [2] § Results › Framework Overview. ↔ HeterogeneousWBrainSimulations.ipynb, lines 19–62 · score 0.68 · Desikan Killiany, ABAGEN toolbox, adaptation parameter, atlas, Gene, subtypes
  3. [3] § Materials and Methods › Computational Simulations. ↔ HeterogeneousWBrainSimulations.ipynb, lines 383–451 · score 0.56 · duration, Heun, stochastic, transients, stimulation, simulations
  4. [4] § Materials and Methods › Data Analysis. ↔ TE.py, lines 12–114 · score 0.51 · TE implementation, history, uncertainty, delay, entropy

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 · 597 lines · 19 KB · no license · 3 matches

  1. # %% [markdown]
  2. # ### This notebook aims to introduce the MF-Adex framework used in our publication.
  3. #
  4. # If you use this code, please cite:
  5. # Dalla Porta, L. (2026). Spatially structured heterogeneity shapes large-scale cortical dynamics in a model of the human cortex. PNAS, 2026
  6. #
  7. # Contact: [email hidden]
  9. # %%
  10. import matplotlib.pylab as plt
  11. import numpy as np
  12. import pickle
  14. from tvb.simulator.lab import *
  15. from tvb.datatypes.connectivity import Connectivity
  17. import os
  19. # %% [markdown]
  20. # ### I. Load receptors map & normalize
  21. # #### In the paper I used the abagen toolbox following the procedure described in Deco et al. 2011 (https://doi.org/10.1126/sciadv.abf4752)
  22. #
  23. # The original values used in the paper are extracted using the abagen toolbox (https://abagen.readthedocs.io/en/stable/) and are commented below.
  24. # For illustration in this notebok I will use random values.
  26. # %%
  27. # Do it only once and save the ACH_clean for future use
  29. # import abagen
  30. # import pandas as pd
  31. # import numpy as np
  32. # from abagen import images
  34. # files = abagen.fetch_microarray(donors='all', data_dir = 'C:\\Users\\Leonardo\\abagen-data\\microarray\\')
  35. # atlas = abagen.fetch_desikan_killiany()
  36. # expression = abagen.get_expression_data(atlas['image'], gene_norm='robust_sigmoid', probe_selection='rnaseq')
  38. # # Define the gene names, here muscarinic subtype 1 (CHRM1) and subtype 2 (CHRM2) were used
  39. # ACH = expression['CHRM1'] + expression['CHRM2']
  40. # # As in Deco et al. 2011, we use only the left hemisphere as it has more samples. We mirror it to the right hems.
  41. # # Make sure to have only Desikan-killiany Cortex regions, labelled accordingly.
  42. # ACH_clean = np.zeros(68)
  43. # ACH_clean[0:34] = ACH[0:34]
  44. # ACH_clean[34:68] = ACH[0:34] # mirror the left side to right side
  47. # # NORMALIZE between 0 and 1
  48. # max_recp = ACH_clean.max()
  49. # min_recp = ACH_clean.min()
  50. # for i in range(len(ACH_clean)):
  51. # ACH_clean[i] = (ACH_clean[i] - min_recp) / (max_recp - min_recp)
  53. # ACH_clean = 1-ACH_clean # invert sign so that the higher the expression, the lower b (the adaptation parameter)
  54. # # Center the distribution around 1, so that it can be more comparable to the homogenous case
  55. # median_ACH = 1 - np.median(ACH_clean)
  56. # for i in range(len(ACH_clean)):
  57. # ACH_clean[i] = ACH_clean[i] + median_ACH
  59. # ACH_clean is now the normalized muscarinic acetylcholine (CHRM1+CHRM2) density that will be used to modulate the adaptation strenght (b) in the model
  61. # For the purpose of this example I will use a normal distribution of values
  62. ACH_clean = np.random.normal(loc=1, scale=0.1, size=68)
  64. # %% [markdown]
  65. # ### II. Define the structure connectivity (SC) following The Virtual Brain (TVB) guidelines
  66. # #### In the paper a specific SC was used (see details in the paper).
  67. # For the sake of illustration, I will use the cocomac dataset. Download it here: https://zenodo.org/records/14992335
  69. # %%
  70. # If you have your own zip file organized as described in TVB, just do:
  71. # path_connectivity = "C://///"
  72. # conn = connectivity.Connectivity.from_file(
  73. # os.path.abspath(path_connectivity + "Connectivity.zip"),
  74. # )
  75. # conn.weights = conn.weights/(np.sum(conn.weights,axis=0)+1e-12)
  76. # conn.speed = np.r_[4.0]
  78. # As example, let's use the cocomac that we just dowloaded
  79. conn = connectivity.Connectivity.from_file("tvb_data/tvb_data/connectivity/connectivity_68.zip")
  80. conn.weights = conn.weights/(np.sum(conn.weights,axis=0)+1e-12)
  81. conn.speed = np.r_[4.0]
  82. conn.configure()
  84. # Visualize
  85. plt.figure()
  86. plt.title("Cocomac SC 68 regions")
  87. plt.xlabel("ROIs")
  88. plt.ylabel("ROIs")
  89. im = plt.imshow(conn.weights)
  90. cbar = plt.colorbar(im)
  91. cbar.set_label("Weights")
  92. plt.show()
  94. # %% [markdown]
  95. # ### III. Define the dynamic model and its parameters
  96. # #### We use de AdEx Mean-Field Model (TVB implementation) with the given parameters
  97. # Pay attention to the b_e variable, it is where the heterogeneity through ACH_clean is introduced following the Eq. described in the paper
  98. # Shortly, if heterog_fac = 0 we have the homogenous case, and heterog_fac = 1 we have the fully heterogenous case. The strenght of adaptation is regulated by adaptation_ACh variable.
  100. # %%
  101. heterog_fac = 1 # between 0 and 1, set the level of heterogeneity
  102. adaptation_ACh = 30 # it is the amount of adaptation dependent on the heterogeneity factor
  103. # Notice a constant value of 10 in the b_e; it sets the baseline adaptation.
  105. adex = models.ZerlautAdaptationSecondOrder(
  106. g_L = np.r_[10.0],
  107. E_L_e = np.r_[-64.0],
  108. E_L_i = np.r_[-65.0],
  109. C_m = np.r_[200.0],
  110. b_e = (10 + ACH_clean*(heterog_fac*adaptation_ACh) + (adaptation_ACh*(1-heterog_fac))),
  111. a_e = np.r_[0.0],
  112. b_i = np.r_[0.0],
  113. a_i = np.r_[0.0],
  114. tau_w_e = np.r_[500.0],
  115. tau_w_i = np.r_[1.0],
  116. E_e = np.r_[0.0],
  117. E_i = np.r_[-80.0],
  118. Q_e = np.r_[1.5],
  119. Q_i = np.r_[5.0],
  120. tau_e = np.r_[5.0],
  121. tau_i = np.r_[5.0],
  122. N_tot = np.r_[10000],
  123. p_connect_e = np.r_[0.05],
  124. p_connect_i = np.r_[0.05],
  125. g = np.r_[0.2],
  126. T = np.r_[20.0],
  127. P_e = np.r_[
  128. [
  129. -0.05017034,
  130. 0.00451531,
  131. -0.00794377,
  132. -0.00208418,
  133. -0.00054697,
  134. 0.00341614,
  135. -0.01156433,
  136. 0.00194753,
  137. 0.00274079,
  138. -0.01066769,
  139. ]
  140. ],
  141. P_i = np.r_[
  142. [
  143. -0.05184978,
  144. 0.0061593,
  145. -0.01403522,
  146. 0.00166511,
  147. -0.0020559,
  148. 0.00318432,
  149. -0.03112775,
  150. 0.00656668,
  151. 0.00171829,
  152. -0.04516385,
  153. ]
  154. ],
  155. external_input_ex_ex = np.r_[0.315 * 1e-3],
  156. external_input_ex_in = np.r_[0.000],
  157. external_input_in_ex = np.r_[0.315 * 1e-3],
  158. external_input_in_in = np.r_[0.000],
  159. K_ext_e = np.r_[400],
  160. K_ext_i = np.r_[0],
  161. tau_OU = np.r_[2.0],
  162. weight_noise = np.r_[2e-4], #1e-4
  163. )
  164. adex.variables_of_interest = ['E', 'I', 'C_ee', 'C_ei', 'C_ii', 'W_e', 'W_i', 'ou_drift']
  166. adex.state_variable_range["E"] = [0.000, 0.000]
  167. adex.state_variable_range["I"] = [0.00, 0.00]
  168. adex.state_variable_range["C_ee"] = [0.0, 0.0]
  169. adex.state_variable_range["C_ei"] = [0.0, 0.0]
  170. adex.state_variable_range["C_ii"] = [0.0, 0.0]
  171. adex.state_variable_range["W_e"] = [100., 100.0] # 100?
  172. adex.state_variable_range["W_i"] = [0.0, 0.0]
  173. adex.state_variable_range["ou_drift"] = [0.0, 0.0]
  175. # %%
  176. # For each region, the value of adaptation is:
  177. if heterog_fac:
  178. print("You have choosen an Heterogenous case!\n")
  179. else:
  180. print("You have choosen an Homogeneous case\n")
  181. print(f"Your mean adaptation value should be approximately: {10+adaptation_ACh}")
  182. print("The adaptation values for each ROI:")
  183. print(adex.b_e)
  184. print(f"The real mean is: {np.mean(adex.b_e)}")
  186. # %% [markdown]
  187. # ### IVa. Set the simulator
  189. # %%
  190. sim = simulator.Simulator(
  191. model=adex,
  192. connectivity=conn,
  193. conduction_speed=conn.speed.item(),
  194. coupling=coupling.Linear(a = np.r_[0.3], b = np.r_[.0]),
  195. integrator=integrators.HeunStochastic(
  196. noise = noise.Additive(nsig=np.r_[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0], noise_seed = 1234567), #variables ['E', 'I', 'C_ee', 'C_ei', 'C_ii', 'W_e', 'W_i', 'ou_drift']
  197. dt = 0.1,
  198. ),
  199. monitors=[monitors.TemporalAverage(period=1.0)],
  200. ).configure()
  202. # %% [markdown]
  203. # ### IVb. Run the simulation
  205. # %%
  206. transient = 2000 # times are in ms
  207. time_simulation = 6000
  209. (out_t, out_d),= sim.run(simulation_length = time_simulation + transient)
  211. # cut transient
  212. out_d = out_d[transient:,:,:,:]
  213. out_t = out_t[transient:]
  215. # store data
  216. data_ = {}
  217. data_["Time"] = out_t*0.001
  218. data_["TimeSeries_chs"] = out_d[:,0,:,0]*1e3
  220. # %% [markdown]
  221. # ### V. Analysis
  223. # %% [markdown]
  224. # #### Functional connectivity (Pearson Correlation)
  226. # %%
  227. FC = np.corrcoef(np.transpose(data_["TimeSeries_chs"][3800:,:]))
  229. FC_aux = FC.copy()
  230. np.fill_diagonal(FC_aux, 0.0)
  231. iu = np.triu_indices_from(FC_aux, k=1) # Upper diagonal only
  232. mean_FC = np.nanmean(FC_aux[iu])
  234. print(f"The <FC> is {mean_FC}")
  236. fig, axs = plt.subplots(ncols=2)
  238. im1 = axs[0].imshow(conn.weights, cmap = 'RdBu_r', vmin = 0, vmax = 0.5)
  239. im2 = axs[1].imshow(FC, cmap = "RdBu_r", vmin = -0.5, vmax = 1)
  241. axs[0].set_title("Structural Connectivity")
  242. axs[1].set_title("Functional Connectivity")
  244. axs[0].set_ylabel("ROIs")
  245. axs[0].set_xlabel("ROIs")
  246. axs[1].set_xlabel("ROIs")
  248. fig.colorbar(im1, ax=axs[0],fraction=0.046, pad=0.04)
  249. fig.colorbar(im2, ax=axs[1], fraction=0.046, pad=0.04)
  251. plt.tight_layout()
  252. plt.show()
  254. # %% [markdown]
  255. # #### Structure-Function coupling
  256. #
  257. # See paper for details. See also Baum et al 2019 (https://www.pnas.org/doi/10.1073/pnas.1912034117) and Fotiadis et al 2024 (https://www.nature.com/articles/s41583-024-00846-6)
  259. # %%
  260. from netneurotools import metrics # Communicability is already implemented, see documentation (https://netneurotools.readthedocs.io/en/latest/)
  261. from scipy.stats import spearmanr
  263. Q = metrics.communicability_wei(conn.weights)
  264. # compute nodwise spearman(FC,Communicability)
  265. n = Q.shape[0]
  266. c = np.full(n, np.nan)
  267. for ii in range(n):
  268. x = Q[ii, :]
  269. y = FC_aux[ii, :]
  270. m = np.isfinite(x) & np.isfinite(y)
  271. m[ii] = False
  272. c[ii] = spearmanr(x[m], y[m]).correlation
  274. print(f"The correlation between structural connectiivyt and functional dynamic patterns is {np.nanmean(c)}")
  276. # %% [markdown]
  277. # #### Power Spectral Density (PSD)
  279. # %%
  280. from scipy.signal import welch
  282. i,j=8,0
  283. dt = data_["Time"][1]-data_["Time"][0]
  284. fs = 1./dt
  286. aux_nfft = 3
  287. aux_nperseg = 1
  289. sig = data_["TimeSeries_chs"][:,4]
  290. f, S = welch(sig, fs, nfft = aux_nfft*fs, nperseg = aux_nperseg * fs)
  292. plt.figure()
  294. plt.loglog(f, S, lw = 2)
  295. plt.xlabel("Frequency (Hz)", fontsize = 14)
  296. plt.ylabel("PSD", fontsize = 14)
  297. plt.gca().spines['right'].set_visible(False)
  298. plt.gca().spines['top'].set_visible(False)
  299. plt.show()
  301. # %% [markdown]
  302. # #### Phase-Lag Index (PLI)
  303. # See paper for details. See also Stam et al. 2007 (https://doi.org/10.1002/hbm.20346)
  305. # %%
  306. from scipy.signal import hilbert
  307. def compute_phase(signal):
  308. return np.angle(hilbert(signal))
  310. data2use = data_["TimeSeries_chs"][3800:,:]
  311. time_len,n_nodes = np.shape(data2use)
  312. matrix_phase = np.zeros((n_nodes,time_len))
  313. PLI_matrix = np.zeros((n_nodes,n_nodes))
  315. # Create phase Matrix
  316. for i in range(n_nodes):
  317. matrix_phase[:][i] = compute_phase(data2use[:,i])
  319. for i in range(n_nodes):
  320. for j in range(i, n_nodes):
  321. phase_diff = matrix_phase[i] - matrix_phase[j]
  322. PLI_aux = np.abs(np.mean(np.sign(phase_diff)))
  323. PLI_matrix[i, j] = PLI_matrix[j, i] = PLI_aux
  325. np.fill_diagonal(PLI_matrix, np.nan)
  326. meanPLI = np.nanmean(PLI_matrix)
  328. print(f"The mean PLI across all ROIs is {meanPLI}")
  330. # %% [markdown]
  331. # #### Correlations between delta power vs adaptation levels/in-degree connectivity
  333. # %%
  334. # compute PSD for each region
  335. dt = data_["Time"][1]-data_["Time"][0]
  336. fs = 1./dt
  338. f = {}
  339. S = {}
  341. aux_nfft = 10
  342. aux_nperseg = 1
  344. for i in range(np.shape(data_["TimeSeries_chs"])[1]):
  345. sig = data_["TimeSeries_chs"][:,i]
  346. f[i], S[i] = welch(sig, fs, nfft = aux_nfft*fs, nperseg = aux_nperseg * fs)
  349. # Filter delta band power
  350. idx_delta = np.where((f[0]>0.5)&(f[0]<3))
  352. # Extract delta power for each ROI
  353. delta_power = np.zeros(np.shape(data_["TimeSeries_chs"])[1])
  354. for i in range(np.shape(data_["TimeSeries_chs"])[1]):
  355. delta_power[i] = np.mean(S[i][idx_delta])
  357. # Extract in-degree connectivity weights for each ROI
  358. weights_incoming = np.zeros(68)
  359. for i in range(np.shape(data_["TimeSeries_chs"])[1]):
  360. weights_incoming[i] = np.mean(conn.weights[i,:])
  362. # %%
  363. fig, axs = plt.subplots(ncols=3,figsize=(12,4))
  365. im = axs[0].imshow(data_["TimeSeries_chs"][3000:5000,:].T,aspect="auto",cmap="inferno")
  366. cbar = plt.colorbar(im)
  367. cbar.set_label("Firing rate (Hz)")
  368. axs[0].set_xlabel("Time")
  369. axs[0].set_ylabel("ROIs")
  371. axs[1].scatter(np.log(delta_power),np.log(ACH_clean), c = "gray")
  372. axs[1].set_ylabel("Log(Adaptation Level)")
  373. axs[1].set_xlabel("Log(DeltaPower)")
  375. axs[2].set_ylabel("Adaptation Level")
  376. axs[2].scatter(np.log(delta_power),np.log(weights_incoming), c = "gray")
  377. axs[2].set_ylabel("Log(<In-degree>)")
  378. axs[2].set_xlabel("Log(DeltaPower)")
  380. plt.tight_layout()
  381. plt.show()
  383. # %% [markdown]
  384. # ### VI. Evoked activity and Transfer Entropy
  385. #
  386. # As described in the manuscript, to compute (TE) we briefly stimulated a single cortical area and computed TE between the stimulated region and all other areas.
  387. # The TE implementation was borrowed from https://github.com/notsebastiano/transfer_entropy
  389. # %%
  390. # First, we will need to modify the simulator to account for a stimuli
  391. nnodes = len(conn.region_labels)
  392. weight = list(np.zeros(nnodes)) # the stimuli strength, initialized to zero
  394. weight[27] = 1e-4 # defined the strength of stimuli and which region will receive it
  396. parameter_stimulus = {
  397. 'onset': 99.0,
  398. "tau": 9.0,
  399. "T": 99.0,
  400. "weights": None,
  401. "variables":[0]
  402. }
  403. parameter_stimulus['onset']= 3000 #onset time of the stimulus [ms]
  404. parameter_stimulus["tau"]= 30 # stimulus duration [ms]
  405. parameter_stimulus["T"]= 1e9 #interstimulus interval [ms]
  406. parameter_stimulus["weights"]= weight
  407. parameter_stimulus["variables"]=[0] #variable to kick; 0 excitatory, 1 inhibitory, 2 std excitatory, 3 covariation of ex and in, 4 std inhibitory,
  408. #5 adaptation excitatory, 6 adaptation inhibitory
  409. adex.stvar = parameter_stimulus['variables']
  411. stim_time = parameter_stimulus['onset']
  412. stim_steps = stim_time*10 #number of steps until stimulus
  414. eqn_t = equations.PulseTrain()
  415. eqn_t.parameters["onset"] = np.array(parameter_stimulus["onset"]) # ms
  416. eqn_t.parameters["tau"] = np.array(parameter_stimulus["tau"]) # ms
  417. eqn_t.parameters["T"] = np.array(parameter_stimulus["T"]) # ms; # 0.02kHz repetition frequency
  419. # Set the stimulation
  420. stimulation = patterns.StimuliRegion(temporal = eqn_t,
  421. connectivity = conn,
  422. weight=np.array(parameter_stimulus['weights']))
  424. # Introduce the stimulus in the Simulator setup
  425. sim = simulator.Simulator(
  426. model=adex,
  427. connectivity=conn,
  428. conduction_speed=conn.speed.item(),
  429. coupling=coupling.Linear(a = np.r_[0.3], b = np.r_[.0]),
  430. integrator=integrators.HeunStochastic(
  431. noise = noise.Additive(nsig=np.r_[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0], noise_seed = 1234567), #variables ['E', 'I', 'C_ee', 'C_ei', 'C_ii', 'W_e', 'W_i', 'ou_drift']
  432. dt = 0.1,
  433. ),
  434. monitors=[monitors.TemporalAverage(period=1.0)],
  435. stimulus = stimulation, # HERE
  436. ).configure()
  438. # Run the model again
  439. transient = 2000 # times are in ms
  440. time_simulation = 6000
  442. (out_t, out_d),= sim.run(simulation_length = time_simulation + transient)
  444. # cut transient
  445. out_d = out_d[transient:,:,:,:]
  446. out_t = out_t[transient:]
  448. # store data
  449. data_ = {}
  450. data_["Time"] = out_t*0.001
  451. data_["TimeSeries_chs"] = out_d[:,0,:,0]*1e3
  453. # %% [markdown]
  454. # #### Compute TE
  455. #
  456. # Notice that in the example stimuli is applied at 3000ms, and later we discard the first 2000ms. Therefore the stimuli should be centered in 1000ms.
  458. # %%
  459. from TE import transfer_entropy # from https://github.com/notsebastiano/transfer_entropy
  461. id1 = 27 # ID stim
  462. inter_i = 1001 # start period
  463. inter_f = 1401 # end period
  465. TF_values = []
  466. for id2 in range(68):
  467. X = data_["TimeSeries_chs"][inter_i:inter_f,id1]
  468. Y = data_["TimeSeries_chs"][inter_i:inter_f,id2]
  469. TE_YX = transfer_entropy(X, Y, delay=5)
  470. TF_values.append(TE_YX)
  472. TF_values[id1] = np.nan
  473. print(f"The mean TE from source to all brain regions is {np.nanmean(TF_values)}")
  475. # %% [markdown]
  476. # ### VII. Extras
  478. # %% [markdown]
  479. # To compute the correlation between receptor maps spin test was perfomed using the Enigma toolbox (https://enigma-toolbox.readthedocs.io/en/latest/index.html)
  481. # %%
  482. # from scipy.stats import spearmanr, pearsonr
  483. # from enigmatoolbox.permutation_testing import spin_test
  485. # data_1 = np.random.rand(68)
  486. # data_2 = np.random.rand(68)
  488. # # Observed correlation
  489. # rho_obs, _ = pearsonr(data_1, data_1)
  491. # # Spin test (use ENIGMA's built-in DK68 surface parcellation)
  492. # p_spin = spin_test(
  493. # data_1,
  494. # data_2,
  495. # surface_name='fsa5',
  496. # parcellation_name='aparc',
  497. # n_rot=10000,
  498. # type='pearson', # IMPORTANT: match your correlation
  499. # null_dist=False,
  500. # ventricles=False
  501. # )
  503. # print(f"Observed rho = {rho_obs:.3f}")
  504. # print(f"Spin-test p = {p_spin:.6f}")
  506. # %% [markdown]
  507. # To compute the spin rotated maps brainspace was used: https://brainspace.readthedocs.io/en/latest/
  509. # %%
  510. # import numpy as np
  511. # import nibabel.freesurfer.io as fsio
  512. # from brainspace.null_models import SpinPermutations
  514. # FS_HOME = "/Applications/freesurfer/8.0.0"
  515. # N_PERM = 1000 # how many spin rotated maps will be generated
  516. # SEED = 42
  518. # # spheres (fsaverage5)
  519. # sphere_lh, _ = fsio.read_geometry(f"{FS_HOME}/subjects/fsaverage5/surf/lh.sphere")
  520. # sphere_rh, _ = fsio.read_geometry(f"{FS_HOME}/subjects/fsaverage5/surf/rh.sphere")
  522. # # annotations
  523. # labels_lh, _, names_lh_raw = fsio.read_annot("lh.aparc_fsaverage5.annot")
  524. # labels_rh, _, names_rh_raw = fsio.read_annot("rh.aparc_fsaverage5.annot")
  525. # labels_lh = labels_lh.astype(int)
  526. # labels_rh = labels_rh.astype(int)
  528. # names_lh = [n.decode("utf-8") if isinstance(n, bytes) else str(n) for n in names_lh_raw]
  529. # names_rh = [n.decode("utf-8") if isinstance(n, bytes) else str(n) for n in names_rh_raw]
  531. # # pick cortical parcel IDs: all labels except "unknown" and "corpuscallosum"
  532. # EXCLUDE = {"unknown", "corpuscallosum"}
  534. # lh_ids = [i for i, n in enumerate(names_lh) if n not in EXCLUDE]
  535. # rh_ids = [i for i, n in enumerate(names_rh) if n not in EXCLUDE]
  537. # # Keep only parcel IDs that actually appear in the vertex labels (avoids empty slices)
  538. # lh_ids = [i for i in lh_ids if np.any(labels_lh == i)]
  539. # rh_ids = [i for i in rh_ids if np.any(labels_rh == i)]
  541. # lh_ids = lh_ids[:34]
  542. # rh_ids = rh_ids[:34]
  544. # parcel_names = [names_lh[i] for i in lh_ids] + [names_rh[i] for i in rh_ids]
  547. # # fit spins
  548. # sp = SpinPermutations(n_rep=N_PERM, random_state=SEED)
  549. # sp.fit(sphere_lh, sphere_rh)
  551. # # load RNA data
  552. # chrm = ACH_clean
  554. # # build vertex maps (per hemisphere)
  555. # vmap_lh = np.zeros(labels_lh.shape, dtype=float)
  556. # vmap_rh = np.zeros(labels_rh.shape, dtype=float)
  558. # for j, pid in enumerate(lh_ids):
  559. # vmap_lh[labels_lh == pid] = chrm[j]
  560. # for j, pid in enumerate(rh_ids):
  561. # vmap_rh[labels_rh == pid] = chrm[34 + j]
  563. # # spin-randomize
  564. # vmap_spins_lh, vmap_spins_rh = sp.randomize(vmap_lh, vmap_rh)
  566. # # average back to DK68
  567. # sp_maps = np.zeros((N_PERM, 68), dtype=float)
  568. # for i in range(N_PERM):
  569. # for j, pid in enumerate(lh_ids):
  570. # sp_maps[i, j] = vmap_spins_lh[i][labels_lh == pid].mean()
  571. # for j, pid in enumerate(rh_ids):
  572. # sp_maps[i, 34 + j] = vmap_spins_rh[i][labels_rh == pid].mean()
  574. # # Spin rotated maps are stored in sp_maps
  576. # %% [markdown]
  577. # Statistics
  579. # %%
  580. # # z-values (z-score)
  581. # # aligned = mean value of the observable
  582. # # mu = mean value of the null model
  583. # # sd = standard deviation of the null model
  584. # delta = aligned - mu
  585. # z = delta / sd if sd > 0 else np.nan
  587. # # Empirical p-value
  588. # # observed = mean value of the observable. It's a single value
  589. # # null = distribution of mean values of the null models
  590. # # z value relative to null
  591. # null_mean = null.mean()
  592. # null_std = null.std(ddof=0)
  593. # # two-sided empirical p value
  594. # obs_dev = abs(observed - null_mean)
  595. # null_dev = abs(null - null_mean)
  596. # p_emp = (np.sum(null_dev >= obs_dev) + 1) / (len(null) + 1)
  597. # print(f"p_emp = {p_emp:.4g}")

HeterogeneousWBrainSimulations.ipynb at commit d458969, no license · at the source

Overview

Authors: Leonardo Dalla Porta1, Jan Fousek2, Alain Destexhe3, Maria V Sanchez-Vives1,4
ORCID iDs: Leonardo Dalla Porta, Jan Fousek, Alain Destexhe, Maria V Sanchez-Vives
  1. Institute of Biomedical Investigations August Pi i Sunyer, Systems Neuroscience, Barcelona 08036, Spain
  2. Central European Institute of Technology, Masaryk University, Brno 65691, Czech Republic
  3. Department for Integrative and Computational Neuroscience, Paris-Saclay University, CNRS, Paris-Saclay Institute of Neuroscience, Saclay 91400, France
  4. Catalan Institution for Research and Advanced Studies, Barcelona 08010, Spain
Institutions: Consorci Institut D'Investigacions Biomediques August Pi I Sunyer (Spain); Central European Institute of Technology – Masaryk University (Czechia); Centre National de la Recherche Scientifique (France); Université Paris-Saclay (France); Institut des Neurosciences Paris-Saclay (France); Institució Catalana de Recerca i Estudis Avançats (Spain)
Journal: Proceedings of the National Academy of Sciences of the United States of America, volume 123, issue 28, article e2532072123
Dates: received 6 November 2025; accepted 10 June 2026; published online 8 July 2026; in print 14 July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1073/pnas.2532072123 · PMID 42418484 · PMCID PMC13367858 · OpenAlex W7167718845
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: PET / SPECT (modality), human (organism), cellular / molecular (subfield)
Methods: Spectral & time-frequency, Statistics, Connectivity, Single-unit activity, calcium imaging, fMRI & imaging
Keywords: hierarchical heterogeneity, brain states, neuromodulation, muscarinic, whole-brain modeling
MeSH: Cerebral Cortex*, Models, Neurological*, Brain Mapping, Humans, Positron-Emission Tomography, Receptors, Muscarinic, Sleep, Wakefulness (* major topic)
Topic: Neural dynamics and brain function (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) (444500/2024-3); EC | Horizon 2020 Framework Programme (945539); Ministerio de Ciencia, Innovación y Universidades (MCIU) (PID2023-152918OB-I00); Government of Catalonia | Agència de Gestió d'Ajuts Universitaris i de Recerca (AGAUR) (2021-SGR-01165); European Research Council (101071900)
Citations: not cited yet (Europe PMC); 115 references in the paper

Abstract

The human brain displays substantial spatial variability in molecular, anatomical, and physiological organization. Yet, how this heterogeneity shapes large-scale neuronal dynamics remains poorly understood. To address this question, we employed a biologically informed large-scale cortical model capable of generating distinct brain states, from awake-like to sleep-like dynamics. Our model was constrained by empirical human structural connectivity (SC) and regional cholinergic muscarinic receptor (CHRM) maps derived from transcriptomic data, together with complementary positron emission tomography (PET)-based receptor maps. These regional maps were implemented as modulators of adaptation-related excitability. We found that modulating excitability according to the spatial organization of CHRM maps significantly impacted large-scale cortical dynamics: It not only facilitated network synchronization but also enhanced information flow between cortical regions. Importantly, these effects were conserved across transcriptomic and PET-derived maps and could not be fully reproduced by multiple null models preserving generic forms of heterogeneity. Moreover, we addressed a particularly intricate dynamic regime characterized by the coexistence of localized sleep-like activity within otherwise awake-like states. We showed that the emergence of these sleep-like slow waves was a byproduct of both regional levels of neuronal adaptation and SC. In summary, our findings highlight the critical role of molecular and anatomical heterogeneity in shaping widespread cortical dynamics, suggesting broad avenues for linking microscale diversity to macroscale function.

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

notsebastiano/transfer_entropy

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 165baefd5621c131aaf07c422bf8cdbf66ca0be2, 28 February 2025
Languages: Python (1), Jupyter (1)
Size: 7 files, 2 scripts
Software Heritage: archived
Found in: the text, “Data Analysis.”
Holds: README, license file, 1 notebook
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (2 files), SciPy (2 files), Matplotlib (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
4 files
At the source: github.com/notsebastiano/transfer_entropy

ldallap/HeterogeneousBrainModel

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: d45896933c6de91b941815625cc012cb688be6d0, 9 July 2026
Languages: Jupyter (1)
Size: 2 files, 1 script
Software Heritage: not archived
Found in: “Data, Materials, and Software Availability”
Holds: README, 1 notebook
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Matplotlib (1 file), netneurotools (1 file), NumPy (1 file), SciPy (1 file), The Virtual Brain (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
2 files
At the source: github.com/ldallap/HeterogeneousBrainModel

The paper's code and data availability statement is in the Data section.

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;
  • 3 scripts, each with its path and the digest of its content;
  • 4 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

No dataset and no data link were found in the paper.

Data, Materials, and Software Availability

The code needed to replicate the main findings of this study has been deposited in Heterogeneous Brain Model (https://github.com/ldallap/HeterogeneousBrainModel) (115). All other data are included in the manuscript and/or SI Appendix (http://www.pnas.org/lookup/doi/10.1073/pnas.2532072123#supplementary-materials).

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

Cite

This paper

Dalla Porta, L., Fousek, J., Destexhe, A., & Sanchez-Vives, M. V. (2026). Spatially structured heterogeneity shapes large-scale cortical dynamics in a model of the human cortex. Proceedings of the National Academy of Sciences of the United States of America, 123(28), e2532072123. https://doi.org/10.1073/pnas.2532072123

BibTeX

@article{dallaporta2026spatially,
author = {Dalla Porta, Leonardo and Fousek, Jan and Destexhe, Alain and Sanchez-Vives, Maria V},
title = {{Spatially structured heterogeneity shapes large-scale cortical dynamics in a model of the human cortex}},
journal = {Proceedings of the National Academy of Sciences of the United States of America},
year = {2026},
month = jul,
volume = {123},
number = {28},
pages = {e2532072123},
publisher = {National Academy of Sciences},
issn = {0027-8424},
doi = {10.1073/pnas.2532072123},
url = {https://doi.org/10.1073/pnas.2532072123},
pmid = {42418484},
pmcid = {PMC13367858}
}

RIS

TY - JOUR
AU - Dalla Porta, Leonardo
AU - Fousek, Jan
AU - Destexhe, Alain
AU - Sanchez-Vives, Maria V
TI - Spatially structured heterogeneity shapes large-scale cortical dynamics in a model of the human cortex
T2 - Proceedings of the National Academy of Sciences of the United States of America
J2 - Proc Natl Acad Sci U S A
PY - 2026
DA - 2026/07/08
VL - 123
IS - 28
SP - e2532072123
SN - 0027-8424
PB - National Academy of Sciences
DO - 10.1073/pnas.2532072123
UR - https://doi.org/10.1073/pnas.2532072123
LA - en
ER -

CSL-JSON

{
"id": "10.1073/pnas.2532072123",
"type": "article-journal",
"title": "Spatially structured heterogeneity shapes large-scale cortical dynamics in a model of the human cortex",
"container-title": "Proceedings of the National Academy of Sciences of the United States of America",
"author": [
{
"family": "Dalla Porta",
"given": "Leonardo"
},
{
"family": "Fousek",
"given": "Jan"
},
{
"family": "Destexhe",
"given": "Alain"
},
{
"family": "Sanchez-Vives",
"given": "Maria V"
}
],
"container-title-short": "Proc Natl Acad Sci U S A",
"volume": "123",
"issue": "28",
"page": "e2532072123",
"DOI": "10.1073/pnas.2532072123",
"PMID": "42418484",
"PMCID": "PMC13367858",
"ISSN": "0027-8424",
"publisher": "National Academy of Sciences",
"URL": "https://doi.org/10.1073/pnas.2532072123",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
8
]
]
}
}

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.1371/journal.pbio.3003916 [code]
Arousal-driven critical roaming reproduces human functional connectivity dynamics.
Journal: PLoS biology
In common: SciPy, Matplotlib, NumPy, 10 references
[2] doi:10.1371/journal.pcbi.1013463 [code]
A multi-frequency whole-brain neural mass model with homeostatic feedback inhibition.
Journal: PLoS computational biology
In common: SciPy, Matplotlib, NumPy, 8 references
[3] doi:10.1186/s12916-026-04903-y [code]
Structural connectome architecture and biological vulnerability shape cortical atrophy in cocaine use disorder.
Journal: BMC medicine
In common: netneurotools, SciPy, Matplotlib, 1 other tool, 7 references
[4] doi:10.1093/sleepadvances/zpag065 [code]
Brain-wide properties of slow waves across vigilance states.
Journal: Sleep advances : a journal of the Sleep Research Society
In common: 8 references
[5] doi:10.1038/s42003-025-09444-3 [code]
Decoupling of neurophysiological activity from structure mirrors global microarchitectural and neuromodulatory trends.
Journal: Communications biology
In common: netneurotools, SciPy, Matplotlib, 1 other tool, cellular / molecular, 6 references
[6] doi:10.1126/sciadv.aef2894 [code]
Human cortical networks trade communication efficiency for computational reliability.
Journal: Science advances
In common: netneurotools, SciPy, Matplotlib, 1 other tool, 5 references
[7] doi:10.1016/j.celrep.2026.117782 [code]
Thermodynamics of consciousness: Non-equilibrium brain dynamics track conscious states.
Journal: Cell reports
In common: SciPy, Matplotlib, NumPy, 6 references
[8] doi:10.1002/hbm.70605 [code]
BrainEnrich: Revealing Biological Insights for Imaging-Derived Phenotypes Through Transcriptomic Enrichment.
Journal: Human brain mapping
In common: netneurotools, SciPy, Matplotlib, 1 other tool, 6 references
[9] doi:10.3390/brainsci16080822
Photopharmacological Cholinergic Modulation of Cortical Activity in Human Brain Slices.
Journal: Brain sciences
In common: 4 references, author Maria V Sanchez-Vives
[10] doi:10.1371/journal.pcbi.1013252 [code]
Dynamic cholinergic signaling differentially desynchronizes cortical microcircuits dependent on modulation rate and network connectivity.
Journal: PLoS computational biology
In common: Matplotlib, NumPy, cellular / molecular, 6 references

Contribute

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

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

Request its removal

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

Discussion, reproductions, activity

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

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

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