OSCR

Focused ultrasound blood-brain barrier opening reveals a paradoxical remote metabolic response in the primate brain.

Code ↔ Paper

3 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 3 matches
  1. [1] § MATERIALS AND METHODS › Statistical analysis ↔ 3-Combined_Monkeys_OEF_Bayesian_template.ipynb, lines 641–699 · score 0.77 · random intercept, PyMC, ROI side, Student, likelihood, hierarchical
  2. [2] § MATERIALS AND METHODS › Quantitative BOLD modeling ↔ 1-Monkey_qBOLD_BBBO_template.ipynb, lines 96–102 · score 0.62 · PyTorch, model parameters, OEF map, fitting
  3. [3] § MATERIALS AND METHODS › LIFU BBBO ↔ 3-Combined_Monkeys_OEF_Bayesian_template.ipynb, lines 18–100 · score 0.55 · right caudate, left caudate, M2, M1, M3, M4

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 · 951 lines · 48 KB · apgl-v3 · 2 matches

  1. # %%
  2. # -----------------------------------------------------------------------------
  3. # 0) IMPORTS
  4. # -----------------------------------------------------------------------------
  5. import os
  6. from pathlib import Path
  7. import numpy as np
  8. import pandas as pd
  9. import nibabel as nib
  10. import matplotlib.pyplot as plt
  11. import seaborn as sns
  12. import matplotlib.gridspec as gridspec
  13. import jax
  14. import pymc as pm
  15. import arviz as az
  16. from IPython.display import display
  18. # -----------------------------------------------------------------------------
  19. # 1) SETTINGS & PATHS for MULTI‐SUBJECT ANALYSIS
  20. # -----------------------------------------------------------------------------
  21. # Set the base directory for all processed data here. Change this to update all paths.
  22. BASE_DATA_DIR = Path("YOUR DIRECTORY")
  23. output_dir = BASE_DATA_DIR / "CombinedMonkeys_OEF_Bayesian"
  24. output_dir.mkdir(parents=True, exist_ok=True) # Ensure the directory exists
  26. subjects = {
  27. 'M4': dict(
  28. SUBFOLDER = "M4",
  29. VOLUMES = {
  30. "Base": ["1_OEF_Base.nii.gz"],
  31. "RCN": ["2_OEF_BBBO_RCN.nii.gz"],
  32. "LCN": ["3_OEF_BBBO_LCN.nii.gz"]
  33. },
  34. MASK_BRAIN = "BrainMask_common.nii.gz",
  35. MASK_ROI_RCN = "2_FUSROI_BBBO_RCN.nii.gz",
  36. MASK_ROI_LCN = "3_FUSROI_BBBO_LCN.nii.gz",
  37. MASK_Caudate_Right = "RightCaudate.nii.gz",
  38. MASK_Caudate_Left = "LeftCaudate.nii.gz",
  39. MASK_Putamen_Right = "RightPutamen.nii.gz",
  40. MASK_Putamen_Left = "LeftPutamen.nii.gz"
  41. ),
  42. 'M3': dict(
  43. SUBFOLDER = "M3",
  44. VOLUMES = {
  45. "Base": ["1_OEF_Base.nii.gz"], # Ctrl1
  46. "RCN": [
  47. "2_OEF_BBBO_RCN.nii.gz", # RCN3
  48. "3_OEF_BBBO_RCN.nii.gz" # RCN4
  49. ],
  50. "LCN": ["4_OEF_BBBO_LCN.nii.gz"] # LCN
  51. },
  52. MASK_BRAIN = "BrainMask_common.nii.gz",
  53. MASK_ROI_RCN = [
  54. "2_FUSROI_BBBO_RCN.nii.gz", # FUSROI for RCN3
  55. "3_FUSROI_BBBO_RCN.nii.gz" # FUSROI for RCN4
  56. ],
  57. MASK_ROI_LCN = "4_FUSROI_BBBO_LCN.nii.gz",
  58. MASK_Caudate_Right = "RightCaudate.nii.gz",
  59. MASK_Caudate_Left = "LeftCaudate.nii.gz",
  60. MASK_Putamen_Right = "RightPutamen.nii.gz",
  61. MASK_Putamen_Left = "LeftPutamen.nii.gz"
  62. ),
  63. 'M1': dict(
  64. SUBFOLDER = "M1",
  65. VOLUMES = {
  66. "Base": [
  67. "1_OEF_Base.nii.gz",
  68. "2_OEF_Base.nii.gz"
  69. ],
  70. "RCN": ["3_OEF_BBBO_RCN.nii.gz"]
  71. },
  72. MASK_BRAIN = "BrainMask_common.nii.gz",
  73. MASK_ROI_RCN = "3_FUSROI_BBBO_RCN.nii.gz",
  74. MASK_Caudate_Right = "RightCaudate.nii.gz",
  75. MASK_Caudate_Left = "LeftCaudate.nii.gz",
  76. MASK_Putamen_Right = "RightPutamen.nii.gz",
  77. MASK_Putamen_Left = "LeftPutamen.nii.gz"
  78. ),
  79. 'M2': dict(
  80. SUBFOLDER = "M2",
  81. VOLUMES = {
  82. "Base": ["1_OEF_Base.nii.gz"],
  83. "RCN": [
  84. "2_OEF_BBBO_RCN.nii.gz", # RCN2
  85. "3_OEF_BBBO_RCN.nii.gz" # RCN3
  86. ]
  87. },
  88. MASK_BRAIN = "BrainMask_common.nii.gz",
  89. MASK_ROI_RCN = [
  90. "2_FUSROI_BBBO_RCN.nii.gz", # FUSROI for RCN3
  91. "3_FUSROI_BBBO_RCN.nii.gz" # FUSROI for RCN4
  92. ],
  93. MASK_Caudate_Right = "RightCaudate.nii.gz",
  94. MASK_Caudate_Left = "LeftCaudate.nii.gz",
  95. MASK_Putamen_Right = "RightPutamen.nii.gz",
  96. MASK_Putamen_Left = "LeftPutamen.nii.gz"
  97. ),
  98. }
  100. TREATMENT_TYPES = ["RCN","LCN"] # Define all possible treatment types
  102. # seaborn style overrides
  103. custom_params = {
  104. 'ytick.left': True, 'xtick.bottom': True,
  105. 'xtick.direction': 'out','xtick.color': 'black',
  106. 'ytick.color': 'black','text.color': 'black',
  107. 'grid.color': 'black'
  108. }
  109. sns.set_theme(
  110. context='notebook', style='white', palette='deep',
  111. font='sans-serif', font_scale=1, color_codes=True,
  112. rc=custom_params
  113. )
  115. # -----------------------------------------------------------------------------
  116. # 2) HELPER FUNCTIONS (bootstrap_diff, cohens_d, load_nifti are unchanged)
  117. # -----------------------------------------------------------------------------
  118. def load_nifti(path):
  119. return nib.load(str(path)).get_fdata()
  121. def bootstrap_diff(T, B, n_boot=20000, ci=95):
  122. T_flat = np.asarray(T).ravel()
  123. B_flat = np.asarray(B).ravel()
  124. T_clean = T_flat[~np.isnan(T_flat)]
  125. B_clean = B_flat[~np.isnan(B_flat)]
  126. if T_clean.size < 2 or B_clean.size < 2 : # Need at least 2 points for meaningful bootstrap
  127. return np.full(n_boot, np.nan), np.nan, np.nan, np.nan
  128. boot_T_means = np.random.choice(T_clean, (n_boot, T_clean.size), replace=True).mean(axis=1)
  129. boot_B_means = np.random.choice(B_clean, (n_boot, B_clean.size), replace=True).mean(axis=1)
  130. boot_diff_means = boot_T_means - boot_B_means
  131. if np.all(np.isnan(boot_diff_means)):
  132. return boot_diff_means, np.nan, np.nan, np.nan
  133. lo, hi = np.nanpercentile(boot_diff_means, [(100-ci)/2, 100-(100-ci)/2])
  134. p_val_numerator = np.nansum(boot_diff_means <= 0) if np.nanmean(boot_diff_means) > 0 else np.nansum(boot_diff_means >= 0)
  135. p_val_denominator = np.sum(~np.isnan(boot_diff_means))
  136. if p_val_denominator == 0: p = np.nan
  137. else: p = 2 * (p_val_numerator / p_val_denominator); p = min(p, 1.0)
  138. return boot_diff_means, lo, hi, p
  140. def cohens_d(x, y):
  141. x_arr = np.asarray(x).ravel(); y_arr = np.asarray(y).ravel()
  142. x_arr = x_arr[~np.isnan(x_arr)]; y_arr = y_arr[~np.isnan(y_arr)]
  143. nx, ny = len(x_arr), len(y_arr)
  144. if nx < 2 or ny < 2: return np.nan
  145. vx, vy = x_arr.var(ddof=1), y_arr.var(ddof=1)
  146. if vx == 0 and vy == 0: return 0.0 if x_arr.mean() == y_arr.mean() else np.inf * np.sign(x_arr.mean() - y_arr.mean())
  147. pooled_sd_numerator = (nx-1)*vx + (ny-1)*vy
  148. pooled_sd_denominator = nx+ny-2
  149. if pooled_sd_denominator == 0 : return np.nan
  150. pooled_sd = np.sqrt(pooled_sd_numerator / pooled_sd_denominator)
  151. if pooled_sd == 0: return np.nan
  152. return (x_arr.mean() - y_arr.mean()) / pooled_sd
  154. # -----------------------------------------------------------------------------
  155. # 3–7) LOOP OVER SUBJECTS → BUILD DATAFRAMES
  156. # -----------------------------------------------------------------------------
  157. all_df_roi = []
  158. all_df_ic = []
  159. all_df_hemi_ic = []
  160. all_df_sum = []
  162. for subj, cfg in subjects.items():
  163. print(f"Processing subject: {subj}")
  164. brain_mask_path = BASE_DATA_DIR / cfg['SUBFOLDER'] / cfg['MASK_BRAIN']
  165. if not brain_mask_path.exists():
  166. print(f"Critical: Brain mask not found for {subj} at {brain_mask_path}. Skipping subject.")
  167. continue
  168. brain_mask = load_nifti(brain_mask_path) > 0
  169. nx_shape = brain_mask.shape[0] # Get shape for flipping ROIs
  171. # Load structural masks
  172. mask_caudate_r = load_nifti(BASE_DATA_DIR / cfg['SUBFOLDER'] / cfg['MASK_Caudate_Right']) > 0 if (BASE_DATA_DIR / cfg['SUBFOLDER'] / cfg['MASK_Caudate_Right']).exists() else np.zeros(brain_mask.shape, dtype=bool)
  173. mask_caudate_l = load_nifti(BASE_DATA_DIR / cfg['SUBFOLDER'] / cfg['MASK_Caudate_Left']) > 0 if (BASE_DATA_DIR / cfg['SUBFOLDER'] / cfg['MASK_Caudate_Left']).exists() else np.zeros(brain_mask.shape, dtype=bool)
  174. mask_putamen_r = load_nifti(BASE_DATA_DIR / cfg['SUBFOLDER'] / cfg['MASK_Putamen_Right']) > 0 if (BASE_DATA_DIR / cfg['SUBFOLDER'] / cfg['MASK_Putamen_Right']).exists() else np.zeros(brain_mask.shape, dtype=bool)
  175. mask_putamen_l = load_nifti(BASE_DATA_DIR / cfg['SUBFOLDER'] / cfg['MASK_Putamen_Left']) > 0 if (BASE_DATA_DIR / cfg['SUBFOLDER'] / cfg['MASK_Putamen_Left']).exists() else np.zeros(brain_mask.shape, dtype=bool)
  177. # Load and z-score all volumes
  178. zvols_all_sessions = {}
  179. raw_vols_all_sessions = {}
  181. for cond_type, file_list in cfg['VOLUMES'].items():
  182. if not isinstance(file_list, list):
  183. print(f"Warning: VOLUMES for {cond_type} in {subj} is not a list. Skipping this condition type.")
  184. continue
  185. for idx, fname in enumerate(file_list):
  186. session_name = f"{cond_type}{idx+1}" if len(file_list) > 1 or cond_type in TREATMENT_TYPES else cond_type
  187. if len(file_list) == 1 and cond_type == "Base": session_name = "Base"
  189. vol_path = BASE_DATA_DIR / cfg['SUBFOLDER'] / fname
  190. if not vol_path.exists():
  191. print(f"Warning: Volume file not found for {subj}, {session_name}: {vol_path}. Skipping this volume.")
  192. continue
  194. vol_data = load_nifti(vol_path)
  195. raw_vols_all_sessions[session_name] = vol_data
  197. vol_brain_voxels = vol_data[brain_mask & ~np.isnan(vol_data)]
  198. if vol_brain_voxels.size > 1:
  199. μ, σ = vol_brain_voxels.mean(), vol_brain_voxels.std()
  200. if σ != 0:
  201. zvols_all_sessions[session_name] = (vol_data - μ) / σ
  202. else:
  203. zvols_all_sessions[session_name] = np.zeros_like(vol_data) if np.all(vol_brain_voxels == μ) else np.full_like(vol_data, np.nan)
  204. else:
  205. zvols_all_sessions[session_name] = np.full_like(vol_data, np.nan)
  206. print(f"Warning: Not enough brain voxels for z-scoring {subj}, {session_name}.")
  208. base_zvol_keys = [k for k in zvols_all_sessions if k.startswith("Base")]
  209. if not base_zvol_keys:
  210. print(f"Critical: No valid Base volumes processed for subject {subj}. Skipping subject.")
  211. continue
  213. mean_base_z_volume = None
  214. valid_base_zvols_list = [zvols_all_sessions[k] for k in base_zvol_keys if k in zvols_all_sessions and not np.all(np.isnan(zvols_all_sessions[k]))]
  216. if not valid_base_zvols_list:
  217. print(f"Critical: All base volumes for {subj} are NaN or missing. Cannot compute mean base. Skipping subject.")
  218. continue
  219. if len(valid_base_zvols_list) > 1:
  220. mean_base_z_volume = np.nanmean(np.stack(valid_base_zvols_list, axis=-1), axis=-1)
  221. else:
  222. mean_base_z_volume = valid_base_zvols_list[0]
  224. if mean_base_z_volume is None or np.all(np.isnan(mean_base_z_volume)):
  225. print(f"Critical: Mean base z-volume is NaN for subject {subj}. Skipping subject.")
  226. continue
  228. mid = nx_shape // 2
  229. mask_hemi_r = brain_mask.copy(); mask_hemi_r[:mid,...] = False; mask_hemi_r &= brain_mask
  230. mask_hemi_l = brain_mask.copy(); mask_hemi_l[mid:,...] = False; mask_hemi_l &= brain_mask
  232. treatment_session_keys = [k for k in zvols_all_sessions if any(k.startswith(tt) for tt in TREATMENT_TYPES) and k in zvols_all_sessions]
  234. subj_roi_rows = []
  235. subj_ic_rows = []
  236. subj_hemi_ic_rows = []
  237. subj_sum_rows = []
  239. for treat_session_key in treatment_session_keys:
  240. treat_type_prefix = next((tt for tt in TREATMENT_TYPES if treat_session_key.startswith(tt)), None)
  241. if treat_type_prefix is None: continue
  243. zvol_treat = zvols_all_sessions[treat_session_key]
  244. if np.all(np.isnan(zvol_treat)):
  245. print(f"Warning: Treatment session {treat_session_key} for {subj} is all NaN. Skipping this session's analysis.")
  246. continue
  248. session_specific_rois = {}
  249. base_roi_data_ipsi = None
  250. mask_fname_to_load_str = None
  252. mask_key_generic = f"MASK_ROI_{treat_type_prefix}"
  253. if mask_key_generic in cfg and cfg[mask_key_generic] is not None:
  254. mask_value = cfg[mask_key_generic]
  255. # Always treat as list
  256. if not isinstance(mask_value, list):
  257. mask_value = [mask_value]
  258. # Determine session index
  259. try:
  260. session_idx = int(treat_session_key.replace(treat_type_prefix, '')) - 1
  261. except ValueError:
  262. session_idx = 0 # Default to first if cannot parse
  263. if 0 <= session_idx < len(mask_value):
  264. mask_fname_to_load_str = mask_value[session_idx]
  265. else:
  266. print(f"Warning: Index {session_idx} out of bounds for {subj}'s {mask_key_generic} (len {len(mask_value)}). Session: {treat_session_key}")
  267. else:
  268. print(f"Notice: No ROI mask key '{mask_key_generic}' defined for {subj}, session {treat_session_key}.")
  270. if mask_fname_to_load_str:
  271. mask_path_to_load = BASE_DATA_DIR / cfg['SUBFOLDER'] / mask_fname_to_load_str
  272. if mask_path_to_load.exists():
  273. loaded_mask = load_nifti(mask_path_to_load) > 0
  274. if np.any(loaded_mask):
  275. base_roi_data_ipsi = loaded_mask & brain_mask
  276. else:
  277. print(f"Warning: Loaded ROI mask is empty for {subj}, {treat_session_key}: {mask_path_to_load}")
  278. else:
  279. print(f"Warning: ROI mask file not found for {subj}, {treat_session_key}: {mask_path_to_load}")
  281. if base_roi_data_ipsi is not None and np.any(base_roi_data_ipsi):
  282. if treat_type_prefix == "LCN":
  283. session_specific_rois[f"{treat_type_prefix} L"] = base_roi_data_ipsi
  284. session_specific_rois[f"{treat_type_prefix} R"] = np.flip(base_roi_data_ipsi, axis=0) & brain_mask
  285. else:
  286. session_specific_rois[f"{treat_type_prefix} R"] = base_roi_data_ipsi
  287. session_specific_rois[f"{treat_type_prefix} L"] = np.flip(base_roi_data_ipsi, axis=0) & brain_mask
  288. else:
  289. if mask_fname_to_load_str :
  290. print(f"Notice: No valid ROI data for {subj}, session {treat_session_key}. ROI-specific analyses will be skipped for this session.")
  292. for hemi_label, hemi_mask in [('Left', mask_hemi_l), ('Right', mask_hemi_r)]:
  293. if not np.any(hemi_mask): continue
  295. dz_vals_hemi = (zvol_treat[hemi_mask] - mean_base_z_volume[hemi_mask]).ravel()
  296. dz_vals_hemi = dz_vals_hemi[~np.isnan(dz_vals_hemi)]
  297. if dz_vals_hemi.size == 0: continue
  299. is_ipsi_hemi = ((treat_type_prefix in ["RCN", "RV"] and hemi_label=='Right') or \
  300. (treat_type_prefix == "LCN" and hemi_label=='Left'))
  301. ipsicontra_label_hemi = 'Ipsilateral' if is_ipsi_hemi else 'Contralateral'
  303. for v_dz_hemi in dz_vals_hemi:
  304. subj_hemi_ic_rows.append({
  305. 'Subject': subj, 'Session': treat_session_key,
  306. 'FUS_Type': treat_type_prefix, 'Hemisphere': hemi_label,
  307. 'IpsiContra': ipsicontra_label_hemi, 'Δz': float(v_dz_hemi)
  308. })
  310. rois_for_session_analysis = {}
  311. for roi_label, roi_mask in session_specific_rois.items():
  312. rois_for_session_analysis[roi_label] = roi_mask
  313. roi_side = roi_label.split()[-1]
  314. caudate_mask_for_side = mask_caudate_r if roi_side == 'R' else mask_caudate_l
  315. putamen_mask_for_side = mask_putamen_r if roi_side == 'R' else mask_putamen_l
  316. rois_for_session_analysis[f"{roi_label} in Caudate"] = roi_mask & caudate_mask_for_side
  317. rois_for_session_analysis[f"{roi_label} in Putamen"] = roi_mask & putamen_mask_for_side
  318. rois_for_session_analysis[f"{roi_label} outside C&P"] = roi_mask & ~caudate_mask_for_side & ~putamen_mask_for_side
  320. for analysis_roi_label, analysis_roi_mask in rois_for_session_analysis.items():
  321. if not np.any(analysis_roi_mask): continue
  323. base_vals_in_roi = mean_base_z_volume[analysis_roi_mask].ravel()
  324. base_vals_in_roi = base_vals_in_roi[~np.isnan(base_vals_in_roi)]
  325. for v_base_roi in base_vals_in_roi:
  326. subj_roi_rows.append({
  327. 'Subject': subj, 'ROI': analysis_roi_label,
  328. 'Session': "MeanBase", 'FUS_Type': "Base", 'z': float(v_base_roi)
  329. })
  330. treat_vals_in_roi = zvol_treat[analysis_roi_mask].ravel()
  331. treat_vals_in_roi = treat_vals_in_roi[~np.isnan(treat_vals_in_roi)]
  332. for v_treat_roi in treat_vals_in_roi:
  333. subj_roi_rows.append({
  334. 'Subject': subj, 'ROI': analysis_roi_label,
  335. 'Session': treat_session_key, 'FUS_Type': treat_type_prefix, 'z': float(v_treat_roi)
  336. })
  338. dz_vals_roi = (zvol_treat[analysis_roi_mask] - mean_base_z_volume[analysis_roi_mask]).ravel()
  339. dz_vals_roi_clean = dz_vals_roi[~np.isnan(dz_vals_roi)]
  341. if dz_vals_roi_clean.size > 0:
  342. roi_side_char = analysis_roi_label.split()[1]
  343. is_ipsi_roi = ((treat_type_prefix in ["RCN", "RV"] and roi_side_char == 'R') or \
  344. (treat_type_prefix == "LCN" and roi_side_char == 'L'))
  345. ipsicontra_label_roi = 'Ipsilateral' if is_ipsi_roi else 'Contralateral'
  347. for v_dz_roi in dz_vals_roi_clean:
  348. subj_ic_rows.append({
  349. 'Subject': subj, 'Session': treat_session_key, 'FUS_Type': treat_type_prefix,
  350. 'ROI_Class': analysis_roi_label, 'Condition': ipsicontra_label_roi, 'Δz': float(v_dz_roi)
  351. })
  353. T_sum, B_sum = treat_vals_in_roi, base_vals_in_roi
  354. if T_sum.size >= 2 and B_sum.size >= 2:
  355. boot_sum, lo_sum, hi_sum, p_sum = bootstrap_diff(T_sum, B_sum)
  356. d_sum = cohens_d(T_sum, B_sum)
  357. subj_sum_rows.append({
  358. 'Subject': subj, 'Session': treat_session_key, 'FUS_Type': treat_type_prefix,
  359. 'ROI': analysis_roi_label, 'Comparison': f"{treat_session_key} vs MeanBase",
  360. 'Δz': T_sum.mean() - B_sum.mean(), 'CI lower': lo_sum, 'CI upper': hi_sum,
  361. 'p-value': p_sum, "Cohen's d": d_sum
  362. })
  363. else:
  364. subj_sum_rows.append({
  365. 'Subject': subj, 'Session': treat_session_key, 'FUS_Type': treat_type_prefix,
  366. 'ROI': analysis_roi_label, 'Comparison': f"{treat_session_key} vs MeanBase",
  367. 'Δz': np.nanmean(T_sum) - np.nanmean(B_sum) if T_sum.size > 0 and B_sum.size > 0 else np.nan,
  368. 'CI lower': np.nan, 'CI upper': np.nan, 'p-value': np.nan, "Cohen's d": np.nan
  369. })
  371. if subj_roi_rows: all_df_roi.append(pd.DataFrame(subj_roi_rows))
  372. if subj_ic_rows: all_df_ic.append(pd.DataFrame(subj_ic_rows))
  373. if subj_hemi_ic_rows: all_df_hemi_ic.append(pd.DataFrame(subj_hemi_ic_rows))
  374. if subj_sum_rows: all_df_sum.append(pd.DataFrame(subj_sum_rows).set_index(['Subject','Session','ROI','Comparison']))
  376. df_roi = pd.concat(all_df_roi, ignore_index=True) if all_df_roi else pd.DataFrame()
  377. if not df_roi.empty:
  378. # Ensure unique entries for MeanBase per ROI for each subject.
  379. mean_base_entries = df_roi[df_roi['Session'] == 'MeanBase']
  380. other_entries = df_roi[df_roi['Session'] != 'MeanBase']
  381. if not mean_base_entries.empty:
  382. mean_base_entries_dedup = mean_base_entries.drop_duplicates(subset=['Subject', 'ROI', 'Session', 'FUS_Type', 'z'])
  383. df_roi = pd.concat([mean_base_entries_dedup, other_entries], ignore_index=True)
  385. df_ic = pd.concat(all_df_ic, ignore_index=True) if all_df_ic else pd.DataFrame()
  386. df_hemi_ic = pd.concat(all_df_hemi_ic, ignore_index=True) if all_df_hemi_ic else pd.DataFrame()
  387. df_sum = pd.concat(all_df_sum) if all_df_sum else pd.DataFrame() # Per-session summary
  389. # -----------------------------------------------------------------------------
  390. # BUILD REGRESSION DESIGN DATAFRAME (SESSION-LEVEL ROI MEAN Δz)
  391. # -----------------------------------------------------------------------------
  392. if not df_sum.empty:
  393. df_sum_reset = df_sum.reset_index() if isinstance(df_sum, pd.DataFrame) else pd.DataFrame()
  394. df_reg = df_sum_reset.dropna(subset=['Δz', 'Subject', 'Session', 'FUS_Type', 'ROI'])
  395. df_reg = df_reg[df_reg['FUS_Type'].isin(TREATMENT_TYPES)].copy()
  396. df_reg['roi_side'] = df_reg['ROI'].str.split().str[1].map({'R': 1, 'L': -1})
  397. df_reg['fus_side'] = df_reg['FUS_Type'].map({'RCN': 1, 'LCN': -1})
  398. df_reg['Contralateral_Effect'] = df_reg['roi_side'] * df_reg['fus_side'] * (-1)
  399. df_reg['subject_cat'] = pd.Categorical(df_reg['Subject']).codes
  401. # -----------------------------------------------------------------------------
  402. # 8) COMBINED PLOT: SESSION-LEVEL ROI & HEMISPHERIC MEAN Δz BY FUS TYPE
  403. # -----------------------------------------------------------------------------
  404. if not df_sum.empty and not df_hemi_ic.empty:
  405. # ROI effects (already session-level in df_reg)
  406. plot_roi_full = df_reg[~df_reg['ROI'].str.contains(" in | outside ")]
  407. plot_roi_full['Laterality'] = (plot_roi_full['Contralateral_Effect']*(-1)).map({1: 'Ipsilateral', -1: 'Contralateral'})
  408. plot_roi_full['PlotGroup'] = 'ROI ' + plot_roi_full['FUS_Type']
  410. # Hemispheric effects: aggregate mean Δz per subject/session/hemisphere/ipsi-contra
  411. hemi_means = (
  412. df_hemi_ic
  413. .groupby(['Subject', 'Session', 'FUS_Type', 'Hemisphere', 'IpsiContra'])
  414. .agg({'Δz': 'mean'})
  415. .reset_index()
  416. )
  417. # Only keep treatment sessions
  418. hemi_means = hemi_means[hemi_means['FUS_Type'].isin(TREATMENT_TYPES)].copy()
  419. hemi_means['PlotGroup'] = 'Hemisphere ' + hemi_means['FUS_Type']
  420. hemi_means = hemi_means.rename(columns={'IpsiContra': 'Laterality'})
  422. # Match columns for concatenation
  423. plot_hemi = hemi_means[['Subject', 'Session', 'FUS_Type', 'Δz', 'Laterality', 'PlotGroup']].copy()
  424. plot_roi_full = plot_roi_full[['Subject', 'Session', 'FUS_Type', 'Δz', 'Laterality', 'PlotGroup']].copy()
  426. # Combine
  427. plot_df = pd.concat([plot_hemi, plot_roi_full], ignore_index=True)
  429. # Set plotting order
  430. plot_group_order = [
  431. "Hemisphere RCN", # leftmost
  432. "ROI RCN",
  433. "ROI LCN",
  434. "Hemisphere LCN" # rightmost
  435. ]
  436. plot_df = plot_df[plot_df['PlotGroup'].isin(plot_group_order)]
  437. plot_df['PlotGroup'] = pd.Categorical(plot_df['PlotGroup'], categories=plot_group_order, ordered=True)
  438. plot_df['Laterality'] = pd.Categorical(plot_df['Laterality'], categories=['Ipsilateral', 'Contralateral'], ordered=True)
  440. # Calculate n for each group/lat for legend
  441. n_dict = {}
  442. for group in plot_group_order:
  443. for lat in ['Ipsilateral', 'Contralateral']:
  444. mask = (plot_df['PlotGroup'] == group) & (plot_df['Laterality'] == lat)
  445. n = mask.sum()
  446. n_dict[(group, lat)] = n
  447. # Build legend labels with n in parentheses
  448. legend_labels = []
  449. for group in plot_group_order:
  450. n_ipsi = n_dict[(group, 'Ipsilateral')]
  451. n_contra = n_dict[(group, 'Contralateral')]
  452. legend_labels.append(f"{group} (n={n_ipsi})")
  454. fig, ax = plt.subplots(figsize=(6, 7))
  456. dodge_offsets = {
  457. "Hemisphere RCN": -0.375,
  458. "ROI RCN": -0.125,
  459. "ROI LCN": 0.125,
  460. "Hemisphere LCN": 0.375
  461. }
  462. # Define marker styles for each group
  463. marker_map = {
  464. "Hemisphere RCN": "X",
  465. "ROI RCN": "o",
  466. "ROI LCN": "s",
  467. "Hemisphere LCN": "^"
  468. }
  470. # Separate out groups with n=2
  471. plot_df_for_pointplot = plot_df.copy()
  472. for group in ["ROI LCN", "Hemisphere LCN"]:
  473. for lat in ['Ipsilateral', 'Contralateral']:
  474. mask = (plot_df['PlotGroup'] == group) & (plot_df['Laterality'] == lat)
  475. sub_df = plot_df[mask]
  476. if len(sub_df) == 2:
  477. x_val = 0 if lat == 'Ipsilateral' else 1
  478. x_val_dodged = x_val + dodge_offsets[group]
  479. # Plot each point individually at the correct x
  480. for y in sub_df['Δz']:
  481. ax.scatter(x_val_dodged, y, marker=marker_map[group], color='black', s=80, zorder=3)
  482. plot_df_for_pointplot = plot_df_for_pointplot[~mask]
  484. # Plot the rest with error bars
  485. if not plot_df_for_pointplot.empty:
  486. sns.pointplot(
  487. data=plot_df_for_pointplot,
  488. x='Laterality',
  489. y='Δz',
  490. hue='PlotGroup',
  491. hue_order=plot_group_order,
  492. markers=['X', 'o', 's', '^'],
  493. palette=['black', 'black', 'black', 'black'],
  494. linestyles=['none'] * 4,
  495. dodge=0.5,
  496. errorbar=('ci', 95),
  497. capsize=0.05,
  498. err_kws={'linewidth': 1.5, 'alpha': 0.7},
  499. ax=ax,
  500. markersize=8
  501. )
  503. ax.axhline(0, ls='--', color='darkgray', zorder=0)
  504. ax.set_title("", fontsize=24)
  505. ax.set_ylabel(r"$\overline{\Delta z}_{\text{OEF}}$ (Treatment - Control)", fontsize=24)
  506. ax.set_xlabel("", fontsize=24)
  507. ax.set_ylim(-0.4, 0.3)
  508. ax.set_yticks([-0.4, 0, 0.3])
  509. ax.set_yticklabels([str(x) for x in [-0.4, 0, 0.3]], fontsize=22)
  510. #ax.tick_params(axis='both', which='major', labelsize=22)
  511. ax.set_xticklabels(["Ipsilateral", "Contralateral"], fontsize=24)
  512. # Set legend with n in parentheses
  513. handles, _ = ax.get_legend_handles_labels()
  514. #ax.legend(handles, legend_labels, title="", loc='upper left', bbox_to_anchor=(1.02, 1), ncol=1, fontsize=10)
  515. ax.legend(handles, legend_labels, title="", loc='lower right', ncol=1, fontsize=20, framealpha=1, edgecolor='black')
  516. plt.tight_layout()
  517. plt.savefig(output_dir / "SESSION_LEVEL_ROI_HEMISPHERIC_MEAN_Delta_z.svg", format='svg')
  518. plt.show()
  519. plt.close(fig)
  520. else:
  521. print("Skipping Combined Plot: No session-level ROI or hemispheric mean Δz data available.")
  523. # -----------------------------------------------------------------------------
  524. # 8a) PLOTS FOR CAUDATE & PUTAMEN SUB-REGIONS
  525. # -----------------------------------------------------------------------------
  526. def plot_subregion(df, region_name):
  527. df_plot = df[df['ROI'].str.contains(f" in {region_name}")].copy()
  528. if not df_plot.empty:
  529. df_plot['Laterality'] = (df_plot['Contralateral_Effect']*(-1)).map({1: 'Ipsilateral', -1: 'Contralateral'})
  530. # Calculate n for each FUS_Type/Laterality for legend
  531. n_dict = {}
  532. for fus_type in ['RCN', 'LCN']:
  533. for lat in ['Ipsilateral', 'Contralateral']:
  534. mask = (df_plot['FUS_Type'] == fus_type) & (df_plot['Laterality'] == lat)
  535. n = mask.sum()
  536. n_dict[(fus_type, lat)] = n
  537. legend_labels = []
  538. for fus_type in ['RCN', 'LCN']:
  539. n_ipsi = n_dict[(fus_type, 'Ipsilateral')]
  540. n_contra = n_dict[(fus_type, 'Contralateral')]
  541. legend_labels.append(f"{fus_type} (n={n_ipsi})")
  542. fig, ax = plt.subplots(figsize=(5, 7))
  544. # --- Custom plotting logic for N=2 ---
  545. marker_map = {'RCN': 'o', 'LCN': 's'}
  546. dodge_offsets = {'RCN': -0.1, 'LCN': 0.1} # Adjust if you change dodge in pointplot
  548. plot_df_for_pointplot = df_plot.copy()
  549. for fus_type in ['RCN', 'LCN']:
  550. for lat in ['Ipsilateral', 'Contralateral']:
  551. mask = (df_plot['FUS_Type'] == fus_type) & (df_plot['Laterality'] == lat)
  552. sub_df = df_plot[mask]
  553. if len(sub_df) == 2:
  554. x_val = 0 if lat == 'Ipsilateral' else 1
  555. x_val_dodged = x_val + dodge_offsets[fus_type]
  556. for y in sub_df['Δz']:
  557. ax.scatter(x_val_dodged, y, marker=marker_map[fus_type], color='black', s=80, zorder=3)
  558. plot_df_for_pointplot = plot_df_for_pointplot[~mask]
  560. # Plot the rest with error bars (N>2) using seaborn pointplot
  561. if not plot_df_for_pointplot.empty:
  562. sns.pointplot(
  563. data=plot_df_for_pointplot,
  564. x='Laterality',
  565. y='Δz',
  566. hue='FUS_Type',
  567. hue_order=['RCN', 'LCN'],
  568. markers=['o', 's'],
  569. palette=['black', 'black'],
  570. linestyles=['none', 'none'],
  571. dodge=0.2,
  572. errorbar=('ci', 95),
  573. err_kws={'linewidth': 1.5, 'alpha': 0.7},
  574. capsize=0.05,
  575. ax=ax,
  576. markersize=8,
  577. )
  578. ax.axhline(0, ls='--', color='darkgray', zorder=0)
  579. ax.set_title(rf"{region_name}", fontsize=24)
  580. ax.set_ylabel(r"$\overline{\Delta z}_{\text{OEF}}$ (Treatment - Control)", fontsize=24)
  581. ax.set_xlabel("", fontsize=24)
  582. ax.set_ylim(-0.4, 0.3)
  583. ax.set_yticks([-0.4, 0, 0.3])
  584. ax.set_yticklabels([str(x) for x in [-0.4, 0, 0.3]], fontsize=22)
  585. ax.set_xticklabels(["Ipsilateral", "Contralateral"], fontsize=24)
  586. handles, _ = ax.get_legend_handles_labels()
  587. ax.legend(handles, legend_labels, title="", loc='lower right', ncol=1, fontsize=20, shadow=False, frameon=False, framealpha=1, edgecolor='black', borderpad=0.4, borderaxespad=0.4)
  588. plt.tight_layout()
  589. if region_name == "Caudate":
  590. plt.savefig(output_dir / "SESSION_LEVEL_Caudate_ROI_HEMISPHERIC_MEAN_Delta_z.svg", format='svg')
  591. elif region_name == "Putamen":
  592. plt.savefig(output_dir / "SESSION_LEVEL_Putamen_ROI_HEMISPHERIC_MEAN_Delta_z.svg", format='svg')
  593. plt.show()
  594. plt.close(fig)
  595. else:
  596. print(f"No data for {region_name} sub-region plot.")
  598. if not df_reg.empty:
  599. plot_subregion(df_reg, "Caudate")
  600. plot_subregion(df_reg, "Putamen")
  602. # -----------------------------------------------------------------------------
  603. # 10) OVERALL SUMMARY TABLE (Δz, CI, p‐value, Cohen’s d) - Aggregated per ROI_Class
  604. # -----------------------------------------------------------------------------
  605. overall_summaries = []
  606. if not df_roi.empty:
  607. df_roi_overall = df_roi.copy()
  608. df_roi_overall['ROI_Class'] = df_roi_overall['ROI']
  609. df_roi_overall['Condition_Type'] = df_roi_overall['FUS_Type']
  611. unique_roi_classes = sorted(df_roi_overall['ROI_Class'].unique())
  613. for roi_cls in unique_roi_classes:
  614. primary_fus_type = roi_cls.split()[0]
  615. if primary_fus_type not in TREATMENT_TYPES : continue
  617. T_all = df_roi_overall.query("ROI_Class == @roi_cls and Condition_Type == @primary_fus_type").z.values
  618. B_all = df_roi_overall.query("ROI_Class == @roi_cls and Condition_Type == 'Base'").z.values
  620. if T_all.size >= 2 and B_all.size >= 2:
  621. # Check for NaNs after selection, as z-scores could be NaN if original data was problematic
  622. T_all_clean = T_all[~np.isnan(T_all)]
  623. B_all_clean = B_all[~np.isnan(B_all)]
  624. if T_all_clean.size >=2 and B_all_clean.size >=2:
  625. boot_ov, lo_ov, hi_ov, p_ov = bootstrap_diff(T_all_clean, B_all_clean)
  626. d_ov = cohens_d(T_all_clean, B_all_clean)
  627. mean_diff_ov = T_all_clean.mean() - B_all_clean.mean()
  628. else:
  629. mean_diff_ov, lo_ov, hi_ov, p_ov, d_ov = np.nan, np.nan, np.nan, np.nan, np.nan
  630. else:
  631. mean_diff_ov, lo_ov, hi_ov, p_ov, d_ov = np.nan, np.nan, np.nan, np.nan, np.nan
  633. overall_summaries.append({
  634. 'ROI_Class': roi_cls,
  635. 'Comparison': f"{primary_fus_type} (pooled) vs Base (pooled)",
  636. 'Δz': mean_diff_ov,
  637. 'CI lower': lo_ov, 'CI upper': hi_ov,
  638. 'p-value': p_ov, "Cohen's d": d_ov
  639. })
  641. # -----------------------------------------------------------------------------
  642. # 12) BAYESIAN REGRESSION (SESSION-LEVEL ROI MEAN Δz, HIERARCHICAL)
  643. # -----------------------------------------------------------------------------
  644. idata = None
  645. df_reg_full_roi = df_reg[~df_reg['ROI'].str.contains(" in | outside ")].copy()
  646. if not df_reg_full_roi.empty and 'Δz' in df_reg_full_roi and not df_reg_full_roi['Δz'].isnull().all() and len(df_reg_full_roi) > 1:
  647. coords = {
  648. "subject": df_reg_full_roi['Subject'].unique(),
  649. "obs_id": np.arange(len(df_reg_full_roi))
  650. }
  651. subject_map = {name: i for i, name in enumerate(df_reg_full_roi['Subject'].unique())}
  652. subject_indices_for_model = df_reg_full_roi['Subject'].map(subject_map).values
  654. with pm.Model(coords=coords) as model_hierarchical:
  655. # Data
  656. delta_z_obs = pm.Data("delta_z_obs", df_reg_full_roi["Δz"].values, dims="obs_id")
  657. fus_side_data = pm.Data("fus_side_data", df_reg_full_roi["fus_side"].values.astype(float), dims="obs_id")
  658. roi_side_data = pm.Data("roi_side_data", df_reg_full_roi["roi_side"].values.astype(float), dims="obs_id")
  659. contralateral_effect_data = pm.Data("contralateral_effect_data", df_reg_full_roi["Contralateral_Effect"].values.astype(float), dims="obs_id")
  660. subject_idx_data = pm.Data("subject_idx_data", subject_indices_for_model, dims="obs_id")
  662. # Hyperpriors for random intercepts
  663. mu_alpha_subj = pm.Normal('mu_alpha_subj', mu=0, sigma=2)
  664. sigma_alpha_subj = pm.HalfNormal('sigma_alpha_subj', sigma=2)
  665. alpha_subj_offset = pm.Normal('alpha_subj_offset', mu=0, sigma=1, dims="subject")
  666. alpha_subj = pm.Deterministic('alpha_subj', mu_alpha_subj + alpha_subj_offset * sigma_alpha_subj, dims="subject")
  668. # Fixed effects
  669. β_fus = pm.Normal('FUS_Target_Side', mu=0, sigma=2)
  670. β_read = pm.Normal('Read_Side', mu=0, sigma=2)
  671. β_int = pm.Normal('Contralateral_Effect', mu=0, sigma=2)
  673. σ_resid = pm.HalfNormal('sigma_resid', sigma=2)
  674. ν_resid = pm.Exponential('nu_resid', 1/30.0)
  676. mu_eq = alpha_subj[subject_idx_data] + \
  677. β_fus * fus_side_data + \
  678. β_read * roi_side_data + \
  679. β_int * contralateral_effect_data
  681. y_obs = pm.StudentT('y', nu=ν_resid, mu=mu_eq, sigma=σ_resid,
  682. observed=delta_z_obs, dims="obs_id")
  684. print("Sampling Bayesian model using CPU (default NUTS sampler)...")
  685. try:
  686. idata = pm.sample(5000, tune=2000, chains=16,
  687. cores=16,
  688. return_inferencedata=True,
  689. idata_kwargs={"log_likelihood":True},
  690. target_accept=0.95,
  691. random_seed=80)
  692. ppc = pm.sample_posterior_predictive(idata, var_names=['y', 'alpha_subj'], random_seed=80, return_inferencedata=True)
  693. idata.extend(ppc)
  694. print("Sampling complete.")
  695. except Exception as e:
  696. print(f"Error during PyMC sampling: {e}")
  697. idata = None
  698. else:
  699. print("Skipping Bayesian Regression: No session-level ROI mean Δz data available.")
  701. # -----------------------------------------------------------------------------
  702. # 12a) BAYESIAN REGRESSION (SUB-REGIONS)
  703. # -----------------------------------------------------------------------------
  704. idata_v2 = None
  705. df_reg_v2 = df_reg[df_reg['ROI'].str.contains(" in | outside ")].copy()
  707. if not df_reg_v2.empty and 'Δz' in df_reg_v2 and not df_reg_v2['Δz'].isnull().all() and len(df_reg_v2) > 1:
  709. def get_subregion_type(roi_name):
  710. if "in Caudate" in roi_name: return "Caudate"
  711. if "in Putamen" in roi_name: return "Putamen"
  712. if "outside C&P" in roi_name: return "Outside"
  713. return "Other"
  714. df_reg_v2['subregion'] = df_reg_v2['ROI'].apply(get_subregion_type)
  716. df_reg_v2['Contralateral_Effect_Caudate'] = df_reg_v2.apply(lambda row: row['Contralateral_Effect'] if row['subregion'] == 'Caudate' else 0, axis=1)
  717. df_reg_v2['Contralateral_Effect_Putamen'] = df_reg_v2.apply(lambda row: row['Contralateral_Effect'] if row['subregion'] == 'Putamen' else 0, axis=1)
  718. df_reg_v2['Contralateral_Effect_Outside'] = df_reg_v2.apply(lambda row: row['Contralateral_Effect'] if row['subregion'] == 'Outside' else 0, axis=1)
  720. coords_v2 = {
  721. "subject": df_reg_v2['Subject'].unique(),
  722. "obs_id": np.arange(len(df_reg_v2))
  723. }
  724. subject_map_v2 = {name: i for i, name in enumerate(df_reg_v2['Subject'].unique())}
  725. subject_indices_for_model_v2 = df_reg_v2['Subject'].map(subject_map_v2).values
  727. with pm.Model(coords=coords_v2) as model_v2:
  728. delta_z_obs = pm.Data("delta_z_obs", df_reg_v2["Δz"].values, dims="obs_id")
  729. fus_side_data = pm.Data("fus_side_data", df_reg_v2["fus_side"].values.astype(float), dims="obs_id")
  730. roi_side_data = pm.Data("roi_side_data", df_reg_v2["roi_side"].values.astype(float), dims="obs_id")
  731. subject_idx_data = pm.Data("subject_idx_data", subject_indices_for_model_v2, dims="obs_id")
  733. int_caudate_data = pm.Data("int_caudate_data", df_reg_v2["Contralateral_Effect_Caudate"].values.astype(float), dims="obs_id")
  734. int_putamen_data = pm.Data("int_putamen_data", df_reg_v2["Contralateral_Effect_Putamen"].values.astype(float), dims="obs_id")
  735. int_outside_data = pm.Data("int_outside_data", df_reg_v2["Contralateral_Effect_Outside"].values.astype(float), dims="obs_id")
  737. mu_alpha_subj = pm.Normal('mu_alpha_subj', mu=0, sigma=2)
  738. sigma_alpha_subj = pm.HalfNormal('sigma_alpha_subj', sigma=2)
  739. alpha_subj_offset = pm.Normal('alpha_subj_offset', mu=0, sigma=1, dims="subject")
  740. alpha_subj = pm.Deterministic('alpha_subj', mu_alpha_subj + alpha_subj_offset * sigma_alpha_subj, dims="subject")
  742. β_fus = pm.Normal('LIFU_Target_Side', mu=0, sigma=2)
  743. β_read = pm.Normal('ROI_Side', mu=0, sigma=2)
  744. β_int_caudate = pm.Normal('Contralateral_Effect_in_Caudate', mu=0, sigma=2)
  745. β_int_putamen = pm.Normal('Contralateral_Effect_in_Putamen', mu=0, sigma=2)
  746. β_int_outside = pm.Normal('Contralateral_Effect_outside_CP', mu=0, sigma=2)
  748. σ_resid = pm.HalfNormal('sigma_resid', sigma=2)
  749. ν_resid = pm.Exponential('nu_resid', 1/30.0)
  751. mu_eq = (alpha_subj[subject_idx_data] +
  752. β_fus * fus_side_data +
  753. β_read * roi_side_data +
  754. β_int_caudate * int_caudate_data +
  755. β_int_putamen * int_putamen_data +
  756. β_int_outside * int_outside_data)
  758. y_obs = pm.StudentT('y', nu=ν_resid, mu=mu_eq, sigma=σ_resid, observed=delta_z_obs, dims="obs_id")
  760. print("Sampling Bayesian sub-region model...")
  761. try:
  762. idata_v2 = pm.sample(5000, tune=2000, chains=16, cores=16, return_inferencedata=True,
  763. idata_kwargs={"log_likelihood": True}, target_accept=0.95, random_seed=81)
  764. ppc_v2 = pm.sample_posterior_predictive(idata_v2, var_names=['y', 'alpha_subj'], random_seed=81, return_inferencedata=True)
  765. idata_v2.extend(ppc_v2)
  766. print("Sub-region model sampling complete.")
  767. except Exception as e:
  768. print(f"Error during PyMC sampling for sub-region model: {e}")
  769. idata_v2 = None
  770. else:
  771. print("Skipping Sub-region Bayesian Regression: No data available.")
  773. # -----------------------------------------------------------------------------
  774. # 13) TRACE & POSTERIOR PAIR PLOTS
  775. # -----------------------------------------------------------------------------
  776. if idata:
  777. var_names_trace = ['mu_alpha_subj', 'sigma_alpha_subj', 'FUS_Target_Side', 'Read_Side', 'Contralateral_Effect', 'sigma_resid', 'nu_resid']
  778. # Check if all variables are in idata before plotting
  779. valid_vars_for_trace = [var for var in var_names_trace if var in idata.posterior]
  780. if valid_vars_for_trace:
  781. # Reverted to original trace plot call
  782. az.plot_trace(idata, var_names=valid_vars_for_trace, compact=True, figsize=(12, len(valid_vars_for_trace)*2))
  783. plt.tight_layout(); plt.show()
  784. else:
  785. print("No valid variables for trace plot found in idata.posterior.")
  787. if 'alpha_subj' in idata.posterior:
  788. az.plot_forest(idata, var_names=['alpha_subj'], combined=True, hdi_prob=0.95, figsize=(8, max(4, len(idata.posterior.subject)*0.3)))
  789. plt.title("Subject-Specific Intercepts (alpha_subj)")
  790. plt.show()
  791. else:
  792. print("'alpha_subj' not found in idata.posterior for forest plot.")
  793. else:
  794. print("Skipping Trace Plots: No idata from Bayesian model.")
  797. #-----------------------------------------------------------------------------
  798. # 15) POSTERIOR KDES of fixed effect coefficients
  799. # -----------------------------------------------------------------------------
  800. if idata:
  801. fixed_effects_kde = ['FUS_Target_Side','Read_Side','Contralateral_Effect']
  802. valid_vars_for_kde_pair = [var for var in fixed_effects_kde if var in idata.posterior]
  804. if valid_vars_for_kde_pair:
  806. df_post_fixed = idata.posterior[valid_vars_for_kde_pair].to_dataframe().reset_index().melt(
  807. id_vars=['chain','draw'], value_vars=valid_vars_for_kde_pair,
  808. var_name='parameter', value_name='value'
  809. )
  810. if not df_post_fixed.empty:
  811. fig, ax = plt.subplots(figsize=(8,7))
  812. # Style maps: different shades of grey/black, line style, and width
  813. style_map = {
  814. 'FUS_Target_Side': {'hatch': None, 'linestyle': ':', 'linewidth': 1, 'label': 'LIFU Target Side', 'text': None},
  815. 'Read_Side': {'hatch': None, 'linestyle': '--', 'linewidth': 1, 'label': 'ROI Side', 'text': None},
  816. 'Contralateral_Effect': {'hatch': '////', 'linestyle': '-', 'linewidth': 2.5, 'label': 'Contralateral Effect', 'text': 'Contralateral'}
  817. }
  818. for param in valid_vars_for_kde_pair:
  819. data = df_post_fixed.query("parameter==@param")['value']
  820. kde = sns.kdeplot(
  821. data=data,
  822. ax=ax,
  823. color='black',
  824. linestyle=style_map[param]['linestyle'],
  825. linewidth=style_map[param]['linewidth'],
  826. fill=False,
  827. label=style_map[param]['label'],
  828. )
  829. if style_map[param]['text']:
  830. x, y = kde.get_lines()[-1].get_data()
  831. # For Caudate, add a light black face shade
  832. if param == 'Contralateral_Effect':
  833. ax.fill_between(x, 0, y, facecolor='black', alpha=0.2, edgecolor='black', hatch=style_map[param]['hatch'])
  834. else:
  835. ax.fill_between(x, 0, y, facecolor='none', edgecolor='black', hatch=style_map[param]['hatch'], alpha=0.2)
  837. ax.axvline(0, ls='--', color='black', lw=1)
  838. ax.set_xlim(-0.1, 0.3)
  839. ax.set_ylim(0, 25)
  840. sns.despine(ax=ax)
  841. ax.set_xlabel('Coefficient value', fontsize=24)
  842. ax.set_ylabel('Density', fontsize=24)
  843. xticks = ax.get_xticks()
  844. ax.set_xticks(xticks[::2])
  845. ax.set_xticklabels([f"{x:.1f}" for x in xticks[::2]], fontsize=22)
  846. ax.set_yticklabels([f"{y:g}" for y in ax.get_yticks()], fontsize=22)
  847. ax.set_title('')
  848. ax.legend(title='', fontsize=22, shadow=False, frameon=True, framealpha=1, edgecolor='black', borderpad=0.2, borderaxespad=0.2)
  849. plt.tight_layout()
  850. plt.savefig(output_dir / "Regression.svg", format='svg')
  851. plt.show()
  852. plt.close(fig)
  853. else:
  854. print("No data for fixed effects KDE plot after melting.")
  855. else:
  856. print("No valid fixed effect variables for KDE plots found in idata.posterior.")
  857. else:
  858. print("Skipping Posterior KDEs: No idata from Bayesian model.")
  860. # -----------------------------------------------------------------------------
  861. # 15a) POSTERIOR KDES of fixed effect coefficients (SUB-REGIONS)
  862. # -----------------------------------------------------------------------------
  863. if idata_v2:
  864. fixed_effects_kde_v2 = ['LIFU_Target_Side', 'ROI_Side', 'Contralateral_Effect_in_Caudate', 'Contralateral_Effect_in_Putamen', 'Contralateral_Effect_outside_CP']
  865. valid_vars_for_kde_v2 = [var for var in fixed_effects_kde_v2 if var in idata_v2.posterior]
  867. if valid_vars_for_kde_v2:
  869. df_post_fixed_v2 = idata_v2.posterior[valid_vars_for_kde_v2].to_dataframe().reset_index().melt(
  870. id_vars=['chain','draw'], value_vars=valid_vars_for_kde_v2,
  871. var_name='parameter', value_name='value'
  872. )
  873. if not df_post_fixed_v2.empty:
  874. # Improved KDE plot for fixed effect coefficients (Sub-region Model) with distinguishable groups (no color)
  876. # Define hatching and line style for each group
  877. group_hatch_map = {
  878. 'Contralateral_Effect_in_Caudate': {'hatch': '////', 'linestyle': ':', 'linewidth': 2.5, 'label': 'Contralateral Effect in Caudate', 'text': 'Caudate'},
  879. 'Contralateral_Effect_in_Putamen': {'hatch': '\\\\\\', 'linestyle': '-', 'linewidth': 2.5, 'label': 'Contralateral Effect in Putamen', 'text': 'Putamen'},
  880. 'Contralateral_Effect_outside_CP': {'hatch': None, 'linestyle': '-.', 'linewidth': 2.5, 'label': 'Contralateral Effect Outside Striatum', 'text': 'Outside'},
  881. 'LIFU_Target_Side': {'hatch': None, 'linestyle': ':', 'linewidth': 1, 'label': 'LIFU Target Side', 'text': None},
  882. 'ROI_Side': {'hatch': None, 'linestyle': '--', 'linewidth': 1, 'label': 'ROI Side', 'text': None},
  883. }
  885. fig, ax = plt.subplots(figsize=(8, 7))
  887. peak_positions = {}
  889. for param in valid_vars_for_kde_v2:
  890. if param in group_hatch_map:
  891. data = df_post_fixed_v2.query("parameter==@param")['value']
  892. # Plot KDE line
  893. kde = sns.kdeplot(
  894. data=data,
  895. ax=ax,
  896. color='black',
  897. linestyle=group_hatch_map[param]['linestyle'],
  898. linewidth=group_hatch_map[param]['linewidth'],
  899. fill=False,
  900. label=group_hatch_map[param]['label'],
  901. )
  902. # Overlay hatched fill for the three main groups
  903. if group_hatch_map[param]['text']:
  904. x, y = kde.get_lines()[-1].get_data()
  905. # For Caudate, add a light black face shade
  906. if param == 'Contralateral_Effect_in_Putamen':
  907. ax.fill_between(x, 0, y, facecolor='black', alpha=0.2, edgecolor='black', hatch=group_hatch_map[param]['hatch'])
  908. else:
  909. ax.fill_between(x, 0, y, facecolor='none', edgecolor='black', hatch=group_hatch_map[param]['hatch'], alpha=0.2)
  910. # Store peak position for annotation
  911. if 'text' in group_hatch_map[param]:
  912. peak_idx = np.argmax(y)
  913. peak_x = x[peak_idx]
  914. peak_y = y[peak_idx]
  915. peak_positions[param] = (peak_x, peak_y)
  917. # Add text annotations for Caudate, Putamen, Outside
  918. for param, (peak_x, peak_y) in peak_positions.items():
  919. text = group_hatch_map[param]['text']
  920. # Offset 'Outside' label to the right
  921. if text == 'Outside':
  922. peak_x = peak_x + 0.015
  923. peak_y = peak_y - 0.3
  924. if text == 'Putamen':
  925. peak_x = peak_x + 0.01
  926. peak_y = peak_y - 0.2
  927. ax.text(peak_x + 0.01, peak_y + 0.02 * ax.get_ylim()[1], text, ha='center', va='bottom', fontsize=18, fontweight='bold', color='black', bbox=dict(facecolor='white', edgecolor='none', alpha=0.7, pad=1.5))
  929. ax.axvline(0, ls='--', color='black', lw=1)
  930. ax.set_xlim(-0.1, 0.3)
  931. ax.set_ylim(0, 25)
  932. sns.despine(ax=ax)
  933. ax.set_xlabel('Coefficient value', fontsize=24)
  934. ax.set_ylabel('Density', fontsize=24)
  935. xticks = ax.get_xticks()
  936. # Use 0.1 precision and skip every other tick
  937. ax.set_xticks(xticks[::2])
  938. ax.set_xticklabels([f"{x:.1f}" for x in xticks[::2]], fontsize=22)
  939. ax.set_yticklabels([f"{y:g}" for y in ax.get_yticks()], fontsize=22)
  940. ax.set_title('')
  941. ax.legend(title='', fontsize=17.5, shadow=False, frameon=True, framealpha=1, edgecolor='black', borderpad=0.2, borderaxespad=0.2)
  942. plt.tight_layout()
  943. plt.savefig(output_dir / "Regression_subregion.svg", format='svg')
  944. plt.show()
  945. plt.close(fig)
  946. else:
  947. print("No data for fixed effects KDE plot (sub-region model) after melting.")
  948. else:
  949. print("No valid fixed effect variables for KDE plots found in idata_v2.posterior.")
  950. else:
  951. print("Skipping Posterior KDEs for Sub-region Model: No idata_v2 from Bayesian model.")

