OSCR

Charting higher-order models of brain function beyond pairwise interactions.

Code ↔ Paper

23 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 23 matches
  1. [1] § Methods › Data sources and preprocessing › Brain maps from neuromaps ↔ Paper/Figure3/Fig3_clean.ipynb, lines 114–141 · score 0.97 · HT1a, HT1b, HT2a, VAChT, mGluR5, CB1
  2. [2] § Results › Higher-order taxonomy ↔ Paper/Figure2/Fig2_analyses_clean.ipynb, lines 812–865 · score 0.69 · classical FC, TScaffold, structural connectivity, PED Syn, PED Red, PhiID
  3. [3] § Results › Higher-order taxonomy ↔ Paper/Figure2/Fig2_analyses_clean.ipynb, lines 556–651 · score 0.65 · complete linkage, TScaffold, PED Syn, PED Red, PhiID, Pscaff
  4. [4] § Methods › Data sources and preprocessing › Brain parcellation ↔ 01_example_surface_plots.ipynb, lines 99–104 · score 0.61 · dorsal attention, ventral attention, limbic, somatomotor, frontoparietal, parcellation
  5. [5] § Methods › Higher-order frameworks › Integrated information decomposition (Phi ID): PhiID Red and PhiID Syn ↔ hoi/metrics/phiid_red.py, lines 27–59 · score 0.60 · Minimum Mutual Information, redundancy atom, MMI, PhiID
  6. [6] § Methods › Higher-order frameworks › Integrated information decomposition (Phi ID): PhiID Red and PhiID Syn ↔ hoi/metrics/phiid_atoms.py, lines 105–142 · score 0.60 · Minimum Mutual Information, redundancy atom, MMI, PhiID
  7. [7] § Methods › Clustering ↔ Paper/Figure2/Fig2_analyses_clean.ipynb, lines 531–553 · score 0.59 · Spearman correlation, hierarchical clustering, complete linkage, dendrogram, nodal, distance
  8. [8] § Methods › Data sources and preprocessing › Brain parcellation ↔ Paper/Figure3/Fig3_clean.ipynb, lines 175–216 · score 0.58 · subcortical regions, DA, VIS, SM, atlas, VA
  9. [9] § Results › Multimodal characterization of topological and informational gradients ↔ Paper/Figure2/Fig2_analyses_clean.ipynb, lines 699–743 · score 0.58 · TScaffold, scaffold frequency, PED Syn, PED Red, PhiID, profiles
  10. [10] § Results › FC decomposition and brain-behavior analysis ↔ Paper/Figure5-6/Fig5_FC_reconstruction.ipynb, lines 245–280 · score 0.57 · Euclidean distance, Task flexibility, centroid, positions, space, synergistic
  11. [11] § Results › Multimodal characterization of topological and informational gradients ↔ Paper/Figure3/Fig3_clean.ipynb, lines 684–734 · score 0.56 · HOI axis, axis correlates, Box, whiskers, Scatter, transporters
  12. [12] § Methods › Higher-order frameworks › Integrated information decomposition (Phi ID): PhiID Red and PhiID Syn ↔ hoi/metrics/phiid_red.py, lines 27–59 · score 0.55 · redundancy atom, mutual information, variables, PhiID
  13. [13] § Methods › Higher-order frameworks › Integrated information decomposition (Phi ID): PhiID Red and PhiID Syn ↔ hoi/metrics/phiid_atoms.py, lines 105–142 · score 0.55 · redundancy atom, mutual information, variables, PhiID
  14. [14] § Results › Brain fingerprinting and task decoding ↔ Paper/Figure2/Fig2_analyses_clean.ipynb, lines 345–418 · score 0.54 · Yeo resting state, resting state network, Box, DA, VIS, SM
  15. [15] § Results › Brain fingerprinting and task decoding ↔ Paper/Figure4/Fig4_tasks.ipynb, lines 40–86 · score 0.54 · Task decoding, PhiID, DA, VIS, SM, VA
  16. [16] § Methods › Higher-order frameworks › Triangles and temporal scaffold ↔ High_order_TS_with_scaffold/persistent_homology_calculation.py, lines 117–156 · score 0.53 · persistent homology, edge weight, cycles, filtration, dimensional, scaffold
  17. [17] § Methods › Higher-order frameworks › Triangles and temporal scaffold ↔ High_order_TS_with_scaffold/simplicial_multivariate.py, lines 36–95 · score 0.53 · simplicial complex, edge weight, violation, violating, filtration, triangles
  18. [18] § Results ↔ Paper/Figure2/Fig2_analyses_clean.ipynb, lines 1–73 · score 0.53 · structural functional cartography, Partial Entropy Decomposition, TScaffold, homological scaffolds, fMRI, PhiID
  19. [19] § Methods › Higher-order frameworks › Partial entropy decomposition (PED): PED Red and PED Syn ↔ Code/04_PED/ped_module.py, lines 36–64 · score 0.52 · Partial Entropy Decomposition, joint entropy, discrete, PED, Syn
  20. [20] § Methods › Higher-order frameworks › Triangles and temporal scaffold ↔ High_order_TS/simplicial_multivariate.py, lines 33–86 · score 0.52 · simplicial complex, edge weight, violation, violating, filtration, triangles
  21. [21] § Methods › Data sources and preprocessing › Brain maps from neuromaps ↔ utils_neuromaps_brain.py, lines 53–202 · score 0.51 · cortical surface, brain maps, neuromaps, Schaefer, parcellated
  22. [22] § Results › FC decomposition and brain-behavior analysis ↔ Paper/Figure5-6/Fig5_FC_reconstruction.ipynb, lines 325–372 · score 0.51 · dopamine transporter, task flexibility, DAT, correlated, FC
  23. [23] § Methods › Higher-order frameworks › Triangles and temporal scaffold ↔ High_order_TS_with_scaffold/simplicial_multivariate.py, lines 36–95 · score 0.51 · persistent homology, edge weight, filtration, dimensional, scaffold, triangle

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 · 932 lines · 38 KB · no license · 6 matches

  1. # %% [markdown]
  2. # # Figure 2 — Higher-order information dynamics in resting-state fMRI
  3. #
  4. # This notebook reproduces Figure 2, which characterises eleven higher-order
  5. # functional connectivity (HOI) measures computed on HCP resting-state fMRI data
  6. # (Schaefer 100-region cortical parcellation + 16 subcortical regions, N = 100 subjects).
  7. #
  8. # **Measures compared**
  9. #
  10. # | Label | Method |
  11. # |---|---|
  12. # | O-Red / O-Syn | O-information redundancy / synergy |
  13. # | PED Red / PED Syn | Partial Entropy Decomposition |
  14. # | PhiID Red / PhiID Syn | Integrated Information Decomposition |
  15. # | Triangles | Local TDA triangle count |
  16. # | TScaffold | Topological scaffold (local TDA) |
  17. # | Fscaff / Pscaff | Homological scaffold – frequency / persistence |
  18. # | FC | Classical pairwise functional connectivity |
  19. #
  20. # **Figure panels produced**
  21. #
  22. # - **(a)** Edge-level connectivity matrices with overlaid brain maps
  23. # - **(b)** Edge-level pairwise Spearman similarity matrix
  24. # - **(c)** Hierarchical clustering dendrogram
  25. # - **(d)** Structure–function cartography scatter plot
  26. #
  27. # > Set the paths in **Section 1 – Configuration** before running.
  29. # %% [markdown]
  30. # ## 1. Imports
  32. # %%
  33. %load_ext autoreload
  34. %autoreload 2
  36. import numpy as np
  37. import h5py
  38. import sys
  39. import glob
  40. import os
  41. import re
  42. import scipy.io as sio
  43. import pandas as pd
  44. import matplotlib as mpl
  45. import matplotlib.pyplot as plt
  46. import matplotlib.colors as mcolors
  47. import seaborn as sns
  49. from collections import OrderedDict
  50. from scipy.spatial.distance import squareform
  51. from scipy.stats import spearmanr, rankdata
  52. from scipy.cluster.hierarchy import dendrogram, linkage
  53. from scipy.cluster.hierarchy import leaves_list as hc_leaves_list
  54. from matplotlib.gridspec import GridSpec
  56. sys.path.append("../utils/")
  57. from utils import (
  58. extract_red_syn_from_Oinfo_optimized,
  59. triangle_projection,
  60. edge_projection,
  61. upper_tri_masking,
  62. )
  63. from utils_neuromaps_brain import normal_view_top
  65. # Typography — Helvetica preferred (Arial fallback), sans-serif math, embedded PDF fonts
  66. mpl.rcParams.update({
  67. "font.family": "sans-serif",
  68. "font.sans-serif": ["Helvetica", "Arial", "DejaVu Sans"],
  69. "mathtext.fontset": "stixsans", # sans-serif math glyphs
  70. "axes.linewidth": 1.5,
  71. "pdf.fonttype": 42, # embed as TrueType in PDF (required by journals)
  72. "ps.fonttype": 42,
  73. })
  76. # %% [markdown]
  77. # ## 2. Configuration — set paths here
  79. # %%
  80. # ── Paths to pre-computed HOI results ────────────────────────────────────────
  81. BASE_PATH = "../../../"
  82. DATA_PATH = f"{BASE_PATH}Results/08_HCP_Diano_preprocessing/REST_subc/"
  83. FC_PATH = f"{BASE_PATH}Dataset/08_HCP_Diano_preprocessing/REST_reordered_SC/"
  85. # Structural connectivity CSV files (copied to data/SC/ — see README)
  86. SC_PATH = "./data/SC/"
  88. # Output folder for figures
  89. FIG_PATH = "./Figures/"
  90. os.makedirs(FIG_PATH, exist_ok=True)
  92. # ── Atlas parameters ──────────────────────────────────────────────────────────
  93. N_ROIS = 116 # 100 cortical (Schaefer) + 16 subcortical
  94. N_SUBC = 16 # subcortical ROIs appended after cortical ones
  96. # ── Colour-map used for all edge matrices / brain maps ────────────────────────
  97. _n, _c = 50, 0.1
  98. _arr = (1 - _c) * plt.get_cmap("BuPu")(np.linspace(0, 1, _n)) + _c * np.ones((_n, 4))
  99. CMAP_STANDARD = mcolors.ListedColormap(_arr)
  101. # ── Measure display names (order determines panels b, c ordering) ─────────────
  102. MEASURES = [
  103. ("classical_FC_proj", "FC"),
  104. ("triangles_proj", "Triangles"),
  105. ("Red", "O-Red"),
  106. ("PED_Red_proj", "PED Red"),
  107. ("phiIDRed_proj", "PhiID Red"),
  108. ("scaffold_freq_proj","Fscaff"),
  109. ("scaffold_pers_proj","Pscaff"),
  110. ("scaffold_proj", "TScaffold"),
  111. ("Syn", "O-Syn"),
  112. ("PED_Syn_proj", "PED Syn"),
  113. ("phiIDSyn_proj", "PhiID Syn"),
  114. ]
  115. LABELS = [lbl for _, lbl in MEASURES]
  117. # ── Plotting style (markers, colours) for structure-function scatter ──────────
  118. MARKER_STYLE = {
  119. "O-Red": ("o", "#E41A1C", 120),
  120. "O-Syn": ("s", "#377EB8", 100),
  121. "Triangles":("^", "#4DAF4A", 140),
  122. "TScaffold":("D", "#984EA3", 100),
  123. "Fscaff": ("v", "#FF7F00", 140),
  124. "Pscaff": ("<", "#FFFF33", 140),
  125. "PED Red": (">", "#A65628", 140),
  126. "PED Syn": ("p", "#F781BF", 150),
  127. "PhiID Syn":("h", "#999999", 150),
  128. "PhiID Red":("8", "#66C2A5", 150),
  129. "FC": ("*", "#FC8D62", 160),
  130. }
  133. # %% [markdown]
  134. # ## 3. Data loading
  136. # %%
  137. def load_and_project_data(proj="nodal",
  138. data_path=DATA_PATH,
  139. fc_path=FC_PATH,
  140. n_rois=N_ROIS,
  141. n_subc=N_SUBC):
  142. """Load and project all HOI measures for one projection type.
  144. Parameters
  145. ----------
  146. proj : 'nodal' | 'edges'
  147. 'nodal' → each measure is a vector of length n_rois (node strength).
  148. 'edges' → each measure is a condensed upper-triangle vector of length
  149. C(n_rois, 2), suitable for squareform().
  150. data_path : str
  151. Directory with sub-folders 01_Red_Syn/, 02_local_TDA/, etc.
  152. fc_path : str
  153. Directory of per-subject time-series CSV files.
  154. n_rois, n_subc : int
  155. Total ROIs and number of subcortical ROIs appended at the end.
  157. Returns
  158. -------
  159. dict keys → measure name, values → ndarray (n_subjects, feature_size)
  160. """
  161. n_cortical = n_rois - n_subc
  163. # Pre-compute triplet → edge index mapping (needed for edge projection)
  164. from itertools import combinations as _comb
  165. all_comb = np.array(list(_comb(np.arange(n_rois), 3)))
  166. if proj == "edges":
  167. indices = np.array([
  168. np.intersect1d(np.where(all_comb == i)[0],
  169. np.where(all_comb == j)[0])
  170. for i, j in _comb(np.arange(n_rois), 2)
  171. ])
  172. else:
  173. indices = None # triangle_projection computes its own nodal indices
  175. print(f"Loading data (proj='{proj}') …")
  177. # ── 1. O-information: Redundancy and Synergy ─────────────────────────────
  178. Oinfo = {}
  179. with h5py.File(f"{data_path}01_Red_Syn/Oinfo.hdf5", "r") as f:
  180. for key in f:
  181. Oinfo[key] = np.ravel(f[key][:])
  183. Red_d, Syn_d = extract_red_syn_from_Oinfo_optimized(
  184. Oinfo, nROIs=n_rois, proj=proj, normalized=False, indices=indices
  185. )
  186. Red = np.array(list(Red_d.values()))
  187. Syn = np.array(list(Syn_d.values()))
  188. print(" ✓ O-information (Red / Syn)")
  190. # ── 2. Local TDA (scaffold and triangles) ────────────────────────────────
  191. def _load_npz_dict(pattern, key):
  192. result = OrderedDict()
  193. for path in sorted(glob.glob(pattern)):
  194. fname = os.path.basename(path).split(".")[0]
  195. day = (re.search(r"REST(\d+)", fname) or type("", (), {"group": lambda s, _: ""})()).group(1)
  196. lr = (re.search(r"(LR|RL)", fname) or type("", (), {"group": lambda s, _: ""})()).group(1)
  197. subj = (re.search(r"(\d{6})", fname) or type("", (), {"group": lambda s, _: ""})()).group(1)
  198. sid = f"{subj}_{lr}_{day}"
  199. result[sid] = np.load(path)[key]
  200. return result
  202. scaffold_d = _load_npz_dict(f"{data_path}02_local_TDA/*scaffold*", "scaffold")
  203. triangles_d = _load_npz_dict(f"{data_path}02_local_TDA/*triangles*", "triangles")
  205. scaffold_proj = np.array([edge_projection(scaffold_d[s], proj=proj, nROIs=n_rois) for s in scaffold_d])
  206. triangles_proj = np.array([triangle_projection(triangles_d[s], proj=proj, nROIs=n_rois, indices=indices) for s in triangles_d])
  207. print(" ✓ Local TDA (scaffold / triangles)")
  209. # ── 3. Homological scaffolds (frequency and persistence) ─────────────────
  210. def _load_npy_dir(folder):
  211. result = OrderedDict()
  212. for fname in sorted(os.listdir(folder)):
  213. if not fname.startswith("."):
  214. result[fname] = np.load(folder + fname)
  215. return result
  217. freq_d = _load_npy_dir(f"{data_path}03_scaffold/freq/")
  218. pers_d = _load_npy_dir(f"{data_path}03_scaffold/pers/")
  220. scaffold_freq_proj = np.array([edge_projection(freq_d[s], proj=proj, nROIs=n_rois) for s in freq_d])
  221. scaffold_pers_proj = np.array([edge_projection(pers_d[s], proj=proj, nROIs=n_rois) for s in pers_d])
  222. print(" ✓ Homological scaffolds (freq / pers)")
  224. # ── 4. PED ───────────────────────────────────────────────────────────────
  225. PED_Red_d, PED_Syn_d = {}, {}
  226. with h5py.File(f"{data_path}04_PED/Red.hdf5", "r") as f:
  227. for key in f:
  228. PED_Red_d[key] = f[key][:]
  229. with h5py.File(f"{data_path}04_PED/Syn.hdf5", "r") as f:
  230. for key in f:
  231. PED_Syn_d[key] = f[key][:]
  233. PED_Red_proj = np.array([triangle_projection(PED_Red_d[s], proj=proj, nROIs=n_rois, indices=indices) for s in PED_Red_d])
  234. PED_Syn_proj = np.array([triangle_projection(PED_Syn_d[s], proj=proj, nROIs=n_rois, indices=indices) for s in PED_Syn_d])
  235. print(" ✓ PED (Red / Syn)")
  237. # ── 5. PhiID ─────────────────────────────────────────────────────────────
  238. phiID_Red_d, phiID_Syn_d = {}, {}
  239. with h5py.File(f"{data_path}05_phiID/Red.hdf5", "r") as f:
  240. for key in f:
  241. phiID_Red_d[key] = f[key][:]
  242. with h5py.File(f"{data_path}05_phiID/Syn.hdf5", "r") as f:
  243. for key in f:
  244. phiID_Syn_d[key] = f[key][:]
  246. # PhiID is stored as a list of (value, …) tuples — extract first element
  247. if proj == "edges":
  248. def _phiID_to_edge(d):
  249. mat = squareform(np.array([row[0] for row in d]))[:n_rois, :n_rois]
  250. return edge_projection(mat, proj=proj, nROIs=n_rois)
  251. phiIDRed_proj = np.array([_phiID_to_edge(phiID_Red_d[s]) for s in phiID_Red_d])
  252. phiIDSyn_proj = np.array([_phiID_to_edge(phiID_Syn_d[s]) for s in phiID_Syn_d])
  253. else:
  254. phiIDRed_proj = np.array([np.mean(squareform(np.array([row[0] for row in phiID_Red_d[s]])), axis=0)
  255. for s in phiID_Red_d])
  256. phiIDSyn_proj = np.array([np.mean(squareform(np.array([row[0] for row in phiID_Syn_d[s]])), axis=0)
  257. for s in phiID_Syn_d])
  258. print(" ✓ PhiID (Red / Syn)")
  260. # ── 6. Classical FC (from time series) ───────────────────────────────────
  261. FC_d = OrderedDict()
  262. for fname in sorted(os.listdir(fc_path)):
  263. token = fname.split(".")[0]
  264. day = (re.search(r"REST(\d+)", token) or type("", (), {"group": lambda s, _: ""})()).group(1)
  265. lr = (re.search(r"(LR|RL)", token) or type("", (), {"group": lambda s, _: ""})()).group(1)
  266. sid = (re.search(r"(\d{6})", token) or type("", (), {"group": lambda s, _: ""})()).group(1)
  267. sid = f"{sid}_{lr}_{day}"
  268. FC_d[sid] = np.corrcoef(np.loadtxt(fc_path + fname))
  270. classical_FC_proj = np.array([edge_projection(FC_d[s], proj=proj, nROIs=n_rois) for s in FC_d])
  271. print(" ✓ Classical FC")
  272. print(f" Done — {len(Red)} scans loaded.")
  274. return {
  275. "Red": Red,
  276. "Syn": Syn,
  277. "triangles_proj": triangles_proj,
  278. "scaffold_proj": scaffold_proj,
  279. "scaffold_freq_proj":scaffold_freq_proj,
  280. "scaffold_pers_proj":scaffold_pers_proj,
  281. "PED_Red_proj": PED_Red_proj,
  282. "PED_Syn_proj": PED_Syn_proj,
  283. "phiIDRed_proj": phiIDRed_proj,
  284. "phiIDSyn_proj": phiIDSyn_proj,
  285. "classical_FC_proj": classical_FC_proj,
  286. }
  289. # %%
  290. # Nodal projection (used for brain maps and Yeo distributions)
  291. data_nodal = load_and_project_data(proj="nodal")
  293. Red_nodal = data_nodal["Red"]
  294. Syn_nodal = data_nodal["Syn"]
  295. triangles_nodal = data_nodal["triangles_proj"]
  296. scaffold_nodal = data_nodal["scaffold_proj"]
  297. scaffold_freq_nodal= data_nodal["scaffold_freq_proj"]
  298. scaffold_pers_nodal= data_nodal["scaffold_pers_proj"]
  299. PED_Red_nodal = data_nodal["PED_Red_proj"]
  300. PED_Syn_nodal = data_nodal["PED_Syn_proj"]
  301. phiIDRed_nodal = data_nodal["phiIDRed_proj"]
  302. phiIDSyn_nodal = data_nodal["phiIDSyn_proj"]
  303. classical_FC_nodal = data_nodal["classical_FC_proj"]
  306. # %%
  307. # Edge projection (used for similarity analysis and structure-function plots)
  308. data_edges = load_and_project_data(proj="edges")
  310. Red_edges = data_edges["Red"]
  311. Syn_edges = data_edges["Syn"]
  312. triangles_edges = data_edges["triangles_proj"]
  313. scaffold_edges = data_edges["scaffold_proj"]
  314. scaffold_freq_edges = data_edges["scaffold_freq_proj"]
  315. scaffold_pers_edges = data_edges["scaffold_pers_proj"]
  316. PED_Red_edges = data_edges["PED_Red_proj"]
  317. PED_Syn_edges = data_edges["PED_Syn_proj"]
  318. phiIDRed_edges = data_edges["phiIDRed_proj"]
  319. phiIDSyn_edges = data_edges["phiIDSyn_proj"]
  320. classical_FC_edges = data_edges["classical_FC_proj"]
  323. # %%
  324. # Load structural connectivity matrices (100 subjects, 116×116 parcellation).
  325. # The CSV files are in parcellation order: first 16 rows/cols = subcortical,
  326. # rows 16–116 = Schaefer cortical parcels. We reorder so that cortical comes
  327. # first (rows 0–99) and subcortical last (rows 100–115) to match the HOI data.
  329. list_SC_edges = {}
  330. for sc_path in sorted(glob.glob(SC_PATH + "SC_100_*nof*csv")):
  331. subj_id = sc_path.split("/")[-1].split("_")[3]
  332. raw = np.loadtxt(sc_path, delimiter=",")
  333. SC = np.zeros((N_ROIS, N_ROIS))
  334. SC[:100, :100] = raw[16:116, 16:116] # cortical–cortical
  335. SC[100:, 100:] = raw[:16, :16] # subcortical–subcortical
  336. SC[100:, :100] = raw[:16, 16:] # subcortical–cortical
  337. SC += SC.T # symmetrise
  338. list_SC_edges[subj_id] = upper_tri_masking(SC)
  340. SC_arr = np.array(list(list_SC_edges.values()))
  341. avg_SC = np.mean(SC_arr, axis=0)
  342. print(f"SC loaded — {SC_arr.shape[0]} subjects, {SC_arr.shape[1]} edges each")
  345. # %% [markdown]
  346. # ## 4. Panel (a) — Brain maps and Yeo network distributions
  347. #
  348. # Each HOI measure is projected to node strength, ranked across regions, and
  349. # visualised on an inflated cortical surface. Distributions across the seven
  350. # Yeo resting-state networks (+ subcortical) are also shown.
  352. # %%
  353. # ── Yeo parcellation ─────────────────────────────────────────────────────────
  354. def load_yeo(filename="../utils/yeo_RS7_Schaefer100S.mat", n_rois=N_ROIS):
  355. """Return yeoOrder, yeoROIs array, and per-network ROI index dict."""
  356. mat = sio.loadmat(filename)
  357. yeo_order = np.array([i - 1 for i in mat["yeoOrder"]])[:n_rois]
  358. yeo_rois = np.array([i[0] - 1 for i in mat["yeoROIs"]])[:n_rois]
  359. yeo_names = ["VIS", "SM", "DA", "VA", "L", "FP", "DMN", "SC"]
  360. yeo_dict = {name: np.where(yeo_rois == k)[0]
  361. for k, name in enumerate(yeo_names)}
  362. # Subcortical: merge Yeo label 7 with any label 8 entries
  363. yeo_dict["SC"] = np.concatenate([yeo_dict["SC"], np.where(yeo_rois == 8)[0]])
  364. return yeo_order, yeo_rois, yeo_dict
  366. yeo_order, yeo_rois, yeo_dict = load_yeo()
  367. # Force subcortical to the 16 appended regions
  368. yeo_dict["SC"] = np.arange(100, N_ROIS)
  370. # ── Helper functions ──────────────────────────────────────────────────────────
  371. def compute_ranks(v):
  372. """Return floor ranks of absolute values."""
  373. return np.floor(rankdata(np.abs(v), method="average"))
  375. def aggregate_yeo_ranks(data, yeo_dict):
  376. """Per-subject rank each nodal vector, then average within each network."""
  377. result = {net: [] for net in yeo_dict}
  378. for subj in data:
  379. r = compute_ranks(subj)
  380. for net in yeo_dict:
  381. result[net].append(np.mean(r[yeo_dict[net]]))
  382. return result
  384. SIZE_LARGE, SIZE_SMALL = 26, 19
  386. def plot_brain_map(data, filename, force_white=False, graymap_rev=True, n_rois=100):
  387. """Plot mean ranked node strength on a cortical surface."""
  388. mean_ranks = np.mean([compute_ranks(v) for v in data], axis=0)[:n_rois]
  389. vmin = np.nanpercentile(mean_ranks, 10)
  390. vmax = np.nanpercentile(mean_ranks, 90)
  391. fig = normal_view_top(
  392. mean_ranks, edges=True, brightness=0.85,
  393. parcellation=100, center_cbar=True, graymap_rev=graymap_rev,
  394. cmap=CMAP_STANDARD, alpha_graymap=1, parcellation_name="schaefer",
  395. exp_form=False, vmin=vmin, vmax=vmax, force_white=force_white,
  396. surftype="inflated",
  397. )
  398. fig.axes[1].set_xlabel(r"$\langle s_i \rangle$" + "\n\npercentile",
  399. labelpad=-70, fontsize=SIZE_LARGE)
  400. fig.axes[1].set_xticks([vmin, vmax])
  401. fig.axes[1].set_xticklabels(["10", "90"], fontsize=SIZE_LARGE)
  402. fig.savefig(f"{FIG_PATH}{filename}.png", dpi=300, bbox_inches="tight", transparent=True)
  403. plt.close(fig)
  405. def plot_yeo_distribution(data, filename):
  406. """Strip + box plot of ranked node strength per Yeo network."""
  407. df = pd.DataFrame(aggregate_yeo_ranks(data, yeo_dict))
  408. fig, ax = plt.subplots(1, 1, figsize=(6, 2.5), dpi=150)
  409. sns.stripplot(data=df, alpha=1, edgecolor="w", linewidth=1, zorder=-10,
  410. ax=ax, palette="muted", size=5.5)
  411. sns.boxplot(data=df, fill=False, width=0.5, ax=ax, showfliers=False)
  412. ax.spines[["top", "right"]].set_visible(False)
  413. ax.set_xticklabels(ax.get_xticklabels(), rotation=0, ha="center", fontsize=SIZE_SMALL)
  414. ax.set_yticks(np.arange(0, 122, 40))
  415. ax.set_yticklabels(np.arange(0, 122, 40), fontsize=SIZE_SMALL)
  416. plt.tight_layout()
  417. plt.savefig(f"{FIG_PATH}{filename}_yeo.png", dpi=300, bbox_inches="tight", transparent=True)
  418. plt.close(fig)
  421. # %%
  422. # Generate brain maps and Yeo distributions for all measures.
  423. # Brain maps use inflated surface; Synergy is plotted with absolute values
  424. # because it is negative-valued (synergy-dominated triplets have O-info < 0).
  426. brain_map_specs = [
  427. # (data, filename, force_white, graymap_rev)
  428. (Red_nodal, "Redundancy", False, True),
  429. (np.abs(Syn_nodal), "Synergy", True, False),
  430. (triangles_nodal, "Triangles", True, False),
  431. (scaffold_nodal, "HO_Scaffold", False, True),
  432. (scaffold_freq_nodal, "Scaffold_freq", True, True),
  433. (scaffold_pers_nodal, "Scaffold_pers", True, True),
  434. (PED_Red_nodal, "PED_Red", True, False),
  435. (PED_Syn_nodal, "PED_Syn", True, False),
  436. (phiIDRed_nodal, "PhiID_Red", True, False),
  437. (phiIDSyn_nodal, "PhiID_Syn", True, True),
  438. (classical_FC_nodal, "Classical_FC", True, True),
  439. ]
  441. for data, fname, fw, gr in brain_map_specs:
  442. plot_brain_map(data, fname, force_white=fw, graymap_rev=gr)
  443. plot_yeo_distribution(data, fname)
  444. print(f" ✓ {fname}")
  446. print("All brain maps saved.")
  449. # %% [markdown]
  450. # ## 5. Panel (a) — Edge-level connectivity matrices
  451. #
  452. # Lower-triangular visualisation of the mean edge-level connectivity matrix for
  453. # each HOI measure. Rows and columns correspond to brain regions.
  455. # %%
  456. # --- Note: Matplotlib rcParams for font may not always have effect, especially for non-PDF outputs or if font is not found ---
  457. def plot_edge_matrix(edge_data, save_path):
  458. """Plot the group-average edge matrix as a lower-triangular image."""
  459. mean_vec = np.mean(edge_data, axis=0)
  460. sq = squareform(mean_vec)
  461. # Mask upper triangle (including diagonal)
  462. mask = np.triu(np.ones_like(sq, dtype=bool))
  463. sq[mask] = np.nan
  465. vmin = np.percentile(mean_vec, 10)
  466. vmax = np.percentile(mean_vec, 90)
  468. # NOTE: Explicitly set font in this plot function for axes
  469. fig, ax = plt.subplots(figsize=(8, 8))
  470. ax.imshow(sq, vmin=vmin, vmax=vmax, cmap=CMAP_STANDARD)
  472. n = sq.shape[0]
  473. ax.plot([-0.25, n - 0.75], [-0.25, n - 0.75], "k-", linewidth=2)
  475. # Thin white grid lines clipped precisely at the diagonal (data coordinates)
  476. for i in range(1, n):
  477. ax.plot([-0.5, i - 0.5], [i - 0.5, i - 0.5], "w-", lw=0.25) # horizontal → diagonal
  478. ax.plot([i - 0.5, i - 0.5], [i - 0.5, n - 0.5], "w-", lw=0.25) # vertical ↓ bottom
  480. ax.spines[["top", "right"]].set_visible(False)
  481. ax.spines[["bottom", "left"]].set_linewidth(2)
  482. ax.set_xticks([])
  483. ax.set_yticks([])
  485. # Try to force Helvetica/Arial via set_fontname on tick labels (if rcParams isn't working)
  486. for label in ax.get_xticklabels() + ax.get_yticklabels():
  487. label.set_fontname('Helvetica') # will fallback if not available
  489. plt.tight_layout()
  490. plt.savefig(save_path, dpi=300, transparent=True, bbox_inches="tight")
  491. plt.close(fig)
  494. # Generate all edge matrices
  495. edge_matrix_specs = [
  496. (Red_edges, "Red_edges"),
  497. (Syn_edges, "Syn_edges"),
  498. (triangles_edges, "Triangles_edges"),
  499. (scaffold_edges, "Scaffold_edges"),
  500. (scaffold_freq_edges, "Scaffold_freq_edges"),
  501. (scaffold_pers_edges, "Scaffold_pers_edges"),
  502. (PED_Red_edges, "PED_Red_edges"),
  503. (PED_Syn_edges, "PED_Syn_edges"),
  504. (phiIDRed_edges, "PhiID_Red_edges"),
  505. (phiIDSyn_edges, "PhiID_Syn_edges"),
  506. (classical_FC_edges, "Classical_FC_edges"),
  507. ]
  509. for data, fname in edge_matrix_specs:
  510. plot_edge_matrix(data, f"{FIG_PATH}{fname}.png")
  511. print(f" ✓ {fname}")
  513. # Standalone colourbar
  514. fig, ax = plt.subplots(figsize=(6, 0.3))
  515. fig.subplots_adjust(bottom=0.5)
  516. norm = mcolors.Normalize(vmin=0, vmax=1)
  517. cbar = plt.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=CMAP_STANDARD),
  518. cax=ax, orientation="horizontal")
  519. cbar.ax.set_xticks([])
  520. cbar.ax.set_yticks([])
  522. # Explicitly set font for colorbar, too (in case global rcParams isn't working)
  523. for label in cbar.ax.get_xticklabels() + cbar.ax.get_yticklabels():
  524. label.set_fontname('Helvetica')
  526. plt.savefig(f"{FIG_PATH}colorbar.png", dpi=300, bbox_inches="tight", transparent=True)
  527. plt.close(fig)
  528. print("Edge matrices and colourbar saved.")
  531. # %% [markdown]
  532. # ## 6. Panels (b) & (c) — Edge-level similarity matrix and hierarchical clustering
  533. #
  534. # Panel (b): Spearman correlation between the group-average edge profiles of each
  535. # measure pair. Panel (c): Complete-linkage dendrogram of those distances.
  537. # %%
  538. # Compute mean edge vector per measure (absolute value for Syn to reflect magnitude)
  539. nodal_all = np.array([
  540. np.mean(np.abs(data_edges[key]), axis=0) if key == "Syn"
  541. else np.mean(data_edges[key], axis=0)
  542. for key, _ in MEASURES
  543. ])
  545. # Spearman correlation matrix (symmetric by construction)
  546. n_meas = len(MEASURES)
  547. corr_matrix = np.zeros((n_meas, n_meas))
  548. for i in range(n_meas):
  549. for j in range(n_meas):
  550. corr_matrix[i, j], _ = spearmanr(nodal_all[i], nodal_all[j])
  551. corr_matrix = (corr_matrix + corr_matrix.T) / 2 # enforce symmetry numerically
  553. print("Spearman correlation matrix computed.")
  556. # %%
  557. # ── Panel (b) + (c): lower-triangle correlation matrix + complete-linkage dendogram ──
  559. # Cluster-group colours for labels and branches
  560. SYNERGY_SET = {"O-Syn", "PED Syn", "PhiID Syn"}
  561. REDUNDANCY_SET= {"FC", "Triangles", "O-Red", "PED Red", "PhiID Red"}
  562. TOPOLOGY_SET = {"Fscaff", "Pscaff", "TScaffold"}
  564. def _group_color(label):
  565. if label in SYNERGY_SET: return "blue"
  566. if label in REDUNDANCY_SET: return "red"
  567. if label in TOPOLOGY_SET: return "purple"
  568. return "black"
  570. # Complete-linkage clustering
  571. dist_sq = np.sqrt(2 * np.clip(1 - corr_matrix, 0, None)) # correlation → distance
  572. Z_complete = linkage(squareform(dist_sq, checks=False), method="complete")
  574. def _leaf_group_color(node_idx, Z, n_leaves):
  575. """Color a dendrogram branch by the group of its descendant leaves."""
  576. def _leaves(idx):
  577. if idx < n_leaves:
  578. return {idx}
  579. l, r = int(Z[idx - n_leaves][0]), int(Z[idx - n_leaves][1])
  580. return _leaves(l) | _leaves(r)
  582. leaf_set = _leaves(node_idx)
  583. groups = {_group_color(LABELS[i]) for i in leaf_set}
  584. groups.discard("black")
  585. return groups.pop() if len(groups) == 1 else "gray"
  587. fig = plt.figure(dpi=150, figsize=(16, 7))
  588. gs = plt.GridSpec(1, 8)
  590. # ─ Left: lower-triangle Spearman matrix ─
  591. ax1 = fig.add_subplot(gs[0, :5])
  592. mask_upper = np.triu(np.ones_like(corr_matrix, dtype=bool), k=0)
  594. import seaborn as sns
  595. sns.heatmap(
  596. np.round(corr_matrix, 2),
  597. annot=True, fmt=".2f", cmap="RdYlBu_r",
  598. mask=mask_upper, ax=ax1, cbar=False,
  599. linewidths=0., square=True,
  600. annot_kws={
  601. "size": 13,
  602. "fontname": "Helvetica"
  603. }
  604. )
  606. ax1.set_yticks(np.arange(n_meas) + 0.5)
  607. ax1.set_yticklabels(LABELS, rotation=0, fontname='Helvetica')
  608. ax1.set_xticks(np.arange(n_meas - 1) + 0.5)
  609. ax1.set_xticklabels(LABELS[:-1], rotation=45, fontname='Helvetica')
  610. plt.setp(ax1.get_yticklabels(), rotation=45, ha="right", rotation_mode="anchor", fontsize=14, fontname="Helvetica")
  611. plt.setp(ax1.get_xticklabels(), rotation=45, ha="right", rotation_mode="anchor", fontsize=14, fontname="Helvetica")
  612. ax1.tick_params(length=7, width=1.5)
  614. bw = 1.5
  615. for i in range(1, n_meas):
  616. ax1.plot([i-1, i], [i, i], "k-", lw=bw, zorder=100, clip_on=False)
  617. ax1.plot([i, i], [i, i+1], "k-", lw=bw, zorder=100, clip_on=False)
  618. ax1.plot([0, 0], [1, n_meas], "k-", lw=bw+0.5, zorder=100, clip_on=False)
  619. ax1.plot([0, n_meas-1], [n_meas, n_meas], "k-", lw=bw+0.5, zorder=100, clip_on=False)
  620. for x in range(1, n_meas):
  621. ax1.axvline(x, ymin=0, ymax=1 - 0.1*x, color="white", lw=bw)
  622. ax1.axhline(x, xmin=0, xmax=0.09*x, color="white", lw=bw)
  623. ax1.set_ylim(n_meas, 1)
  625. # ─ Right: complete-linkage dendrogram ─
  626. ax2 = fig.add_subplot(gs[0, 5:])
  627. n_leaves = len(LABELS)
  628. dendro = dendrogram(
  629. Z_complete, labels=LABELS, orientation="top",
  630. color_threshold=0,
  631. link_color_func=lambda k: _leaf_group_color(k, Z_complete, n_leaves),
  632. ax=ax2, leaf_rotation=45,
  633. )
  634. ax2.spines[["top", "left", "right"]].set_visible(False)
  635. ax2.spines["bottom"].set_linewidth(1.5)
  636. ax2.set_yticks([])
  637. for tick in ax2.get_xticklabels():
  638. tick.set_color(_group_color(tick.get_text()))
  639. # Set font to Helvetica, fallback to Arial if not available
  640. try:
  641. tick.set_fontname("Helvetica")
  642. except Exception:
  643. tick.set_fontname("Arial")
  644. ax2.tick_params(length=10, width=1)
  645. plt.xticks(rotation=45, ha="right", fontsize=13, fontname="Helvetica")
  647. plt.subplots_adjust(wspace=0)
  648. plt.savefig(f"{FIG_PATH}Fig2_bc_similarity_dendrogram_complete.pdf",
  649. dpi=300, bbox_inches="tight", transparent=True)
  650. plt.savefig(f"{FIG_PATH}Fig2_bc_similarity_dendrogram_complete.svg")
  651. print("Panel (b)+(c) [complete linkage] saved.")
  654. # %%
  655. # ── Alternative: Ward-linkage clustering with ordered heatmap ─────────────────
  656. # This version reorders the correlation matrix rows/columns according to the
  657. # Ward dendrogram, making cluster structure visible directly in the heatmap.
  659. Z_ward = linkage(squareform(dist_sq, checks=False), method="ward", optimal_ordering=True)
  660. idx = hc_leaves_list(Z_ward)
  662. # Reorder matrix and labels by dendrogram leaf order
  663. corr_ordered = corr_matrix[np.ix_(idx[::-1], idx[::-1])]
  664. labels_ordered= [LABELS[i] for i in idx[::-1]]
  666. fig = plt.figure(figsize=(16, 7), dpi=150)
  667. gs = GridSpec(1, 20, figure=fig)
  669. ax_hm = fig.add_subplot(gs[0, :19])
  670. sns.heatmap(corr_ordered, annot=True, fmt=".2f", cmap="RdYlBu_r",
  671. ax=ax_hm, cbar=False, linewidths=0.5, square=True,
  672. annot_kws={"size": 12})
  673. ax_hm.set_xticks(np.arange(n_meas) + 0.5)
  674. ax_hm.set_yticks(np.arange(n_meas) + 0.5)
  675. ax_hm.set_xticklabels(labels_ordered, rotation=45, ha="right", fontsize=13)
  676. ax_hm.set_yticklabels(labels_ordered, rotation=0, fontsize=13)
  677. for tick in ax_hm.get_xticklabels():
  678. tick.set_color(_group_color(tick.get_text()))
  679. for tick in ax_hm.get_yticklabels():
  680. tick.set_color(_group_color(tick.get_text()))
  682. # Ward dendrogram on the right
  683. _ward_colors = (["darkred"]*11 + ["red"]*2 + ["purple"] + ["red"]*2 +
  684. ["blue"]*2 + ["purple"] + ["gray"]*5)
  685. ax_dg = fig.add_subplot(gs[0, 19:])
  686. dendrogram(Z_ward, orientation="right", labels=LABELS, leaf_font_size=13,
  687. ax=ax_dg, color_threshold=2.5,
  688. link_color_func=lambda k: _ward_colors[k])
  689. ax_dg.spines[["top", "right", "bottom"]].set_visible(False)
  690. ax_dg.set_yticks([])
  691. ax_dg.set_xticks([])
  693. plt.subplots_adjust(wspace=-0.8775)
  694. plt.savefig(f"{FIG_PATH}Fig2_bc_similarity_dendrogram_ward.pdf",
  695. dpi=300, bbox_inches="tight", transparent=True)
  696. print("Panel (b)+(c) [Ward linkage] saved.")
  699. # %% [markdown]
  700. # ## 7. Panel (d) — Structure–function cartography
  701. #
  702. # Each HOI measure is characterised by two Spearman correlations computed from
  703. # its group-average edge-level profile:
  704. # - **x-axis**: correlation with mean structural connectivity (SC)
  705. # - **y-axis**: correlation with mean functional connectivity (FC)
  707. # %%
  708. # Average edge-level profiles across subjects (used for scatter)
  709. ROIs_used = N_ROIS
  711. def _avg_edge(arr):
  712. """Mean over subjects of the upper-triangle restricted to ROIs_used × ROIs_used."""
  713. return np.mean([upper_tri_masking(squareform(arr[i])[:ROIs_used, :ROIs_used])
  714. for i in range(len(arr))], axis=0)
  716. avg_FC = _avg_edge(classical_FC_edges)
  717. avg_Red = _avg_edge(Red_edges)
  718. avg_Syn = _avg_edge(np.abs(Syn_edges)) # absolute value for magnitude
  719. avg_tri = _avg_edge(triangles_edges)
  720. avg_scaf = _avg_edge(scaffold_edges)
  721. avg_freq = _avg_edge(scaffold_freq_edges)
  722. avg_pers = _avg_edge(scaffold_pers_edges)
  723. avg_PEDR = _avg_edge(PED_Red_edges)
  724. avg_PEDS = _avg_edge(PED_Syn_edges)
  725. avg_phiS = _avg_edge(phiIDSyn_edges)
  726. avg_phiR = _avg_edge(phiIDRed_edges)
  728. SF = {
  729. "O-Red": (spearmanr(avg_SC, avg_Red).statistic, spearmanr(avg_FC, avg_Red).statistic),
  730. "O-Syn": (spearmanr(avg_SC, avg_Syn).statistic, spearmanr(avg_FC, avg_Syn).statistic),
  731. "Triangles":(spearmanr(avg_SC, avg_tri).statistic, spearmanr(avg_FC, avg_tri).statistic),
  732. "TScaffold":(spearmanr(avg_SC, avg_scaf).statistic, spearmanr(avg_FC, avg_scaf).statistic),
  733. "Fscaff": (spearmanr(avg_SC, avg_freq).statistic, spearmanr(avg_FC, avg_freq).statistic),
  734. "Pscaff": (spearmanr(avg_SC, avg_pers).statistic, spearmanr(avg_FC, avg_pers).statistic),
  735. "PED Red": (spearmanr(avg_SC, avg_PEDR).statistic, spearmanr(avg_FC, avg_PEDR).statistic),
  736. "PED Syn": (spearmanr(avg_SC, avg_PEDS).statistic, spearmanr(avg_FC, avg_PEDS).statistic),
  737. "PhiID Syn":(spearmanr(avg_SC, avg_phiS).statistic, spearmanr(avg_FC, avg_phiS).statistic),
  738. "PhiID Red":(spearmanr(avg_SC, avg_phiR).statistic, spearmanr(avg_FC, avg_phiR).statistic),
  739. "FC": (spearmanr(avg_SC, avg_FC).statistic, spearmanr(avg_FC, avg_FC).statistic),
  740. }
  742. for name, (r_SC, r_FC) in SF.items():
  743. print(f" {name:12s} SC: {r_SC:+.3f} FC: {r_FC:+.3f}")
  746. # %%
  747. fig, ax = plt.subplots(figsize=(6.75, 6), dpi=150)
  748. for spine in ax.spines.values():
  749. spine.set_linewidth(1.5)
  751. for method, (r_SC, r_FC) in SF.items():
  752. mk, col, sz = MARKER_STYLE[method]
  753. ax.scatter(
  754. r_SC, r_FC, label=method, marker=mk, s=sz, c=col,
  755. edgecolor="black", linewidth=1
  756. )
  758. ax.axvline(0, color="black", lw=1, ls=":", alpha=0.1)
  759. ax.axhline(0, color="black", lw=1, ls=":", alpha=0.1)
  760. ax.set_xlabel(
  761. "Structural connectivity", fontsize=20, labelpad=5,
  762. fontname="Helvetica"
  763. )
  764. ax.set_ylabel(
  765. "Functional connectivity", fontsize=20, labelpad=0,
  766. fontname="Helvetica"
  767. )
  768. ax.tick_params(width=1.5, length=6, labelsize=16)
  769. for label in (ax.get_xticklabels() + ax.get_yticklabels()):
  770. label.set_fontname("Helvetica")
  771. ax.set_xlim(-1.15, 1.15)
  772. ax.set_ylim(-1.15, 1.15)
  773. leg = ax.legend(
  774. loc="upper left", bbox_to_anchor=(-0.005, 1.005), ncol=2,
  775. columnspacing=1., handletextpad=0.0, labelspacing=1,
  776. fontsize=11.5, markerscale=1.1, frameon=True, fancybox=True, shadow=True
  777. )
  778. for text in leg.get_texts():
  779. text.set_fontname("Helvetica")
  780. plt.tight_layout()
  781. plt.savefig(
  782. f"{FIG_PATH}Fig2d_structure_function_edge.pdf",
  783. dpi=300, bbox_inches="tight", transparent=True
  784. )
  785. print("Panel (d) saved.")
  788. # %% [markdown]
  789. # ## 8. Supplementary — Nodal-level analyses
  790. #
  791. # Repeats the similarity matrix, dendrogram, and structure–function scatter at
  792. # the nodal (node-strength) projection level instead of edge level.
  794. # %%
  795. # Nodal SC: per-subject degree-strength (mean SC per node), then group average
  796. list_SC_nodal = {}
  797. for sc_path in sorted(glob.glob(SC_PATH + "santo_100_*nof*csv")):
  798. subj_id = sc_path.split("/")[-1].split("_")[3]
  799. raw = np.loadtxt(sc_path, delimiter=",")
  800. SC = np.zeros((N_ROIS, N_ROIS))
  801. SC[:100, :100] = raw[16:116, 16:116]
  802. SC[100:, 100:] = raw[:16, :16]
  803. SC[100:, :100] = raw[:16, 16:]
  804. SC += SC.T
  805. list_SC_nodal[subj_id] = np.mean(SC, axis=0)
  807. SC_nodal_arr = np.array(list(list_SC_nodal.values()))
  808. avg_SC_nodal = np.mean(SC_nodal_arr, axis=0)
  809. print(f"Nodal SC loaded — {SC_nodal_arr.shape[0]} subjects, {SC_nodal_arr.shape[1]} nodes each")
  812. # %%
  813. # Nodal structure-function correlations
  814. def _avg_nodal(arr):
  815. """Group-average nodal profile, restricted to ROIs_used."""
  816. return np.mean([np.mean(squareform(arr[i])[:ROIs_used, :ROIs_used], axis=0)
  817. for i in range(len(arr))], axis=0)
  819. avg_FC_n = _avg_nodal(classical_FC_edges)
  820. avg_Red_n = _avg_nodal(Red_edges)
  821. avg_Syn_n = _avg_nodal(np.abs(Syn_edges))
  822. avg_tri_n = _avg_nodal(triangles_edges)
  823. avg_scaf_n = _avg_nodal(scaffold_edges)
  824. avg_freq_n = _avg_nodal(scaffold_freq_edges)
  825. avg_pers_n = _avg_nodal(scaffold_pers_edges)
  826. avg_PEDR_n = _avg_nodal(PED_Red_edges)
  827. avg_PEDS_n = _avg_nodal(PED_Syn_edges)
  828. avg_phiS_n = _avg_nodal(phiIDSyn_edges)
  829. avg_phiR_n = _avg_nodal(phiIDRed_edges)
  831. SF_nodal = {
  832. "O-Red": (spearmanr(avg_SC_nodal, avg_Red_n).statistic, spearmanr(avg_FC_n, avg_Red_n).statistic),
  833. "O-Syn": (spearmanr(avg_SC_nodal, avg_Syn_n).statistic, spearmanr(avg_FC_n, avg_Syn_n).statistic),
  834. "Triangles":(spearmanr(avg_SC_nodal, avg_tri_n).statistic, spearmanr(avg_FC_n, avg_tri_n).statistic),
  835. "TScaffold":(spearmanr(avg_SC_nodal, avg_scaf_n).statistic, spearmanr(avg_FC_n, avg_scaf_n).statistic),
  836. "Fscaff": (spearmanr(avg_SC_nodal, avg_freq_n).statistic, spearmanr(avg_FC_n, avg_freq_n).statistic),
  837. "Pscaff": (spearmanr(avg_SC_nodal, avg_pers_n).statistic, spearmanr(avg_FC_n, avg_pers_n).statistic),
  838. "PED Red": (spearmanr(avg_SC_nodal, avg_PEDR_n).statistic, spearmanr(avg_FC_n, avg_PEDR_n).statistic),
  839. "PED Syn": (spearmanr(avg_SC_nodal, avg_PEDS_n).statistic, spearmanr(avg_FC_n, avg_PEDS_n).statistic),
  840. "PhiID Syn":(spearmanr(avg_SC_nodal, avg_phiS_n).statistic, spearmanr(avg_FC_n, avg_phiS_n).statistic),
  841. "PhiID Red":(spearmanr(avg_SC_nodal, avg_phiR_n).statistic, spearmanr(avg_FC_n, avg_phiR_n).statistic),
  842. "FC": (spearmanr(avg_SC_nodal, avg_FC_n).statistic, spearmanr(avg_FC_n, avg_FC_n).statistic),
  843. }
  845. fig, ax = plt.subplots(figsize=(6.75, 6), dpi=150)
  846. for spine in ax.spines.values():
  847. spine.set_linewidth(1.5)
  848. for method, (r_SC, r_FC) in SF_nodal.items():
  849. mk, col, sz = MARKER_STYLE[method]
  850. ax.scatter(r_SC, r_FC, label=method, marker=mk, s=sz,
  851. c=col, edgecolor="black", linewidth=1)
  852. ax.axvline(0, color="black", lw=1, ls=":", alpha=0.1)
  853. ax.axhline(0, color="black", lw=1, ls=":", alpha=0.1)
  854. ax.set_xlabel("Structural connectivity", fontsize=20, labelpad=10)
  855. ax.set_ylabel("Functional connectivity", fontsize=20, labelpad=5)
  856. ax.tick_params(width=1.5, length=6, labelsize=16)
  857. ax.set_xlim(-1.15, 1.15)
  858. ax.set_ylim(-1.15, 1.15)
  859. ax.legend(loc="upper left", bbox_to_anchor=(-0.005, 1.005), ncol=2,
  860. columnspacing=1., handletextpad=0.0, labelspacing=1,
  861. fontsize=12.5, markerscale=1.1, frameon=True, fancybox=True, shadow=True)
  862. plt.tight_layout()
  863. plt.savefig(f"{FIG_PATH}Fig2_structure_function_nodal.pdf",
  864. dpi=300, bbox_inches="tight", transparent=True)
  865. print("Nodal structure-function scatter saved.")
  868. # %%
  869. # Nodal similarity matrix + complete-linkage dendrogram
  870. nodal_nodal_all = np.array([
  871. np.mean(np.abs(data_nodal[key]), axis=0) if key == "Syn"
  872. else np.mean(data_nodal[key], axis=0)
  873. for key, _ in MEASURES
  874. ])
  876. corr_nodal = np.zeros((n_meas, n_meas))
  877. for i in range(n_meas):
  878. for j in range(n_meas):
  879. corr_nodal[i, j], _ = spearmanr(nodal_nodal_all[i], nodal_nodal_all[j])
  880. corr_nodal = (corr_nodal + corr_nodal.T) / 2
  882. dist_nodal = np.sqrt(2 * np.clip(1 - corr_nodal, 0, None))
  883. Z_nod_comp = linkage(squareform(dist_nodal, checks=False), method="complete")
  885. fig = plt.figure(dpi=150, figsize=(16, 7))
  886. gs = plt.GridSpec(1, 8)
  888. ax1 = fig.add_subplot(gs[0, :5])
  889. mask_upper = np.triu(np.ones_like(corr_nodal, dtype=bool), k=0)
  890. sns.heatmap(np.round(corr_nodal, 2),
  891. annot=True, fmt=".2f", cmap="RdYlBu_r", mask=mask_upper,
  892. ax=ax1, cbar=False, linewidths=0., square=True, annot_kws={"size": 12})
  893. ax1.set_yticks(np.arange(n_meas) + 0.5)
  894. ax1.set_yticklabels(LABELS, rotation=0)
  895. ax1.set_xticks(np.arange(n_meas - 1) + 0.5)
  896. ax1.set_xticklabels(LABELS[:-1], rotation=45)
  897. plt.setp(ax1.get_yticklabels(), rotation=45, ha="right", rotation_mode="anchor", fontsize=13)
  898. plt.setp(ax1.get_xticklabels(), rotation=45, ha="right", rotation_mode="anchor", fontsize=13)
  899. ax1.tick_params(length=7, width=1.5)
  900. for i in range(1, n_meas):
  901. ax1.plot([i-1, i], [i, i], "k-", lw=1.5, zorder=100, clip_on=False)
  902. ax1.plot([i, i], [i, i+1], "k-", lw=1.5, zorder=100, clip_on=False)
  903. ax1.plot([0, 0], [1, n_meas], "k-", lw=2, zorder=100, clip_on=False)
  904. ax1.plot([0, n_meas-1], [n_meas, n_meas], "k-", lw=2, zorder=100, clip_on=False)
  905. for x in range(1, n_meas):
  906. ax1.axvline(x, ymin=0, ymax=1-0.1*x, color="white", lw=1.5)
  907. ax1.axhline(x, xmin=0, xmax=0.09*x, color="white", lw=1.5)
  908. ax1.set_ylim(n_meas, 1)
  910. ax2 = fig.add_subplot(gs[0, 5:])
  911. dendro = dendrogram(
  912. Z_nod_comp, labels=LABELS, orientation="top",
  913. color_threshold=0,
  914. #link_color_func=lambda k: _leaf_group_color(k, Z_nod_comp, n_leaves),
  915. ax=ax2, leaf_rotation=45,
  916. )
  917. ax2.spines[["top", "left", "right"]].set_visible(False)
  918. ax2.spines["bottom"].set_linewidth(1.5)
  919. ax2.set_yticks([])
  920. for tick in ax2.get_xticklabels():
  921. tick.set_color(_group_color(tick.get_text()))
  922. ax2.tick_params(length=10, width=1)
  923. plt.xticks(rotation=45, ha="right", fontsize=13)
  925. plt.subplots_adjust(wspace=0)
  926. plt.savefig(f"{FIG_PATH}Fig2_nodal_similarity_dendrogram.pdf",
  927. dpi=300, bbox_inches="tight", transparent=True)
  928. print("Nodal similarity + dendrogram saved.")
  931. # %%

Fig2_analyses_clean.ipynb at commit ce1a32c, no license · at the source

Overview

Authors: Andrea Santoro1, Matteo Neri2, Simone Poetto3,4, Davide Orsenigo5, Matteo Diano5, Marilyn Gatica6,7, Giovanni Petri7,8
ORCID iDs: Andrea Santoro, Matteo Neri, Davide Orsenigo, Marilyn Gatica, Giovanni Petri
  1. ISI Foundation,Torino, Italy
  2. Institut de Neurosciences de la Timone UMR 7289, Aix Marseille Université, CNRS,Marseille, France
  3. Nicolaus Copernicus University,Toruń, Poland
  4. Intesa Sanpaolo A.I. Research, Torino, Italy
  5. Department of Psychology and Neuroscience Institute of Turin, University of Torino,Torino, Italy
  6. Centre for Apprenticeships, Northeastern University London,London, UK
  7. NPLab, Network Science Institute, Northeastern University London,London, UK
  8. Department of Physics, Northeastern University,Boston, MA, USA
Institutions: Institute for Scientific Interchange (Italy); Institut de Neurosciences de la Timone (France); Nicolaus Copernicus University (Poland); Neuroscience Institute (Italy); Northeastern University London (United Kingdom); Northeastern University (United States)
Journal: Nature communications, volume 17, issue 1, article 9207
Dates: received 9 August 2025; accepted 16 July 2026; published online 29 July 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1038/s41467-026-75959-w · PMID 42660887 · PMCID PMC13522477 · OpenAlex W7171682201
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: human (organism)
Methods: Spectral & time-frequency, Connectivity, Statistics, Smoothing, state filtering, decompositions, Machine learning, Preprocessing, Graphs, fMRI & imaging
Keywords: Network models, Applied mathematics
MeSH: Brain*, Connectome*, Models, Neurological*, Adult, Brain Mapping, Female, Humans, Magnetic Resonance Imaging, Male, Nerve Net (* major topic)
Topic: Functional Brain Connectivity Studies (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: EC | EC Seventh Framework Programm | FP7 People: Marie-Curie Actions (FP7-PEOPLE - Specific Programme "People" Implementing the Seventh Framework Programme of the European Community for Research, Technological Development and Demonstration Activities (2007 to 2013)) (101120085); EC | EC Seventh Framework Programm | FP7 Ideas: European Research Council (FP7-IDEAS-ERC - Specific Programme: "Ideas" Implementing the Seventh Framework Programme of the European Community for Research, Technological Development and Demonstration Activities (2007 to 2013)) (101171380)
Citations: cited by 1 paper (Europe PMC); 181 references in the paper

Abstract

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

Repositories

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

andresantoro/RHOSTS

License: GPL-3.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: f42f531c4434ddc5a7b1144ceadc7b246a31a3f8, 3 October 2024
Languages: Python (8), Shell (7), Julia (3), Jupyter (1)
Size: 42 files, 19 scripts
Software Heritage: not archived
Found in: the text, “Triangles and temporal scaffold”
Holds: README, license file, environment (requirements.txt, Containers/Dockerfile, Containers/RHOSTS.def, High_order_TS/requirements.txt, High_order_TS_with_scaffold/requirements.txt, Julia_MTS_optimized/SimplicialTS/Manifest.toml, Julia_MTS_optimized/SimplicialTS/Project.toml), 1 notebook
Not found: CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (5 files), h5py (3 files), NetworkX (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
21 files
At the source: github.com/andresantoro/RHOSTS

Zenodo 20260959

License: CC-BY-4.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Size: 16 files
Software Heritage: not checked
Found in: “Data availability”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
1 file
At the source:

nplresearch/HOI_lenses_analysis

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: ce1a32ccba63dcc20c495b77be7352ad2993575d, 18 May 2026
Languages: Python (32), Shell (19), MATLAB (12), Jupyter (6), Julia (4)
Size: 128 files, 73 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, environment (Container/requirements_julia.txt, Container/requirements_python.txt, Code/02_local_higher_order_TDA/SimplicialTS/Manifest.toml, Code/02_local_higher_order_TDA/SimplicialTS/Project.toml), 6 notebooks
Not found: license file, CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (34 files), Matplotlib (23 files), SciPy (20 files), h5py (10 files), BrainSpace (8 files), pandas (8 files), seaborn (8 files), NiBabel (5 files), neuromaps (4 files), Nilearn (4 files), scikit-learn (4 files), Statistics and Machine Learning Toolbox (3 files), NetworkX (3 files), statsmodels (3 files), fdr_bh (Benjamini-Hochberg FDR) (1 file), Numba (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
74 files
At the source: github.com/nplresearch/HOI_lenses_analysis

brainets/hoi

License: BSD-3-Clause
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 4db2fbd701d40d33fe2375f3df5a4a4dc94c6b83, 30 January 2026
Languages: Python (58)
Size: 109 files, 58 scripts
Software Heritage: archived
Found in: “Code availability”
Holds: README, license file, CITATION.cff, environment (pyproject.toml, requirements.txt, setup.cfg, setup.py), tests, continuous integration, documentation
Tools: NumPy (44 files), JAX (25 files), Matplotlib (15 files), scikit-learn (4 files), pandas (2 files), xarray (2 files)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
60 files
At the source: github.com/brainets/hoi

andresantoro/EasyBrainSurfacePlots

License: GPL-3.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: c709d995f97c8f8e02055f7d2862b6aa0606ea89, 22 June 2023
Languages: Jupyter (1), Python (1)
Size: 10 files, 2 scripts
Software Heritage: not archived
Found in: the references
Holds: README, license file, 1 notebook
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Matplotlib (2 files), NumPy (2 files), BrainSpace (1 file), neuromaps (1 file), SciPy (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/andresantoro/EasyBrainSurfacePlots

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-75959-w.

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:

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

  • no repository, dataset or request procedure was recognized in it

Read it in the paper: doi.org/10.1038/s41467-026-75959-w.

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, 7 authors, 2 keywords, 10 MeSH terms, 2 funders, 154 references.

Cite

This paper

Santoro, A., Neri, M., Poetto, S., Orsenigo, D., Diano, M., Gatica, M., & Petri, G. (2026). Charting higher-order models of brain function beyond pairwise interactions. Nature communications, 17(1), 9207. https://doi.org/10.1038/s41467-026-75959-w

BibTeX

@article{santoro2026charting,
author = {Santoro, Andrea and Neri, Matteo and Poetto, Simone and Orsenigo, Davide and Diano, Matteo and Gatica, Marilyn and Petri, Giovanni},
title = {{Charting higher-order models of brain function beyond pairwise interactions}},
journal = {Nature communications},
year = {2026},
month = jul,
volume = {17},
number = {1},
pages = {9207},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-75959-w},
url = {https://doi.org/10.1038/s41467-026-75959-w},
pmid = {42660887},
pmcid = {PMC13522477}
}

RIS

TY - JOUR
AU - Santoro, Andrea
AU - Neri, Matteo
AU - Poetto, Simone
AU - Orsenigo, Davide
AU - Diano, Matteo
AU - Gatica, Marilyn
AU - Petri, Giovanni
TI - Charting higher-order models of brain function beyond pairwise interactions
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/07/29
VL - 17
IS - 1
SP - 9207
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-75959-w
UR - https://doi.org/10.1038/s41467-026-75959-w
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-75959-w",
"type": "article-journal",
"title": "Charting higher-order models of brain function beyond pairwise interactions",
"container-title": "Nature communications",
"author": [
{
"family": "Santoro",
"given": "Andrea"
},
{
"family": "Neri",
"given": "Matteo"
},
{
"family": "Poetto",
"given": "Simone"
},
{
"family": "Orsenigo",
"given": "Davide"
},
{
"family": "Diano",
"given": "Matteo"
},
{
"family": "Gatica",
"given": "Marilyn"
},
{
"family": "Petri",
"given": "Giovanni"
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "9207",
"DOI": "10.1038/s41467-026-75959-w",
"PMID": "42660887",
"PMCID": "PMC13522477",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-75959-w",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
29
]
]
}
}

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.1016/j.patter.2026.101619 [code]
Sampling bias corrections for discrete and Gaussian partial information decompositions.
Journal: Patterns (New York, N.Y.)
In common: Nilearn, NiBabel, Statistics and Machine Learning Toolbox, 6 other tools, 17 references, author Davide Orsenigo
[2] doi:10.1371/journal.pone.0348005 [code]
THOI: An efficient and accessible library for computing higher-order interactions enhanced by batch-processing.
Journal: PloS one
In common: NetworkX, seaborn, scikit-learn, 4 other tools, 15 references
[3] doi:10.1038/s41593-026-02205-3 [code]
Competitive interactions shape mammalian brain network dynamics and computation.
Journal: Nature neuroscience
In common: NetworkX, NiBabel, Statistics and Machine Learning Toolbox, 6 other tools, 15 references
[4] doi:10.1038/s41531-026-01354-3 [code]
Neuromodulation-induced normalization of cortical metastable dynamics signatures in Parkinson's disease.
Journal: NPJ Parkinson's disease
In common: neuromaps, BrainSpace, Nilearn, 9 other tools, 9 references
[5] doi:10.1002/hbm.70485 [code]
Exploring the Role of the Rich Club in Network Control of Neurocognitive States.
Journal: Human brain mapping
In common: BrainSpace, Nilearn, NiBabel, 5 other tools, 12 references
[6] doi:10.1038/s41467-026-74466-2 [code]
Neuromorphic hierarchical modular reservoirs.
Journal: Nature communications
In common: neuromaps, Numba, Nilearn, 9 other tools, 8 references
[7] doi:10.1038/s41467-026-71270-w [code]
Spatiotemporal dynamics of the human cortical functional hierarchy across the lifespan.
Journal: Nature communications
In common: BrainSpace, Nilearn, NetworkX, 9 other tools, 9 references
[8] doi:10.1038/s41398-026-04025-2 [code]
Brain energetic landscapes shape state dysregulation in major depressive disorder: a morphological network controllability perspective.
Journal: Translational psychiatry
In common: neuromaps, BrainSpace, Nilearn, 9 other tools, 8 references
[9] 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: neuromaps, Nilearn, statsmodels, 7 other tools, 12 references
[10] doi:10.21203/rs.3.rs-9326213/v1 [code]
Multi-task fMRI outperforms resting-state fMRI for revealing task-invariant organization of the human brain
Journal: Research Square (preprint)
In common: neuromaps, Numba, Nilearn, 8 other tools, 9 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.