3-Combined_Monkeys_OEF_Bayesian_template.ipynb, under apgl-v3 · at the source

Overview

Authors: Soroosh Sanatkhani1, Dong Liu1, Fabian Munoz1, Jack Grinband2,3, Elisa E Konofagou3,4, Vincent P Ferrera1,2,5
ORCID iDs: Soroosh Sanatkhani, Jack Grinband, Elisa E Konofagou, Vincent P Ferrera
  1. Zuckerman Mind Brain Behavior Institute, Columbia University, New York, NY, USA
  2. Department of Psychiatry, Columbia University, New York, NY, USA
  3. Department of Radiology, Columbia University, New York, NY, USA
  4. Department of Biomedical Engineering, Columbia University, New York, NY, USA
  5. Department of Neuroscience, Columbia University, New York, NY, USA
Institutions: Mortimer B. Zuckerman Mind Brain Behavior Institute (United States); Columbia University (United States)
Journal: Science advances, volume 12, issue 29, article eaed4944
Dates: received 29 October 2025; accepted 2 June 2026; published online 17 July 2026; in print July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1126/sciadv.aed4944 · PMID 42467769 · PMCID PMC13378551 · OpenAlex W7169524987
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: structural MRI / diffusion (modality), other (modality), non-human primate (organism), clinical / translational (subfield)
Methods: Statistics, fMRI & imaging
MeSH: Blood-Brain Barrier*, Brain*, Ultrasonic Waves*, Animals, Macaca mulatta, Magnetic Resonance Imaging, Male, Oxygen, Ultrasonography (* major topic)
Topic: Ultrasound and Hyperthermia Applications (Biomedical Engineering, Engineering), according to OpenAlex
Funding: NIMH NIH HHS (R01 MH133020)
Citations: not cited yet (Europe PMC); 66 references in the paper

Abstract

Low-intensity focused ultrasound (LIFU) is a promising technique for opening the blood-brain barrier (BBB) for drug delivery, but its physiological consequences in remote brain regions remain a major blind spot for clinical safety and efficacy. To address this gap, we performed the quantitative mapping of brain metabolism following a focal LIFU-induced BBB opening in a nonhuman primate model, using quantitative BOLD MRI to measure the oxygen extraction fraction (OEF). We report a paradoxical response: While the targeted striatum showed no significant metabolic changes, we observed a profound and spatially specific increase in OEF in the homologous contralateral striatum, an effect predominantly driven by the putamen. These findings demonstrate that focal BBB opening is not merely a localized vascular event but a potent neuromodulatory intervention that induces metabolic stress in distant, untreated brain regions, a discovery with critical implications for the safe and effective clinical translation of all focal brain therapies.

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

Repository

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

Zenodo 19188809

License: apgl-v3
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Languages: Jupyter (4)
Size: 4 files, 4 scripts
Software Heritage: not checked
Found in: “Data, code, and materials availability:”
Holds: 4 notebooks
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Matplotlib (4 files), NumPy (4 files), NiBabel (3 files), pandas (3 files), seaborn (3 files), ArviZ (2 files), PyMC (2 files), SciPy (2 files), statsmodels (2 files), AFNI (1 file), CuPy (1 file), JAX (1 file), Numba (1 file), PyTorch (1 file), scikit-image (1 file), TensorFlow (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
4 files
At the source: doi.org/10.5281/zenodo.19188809

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:

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

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

Data

Datasets cited

Data, code, and materials availability

All data and code needed to evaluate and reproduce the results in the paper are present in the paper and/or the Supplementary Materials. The raw MRI datasets supporting the qBOLD analysis have been deposited in a permanent, independent repository and can be accessed at https://doi.org/10.5281/zenodo.19186774. The computational pipeline, consisting of the Jupyter notebooks required to generate the OEF maps and reproduce the statistical analyses, is archived and is available at https://doi.org/10.5281/zenodo.19188809.

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

Versions

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

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 6 authors, 9 MeSH terms, 1 funder, 63 references.

Cite

This paper

Sanatkhani, S., Liu, D., Munoz, F., Grinband, J., Konofagou, E. E., & Ferrera, V. P. (2026). Focused ultrasound blood-brain barrier opening reveals a paradoxical remote metabolic response in the primate brain. Science advances, 12(29), eaed4944. https://doi.org/10.1126/sciadv.aed4944

BibTeX

@article{sanatkhani2026focused,
author = {Sanatkhani, Soroosh and Liu, Dong and Munoz, Fabian and Grinband, Jack and Konofagou, Elisa E and Ferrera, Vincent P},
title = {{Focused ultrasound blood-brain barrier opening reveals a paradoxical remote metabolic response in the primate brain}},
journal = {Science advances},
year = {2026},
month = jul,
volume = {12},
number = {29},
pages = {eaed4944},
publisher = {American Association for the Advancement of Science},
issn = {2375-2548},
doi = {10.1126/sciadv.aed4944},
url = {https://doi.org/10.1126/sciadv.aed4944},
pmid = {42467769},
pmcid = {PMC13378551}
}

RIS

TY - JOUR
AU - Sanatkhani, Soroosh
AU - Liu, Dong
AU - Munoz, Fabian
AU - Grinband, Jack
AU - Konofagou, Elisa E
AU - Ferrera, Vincent P
TI - Focused ultrasound blood-brain barrier opening reveals a paradoxical remote metabolic response in the primate brain
T2 - Science advances
J2 - Sci Adv
PY - 2026
DA - 2026/07/17
VL - 12
IS - 29
SP - eaed4944
SN - 2375-2548
PB - American Association for the Advancement of Science
DO - 10.1126/sciadv.aed4944
UR - https://doi.org/10.1126/sciadv.aed4944
LA - en
ER -

CSL-JSON

{
"id": "10.1126/sciadv.aed4944",
"type": "article-journal",
"title": "Focused ultrasound blood-brain barrier opening reveals a paradoxical remote metabolic response in the primate brain",
"container-title": "Science advances",
"author": [
{
"family": "Sanatkhani",
"given": "Soroosh"
},
{
"family": "Liu",
"given": "Dong"
},
{
"family": "Munoz",
"given": "Fabian"
},
{
"family": "Grinband",
"given": "Jack"
},
{
"family": "Konofagou",
"given": "Elisa E"
},
{
"family": "Ferrera",
"given": "Vincent P"
}
],
"container-title-short": "Sci Adv",
"volume": "12",
"issue": "29",
"page": "eaed4944",
"DOI": "10.1126/sciadv.aed4944",
"PMID": "42467769",
"PMCID": "PMC13378551",
"ISSN": "2375-2548",
"publisher": "American Association for the Advancement of Science",
"URL": "https://doi.org/10.1126/sciadv.aed4944",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
17
]
]
}
}

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

Similar papers

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

[1] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: ArviZ, PyMC, CuPy, 11 other tools
[2] doi:10.1038/s41467-026-72940-5 [code]
Cerebellar growth is associated with domain-specific cerebral maturation and socio-linguistic behavior.
Journal: Nature communications
In common: ArviZ, PyMC, TensorFlow, 7 other tools, 1 reference
[3] doi:10.1162/imag.a.1198 [code]
MEPrep: A robust pipeline for multi-echo fMRI denoising and preprocessing.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: AFNI, scikit-image, NiBabel, 5 other tools, author Jack Grinband
[4] doi:10.1371/journal.pcbi.1014340 [code]
pyhgf: A neural network library for predictive coding.
Journal: PLoS computational biology
In common: ArviZ, PyMC, JAX, 5 other tools, 1 reference
[5] doi:10.1371/journal.pone.0346575 [code]
Statistically valid explainable black-box machine learning: applications in sex classification across species using brain imaging.
Journal: PloS one
In common: JAX, Numba, TensorFlow, 7 other tools, non-human primate, structural MRI / diffusion
[6] doi:10.1038/s41592-026-03057-2 [code]
CREsted: modeling genomic and synthetic cell-type-specific enhancers across tissues and species.
Journal: Nature methods
In common: CuPy, Numba, TensorFlow, 7 other tools
[7] doi:10.1162/imag.a.1269 [code]
From early to contemporary normative modeling: Mapping individual differences in neurophysiological signals.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: ArviZ, PyMC, statsmodels, 6 other tools
[8] doi:10.1016/j.isci.2026.116825 [code]
Social hierarchy shapes behavioral and transcriptional responses to chronic stress and ketamine in male mice.
Journal: iScience
In common: CuPy, Numba, scikit-image, 7 other tools
[9] doi:10.1038/s41591-026-04287-9 [code]
An international mega-analysis of psychedelic drug effects on brain circuit function.
Journal: Nature medicine
In common: ArviZ, PyMC, NiBabel, 5 other tools, 1 reference
[10] doi: [code]
Naturalistic behavior and self-generated neural activity predictive of self-correction
Journal: bioRxiv : the preprint server for biology
In common: CuPy, JAX, scikit-image, 6 other tools

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.