OSCR

Restoring access to long-term social recognition memories disrupted by sleep deprivation.

Code ↔ Paper

6 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 6 matches
  1. [1] § MATERIALS AND METHODS › Statistical analysis ↔ src/stats.py, lines 101–150 · score 0.69 · EGFP positive cells, dual engram, cFos, engram reactivation, mCherry, overlap
  2. [2] § MATERIALS AND METHODS › Statistical analysis ↔ src/plots.py, lines 99–128 · score 0.68 · EGFP positive cells, dual engram, cFos, engram reactivation, mCherry, overlap
  3. [3] § RESULTS › Optogenetic stimulation of DG engram cells enables persistent retrieval of social memory several days later by increasing engram reactivation ↔ run_analysis.py, lines 118–178 · score 0.60 · positive cells, engram reactivation, cFos, mCherry, blade, colocalization
  4. [4] § RESULTS › Optogenetic stimulation of DG engram cells enables persistent retrieval of social memory several days later by increasing engram reactivation ↔ src/data_loader.py, lines 68–88 · score 0.60 · positive cells, engram reactivation, cFos, mCherry, blade, colocalization
  5. [5] § MATERIALS AND METHODS › Cell counting ↔ run_analysis.py, lines 118–178 · score 0.60 · EGFP positive cells, cFos, mCherry, blades, colocalization, model
  6. [6] § MATERIALS AND METHODS › Cell counting ↔ src/data_loader.py, lines 90–114 · score 0.59 · EGFP positive cells, cFos, mCherry, blades, colocalization, animal

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

Python · 515 lines · 25 KB · MIT · 2 matches

  1. """
  2. run_analysis.py — V5.1
  3. =======================
  4. Pipeline reproducible V5.1. Genera exactamente los outputs del manuscrito.
  6. Correcciones V5.1 respecto a V5:
  7. - bootstrap_rmse_difference: predicciones fijas, índices remuestreados (N=1000)
  8. - power_formation_contrast: parámetros por marcador (EGFP/mCherry), N_sim=5000
  9. - Model averaging: agregado (acc/mix/sto fixR) con varianza incondicional
  10. - AICc individual corregido: k = k_mech + 8, n = 275
  11. - Ceros declarados; sensibilidad pseudocuenta
  12. - version = V5.1 en CONFIG y run_metadata
  14. Para publicación: N_STARTS=100, N_BOOT=1000, N_SIM_REC=500
  15. """
  17. import sys, json, time, shutil, hashlib, datetime
  18. import numpy as np, pandas as pd
  19. from pathlib import Path
  21. sys.path.insert(0, str(Path(__file__).parent))
  23. from src.config import TABLE_DIR, MODEL_DIR, FIG_DIR, OUTPUT_DIR, DATA_DIR
  24. from src import data_loader, stats, fitting, validation, plots
  25. from src.model import (identifiability_proposition, cellular_analytic_estimate,
  26. omega_lambda_compatibility_map, latent_states_summary,
  27. parameter_spec, test_equivalence, model_averaged_estimates,
  28. power_formation_contrast, power_decay_detection,
  29. cellular_individual_nll)
  31. # ══════════════════════════════════════════════════════════════════════════════
  32. # CONFIGURACIÓN
  33. # ══════════════════════════════════════════════════════════════════════════════
  34. CONFIG = {
  35. "N_STARTS": 50,
  36. "N_BOOT": 500,
  37. "N_SIM_REC": 100,
  38. "CELL_WEIGHT": 1.0,
  39. "SEED": 12345,
  40. "version": "V5.1",
  41. }
  42. # ══════════════════════════════════════════════════════════════════════════════
  45. def _s(obj):
  46. if isinstance(obj, (bool, np.bool_)): return bool(obj)
  47. if isinstance(obj, np.integer): return int(obj)
  48. if isinstance(obj, (float, np.floating)): return float(obj)
  49. if isinstance(obj, np.ndarray): return obj.tolist()
  50. return obj
  53. def _verify_sync(config):
  54. meta = json.loads((MODEL_DIR / "run_metadata.json").read_text())
  55. for key in ("N_STARTS", "N_BOOT", "N_SIM_REC", "SEED", "version"):
  56. assert str(meta.get(key)) == str(config.get(key)), \
  57. f"DESINCRONIZACIÓN: {key} metadata={meta.get(key)} config={config.get(key)}"
  58. print(" ✓ Sincronización verificada")
  61. def _aicc_individual_corrected(fits_ind, n_obs_correct=275):
  62. """
  63. AICc individual con k = k_mech + 8 (cuatro intercepts + cuatro varianzas)
  64. y n = 275 (189 beh + 44 form + 42 reac).
  65. """
  66. rows = []
  67. for key, fit in fits_ind.items():
  68. k_mech = len(fit["names"])
  69. k_corr = k_mech + 8
  70. ll = -fit["result"].fun
  71. aic = 2*k_corr - 2*ll
  72. aicc = aic + (2*k_corr*(k_corr+1)) / max(n_obs_correct - k_corr - 1, 1)
  73. bic = k_corr * np.log(n_obs_correct) - 2*ll
  74. rows.append({"key": key, "k_mech": k_mech, "k_corr": k_corr,
  75. "logLik": ll, "AICc_corrected": aicc, "BIC_corrected": bic})
  76. df = pd.DataFrame(rows).sort_values("AICc_corrected").reset_index(drop=True)
  77. d = df["AICc_corrected"] - df["AICc_corrected"].min()
  78. w = np.exp(-0.5 * d)
  79. df["delta_AICc"] = d; df["AICc_weight"] = w / w.sum()
  80. return df
  83. def _model_averaging_unconditional(fits_fixR, ic_fixR, boot_raw_df):
  84. """
  85. Model averaging con varianza incondicional = within + between.
  86. Within-model variance: bootstrap de access_fixR (solo el disponible).
  87. Between-model variance: dispersión entre estimaciones ponderadas.
  88. """
  89. rows = []
  90. for key, fit in fits_fixR.items():
  91. w = ic_fixR.loc[ic_fixR["key"] == key, "AICc_weight"].values
  92. if len(w) == 0: continue
  93. rows.append({"key": key, "weight": float(w[0]),
  94. "d_E": fit["parameters"].get("d_E", 0.0),
  95. "d_C": fit["parameters"].get("d_C", 0.0)})
  96. df_w = pd.DataFrame(rows)
  97. d_C_MA = float((df_w["weight"] * df_w["d_C"]).sum())
  98. d_E_MA = float((df_w["weight"] * df_w["d_E"]).sum())
  99. # within-model (bootstrap access_fixR)
  100. var_within_C = float(boot_raw_df["d_C"].var(ddof=1)) if "d_C" in boot_raw_df else 0.0
  101. # between-model
  102. var_between_C = float((df_w["weight"] * (df_w["d_C"] - d_C_MA)**2).sum())
  103. var_total_C = var_within_C + var_between_C
  104. sd_total_C = var_total_C**0.5
  105. ci_lo = d_C_MA - 1.96 * sd_total_C
  106. ci_hi = d_C_MA + 1.96 * sd_total_C
  107. return {
  108. "d_E_MA": d_E_MA, "d_C_MA": d_C_MA,
  109. "var_within_d_C": var_within_C, "var_between_d_C": var_between_C,
  110. "var_total_d_C": var_total_C, "sd_total_d_C": sd_total_C,
  111. "ci_lo_unconditional": float(ci_lo), "ci_hi_unconditional": float(ci_hi),
  112. "note": ("Aggregate-model MA: access/mixed/storage-fixR weighted by AICc. "
  113. "Unconditional CI = within (bootstrap access_fixR) + between models."),
  114. "details": df_w.to_dict(orient="records"),
  115. }
  118. def _pseudocount_sensitivity(data, best_fit_ind):
  119. """Sensibilidad de d_C a ceros celulares con pseudocuenta log(v+0.5)."""
  120. from scipy.optimize import minimize
  121. cell_raw = data["cell_raw"]
  122. p12 = cell_raw[(cell_raw["protocol_id"] == "P12_dual_engram_tagging")
  123. & (cell_raw["blade"] == "whole_DG")]
  124. reac_specs = [("engram_reactivation_colocalization", "EGFP+cFos", 0),
  125. ("engram_reactivation_colocalization", "mCherry+cFos", 1)]
  126. form_specs = [("EGFP_positive_cells", "EGFP", 0),
  127. ("mCherry_positive_cells", "mCherry", 1)]
  129. log_yF = []; sdF = []; mF = []
  130. for meas, marker, midx in form_specs:
  131. d = p12[(p12["measurement"] == meas) & (p12["marker"] == marker)]
  132. for _, row in d.iterrows():
  133. v = row["value"]
  134. if v <= 0: continue
  135. log_yF.append(np.log(v)); sdF.append(1 if row["sleep_status"] == "SD" else 0)
  136. mF.append(midx)
  138. log_yR_p = []; sdR_p = []; mR_p = []
  139. for meas, marker, midx in reac_specs:
  140. d = p12[(p12["measurement"] == meas) & (p12["marker"] == marker)]
  141. for _, row in d.iterrows():
  142. log_yR_p.append(np.log(row["value"] + 0.5)) # pseudocuenta
  143. sdR_p.append(1 if row["sleep_status"] == "SD" else 0)
  144. mR_p.append(midx)
  146. cell_pseudo = {
  147. "log_y_F_arr": np.array(log_yF), "sd_flag_F_arr": np.array(sdF), "m_idx_F_arr": np.array(mF),
  148. "log_y_R_arr": np.array(log_yR_p), "sd_flag_R_arr": np.array(sdR_p), "m_idx_R_arr": np.array(mR_p),
  149. "n_F": len(log_yF), "n_R": len(log_yR_p),
  150. }
  152. model = best_fit_ind["model"]; decay = best_fit_ind["decay"]; fix_r = best_fit_ind["fix_rescues"]
  153. names, lower, upper = parameter_spec(model, decay, fix_r)
  154. bounds = list(zip(lower, upper))
  155. rng = np.random.default_rng(42)
  157. from src.model import cellular_individual_nll_theta
  158. best = None
  159. for _ in range(20):
  160. x0 = rng.uniform(lower, upper)
  161. res = minimize(cellular_individual_nll_theta, x0=x0,
  162. args=(data["cal_arr"], model, cell_pseudo, decay, fix_r),
  163. method="L-BFGS-B", bounds=bounds,
  164. options={"maxiter": 4000, "ftol": 1e-12})
  165. if best is None or res.fun < best.fun: best = res
  167. params_p = dict(zip(names, best.x))
  168. return {
  169. "method": "pseudocount_log_v_plus_0.5",
  170. "n_R_with_pseudo": len(log_yR_p),
  171. "n_zeros_included": 2,
  172. "d_C_pseudo": float(params_p.get("d_C", float("nan"))),
  173. "d_C_original": float(best_fit_ind["parameters"].get("d_C", float("nan"))),
  174. "difference": float(abs(params_p.get("d_C", 0) - best_fit_ind["parameters"].get("d_C", 0))),
  175. "note": ("Two reactivation values were exactly zero in the original data. "
  176. "The primary analysis excludes them (log-transformation). "
  177. "This sensitivity uses log(v+0.5)."),
  178. }
  181. def main():
  182. t0 = time.time()
  183. print("=" * 65)
  184. print(" Sleep Memory Model V5.1 — Pipeline reproducible")
  185. print("=" * 65)
  187. print("\n[0] Limpiando output/...")
  188. if OUTPUT_DIR.exists(): shutil.rmtree(OUTPUT_DIR)
  189. for d in (FIG_DIR, TABLE_DIR, MODEL_DIR): d.mkdir(parents=True, exist_ok=True)
  191. # [1] Metadatos
  192. print("\n[1] Metadatos...")
  193. data_hash = fitting.compute_data_hash(
  194. DATA_DIR / "behavior_paired.csv", DATA_DIR / "cellular_long.csv")
  195. code_hash = fitting.compute_code_hash()
  196. fitting.save_run_metadata(CONFIG, data_hash, code_hash)
  197. print(f" data={data_hash} code={code_hash}")
  199. # [2] Datos
  200. print("\n[2] Datos...")
  201. data = data_loader.load_all()
  202. audit = data["spi_audit"]
  203. print(f" SPI audit: {'✓' if audit['passed'] else '✗'} "
  204. f"Cal={len(data['cal'])} Val={len(data['val'])}")
  205. print(f" cell_ind: n_F={data['cell_ind']['n_F']} n_R={data['cell_ind']['n_R']} "
  206. f"(2 zeros excluidos de reactivación)")
  208. # [3] Equivalencia
  209. print("\n[3] Test de equivalencia...")
  210. eq = test_equivalence(data["cal_arr"])
  211. assert eq["all_passed"], "EQUIVALENCIA FALLIDA"
  212. print(" ✓ storage ≡ access en conducta (max_diff = 0)")
  214. # [4] Estadísticas descriptivas
  215. print("\n[4] Estadística descriptiva...")
  216. stats.behavior_summary(data["memory"]).to_csv(TABLE_DIR / "table1_behavior_summary.csv", index=False)
  217. stats.one_sample_spi_tests(data["memory"]).to_csv(TABLE_DIR / "table1b_spi_vs_zero.csv", index=False)
  218. cell_cont = stats.cellular_contrasts(data["cell_raw"])
  219. cell_cont.to_csv(TABLE_DIR / "table2_cellular_contrasts.csv", index=False)
  220. stats.exploration_confound_tests(data["memory"]).to_csv(TABLE_DIR / "tableS1_exploration_confound.csv", index=False)
  222. cell_est = cellular_analytic_estimate(data["cell_obs"])
  223. print(f" {cell_est['message']}")
  224. (TABLE_DIR / "cell_analytic_estimate.json").write_text(
  225. json.dumps({k: _s(v) for k, v in cell_est.items()}, indent=2))
  226. prop = identifiability_proposition(data["cell_obs"])
  227. (TABLE_DIR / "identifiability_proposition.json").write_text(
  228. json.dumps({k: _s(v) for k, v in prop.items()
  229. if k not in ("d_E_equipotential", "d_C_equipotential")}, indent=2))
  231. # [5] Ajuste familia simétrica (agregado — sensibilidad/secundario)
  232. print(f"\n[5] Familia simétrica (N_STARTS={CONFIG['N_STARTS']})...")
  233. fits, ic = fitting.run_fitting(data["cal_arr"], data["cell_obs"],
  234. n_starts=CONFIG["N_STARTS"],
  235. cell_weight=CONFIG["CELL_WEIGHT"],
  236. seed=CONFIG["SEED"], verbose=True)
  237. best_key = fitting.best_fit_key(ic)
  238. best_fit = fits[best_key]
  239. fitting.parameter_table(best_fit).to_csv(TABLE_DIR / "table3_parameters.csv", index=False)
  240. ic.to_csv(TABLE_DIR / "table3_model_comparison.csv", index=False)
  241. print(f"\n Mejor agregado: {best_key}")
  243. # Solo conducta
  244. validation.sensitivity_behavior_only(data["cal_arr"], data["cell_obs"],
  245. n_starts=CONFIG["N_STARTS"] // 2,
  246. seed=CONFIG["SEED"] + 10)["ic"].to_csv(
  247. TABLE_DIR / "tableS3_sensitivity_beh_only.csv", index=False)
  249. # [5b] Modelo celular individual (PRIMARIO)
  250. print("\n[5b] Modelo celular individual (PRIMARIO)...")
  251. fits_ind = {}
  252. for model in ("storage", "access", "mixed"):
  253. fits_ind[f"{model}_IND"] = fitting.fit_model_individual(
  254. model, data["cal_arr"], data["cell_ind"],
  255. n_starts=CONFIG["N_STARTS"], decay=False, fix_rescues=True,
  256. seed=CONFIG["SEED"], verbose=True)
  258. # AICc corregido: k = k_mech + 8, n = 275
  259. ic_ind_corr = _aicc_individual_corrected(fits_ind, n_obs_correct=275)
  260. ic_ind_corr.to_csv(TABLE_DIR / "table3b_model_comparison_individual_corrected.csv", index=False)
  261. print("\n AICc individual CORREGIDO (k+8, n=275):")
  262. print(ic_ind_corr[["key", "k_corr", "AICc_corrected", "delta_AICc", "AICc_weight"]].to_string(index=False))
  264. best_fit_ind = fits_ind["access_IND"]
  265. print(f"\n d_C individual MLE: {best_fit_ind['parameters'].get('d_C', 0):.3f}")
  267. # [5c] Sensibilidades
  268. print("\n[5c] Sensibilidades...")
  269. sw = fitting.sensitivity_cell_weight(data["cal_arr"], data["cell_obs"],
  270. n_starts=CONFIG["N_STARTS"] // 2, seed=77)
  271. sw.to_csv(TABLE_DIR / "tableS_sensitivity_cell_weight.csv", index=False)
  272. se_sens = fitting.sensitivity_cell_error(data["cal_arr"], data["cell_obs"],
  273. n_starts=CONFIG["N_STARTS"] // 2, seed=88)
  274. se_sens.to_csv(TABLE_DIR / "tableS_sensitivity_cell_error.csv", index=False)
  275. fitting.baselines_beh_aic(data["cal_arr"], len(data["cal_arr"]["y"])).to_csv(
  276. TABLE_DIR / "tableS_baselines_beh_aic.csv", index=False)
  278. # Pseudocuenta sensibilidad
  279. print(" Pseudocuenta sensibilidad...")
  280. zero_sens = _pseudocount_sensitivity(data, best_fit_ind)
  281. (TABLE_DIR / "tableS_zero_sensitivity.json").write_text(json.dumps(zero_sens, indent=2))
  282. print(f" d_C original={zero_sens['d_C_original']:.3f} "
  283. f"d_C_pseudo={zero_sens['d_C_pseudo']:.3f} diff={zero_sens['difference']:.3f}")
  285. # [5d] Potencia con parámetros reales por marcador, N_sim=5000
  286. print("\n[5d] Potencia (N_sim=5000, parámetros por marcador)...")
  287. power_rows = []
  288. for d_E in [0.082, 0.10, 0.15, 0.20, 0.30, 0.40]:
  289. for n in [8, 10, 12, 20]:
  290. for marker in ["EGFP", "mCherry"]:
  291. p = power_formation_contrast(n, d_E, marker=marker,
  292. n_sim=5000, seed=42)
  293. power_rows.append({"d_E_true": d_E, "n_per_group": n,
  294. "marker": marker, "power": round(p, 3)})
  295. power_df = pd.DataFrame(power_rows)
  296. power_df.to_csv(TABLE_DIR / "tableS_power_formation.csv", index=False)
  297. print(" Potencia n=10:")
  298. print(power_df[power_df["n_per_group"] == 10].pivot(
  299. index="d_E_true", columns="marker", values="power").to_string())
  301. pd_rows = []
  302. for muq in [0.05, 0.10, 0.20, 0.30]:
  303. p = power_decay_detection(muq, n_per_group=10, n_sim=2000, seed=43)
  304. pd_rows.append({"mu_q_true": muq, "power": round(p, 3)})
  305. pd.DataFrame(pd_rows).to_csv(TABLE_DIR / "tableS_power_decay.csv", index=False)
  307. # [6] Calibración y validación
  308. print("\n[6] Calibración y validación...")
  309. cal_out, cal_diag, cal_group = validation.calibration_diagnostics(
  310. data["cal"], data["cal_arr"], best_fit)
  311. val_out, val_metrics, group_metrics = validation.external_validation(
  312. data["val"], data["val_arr"], best_fit)
  313. print(f" Cal RMSE={cal_diag['rmse']:.4f} "
  314. f"Val RMSE={val_metrics['rmse']:.4f} R²={val_metrics['r2']:.4f}")
  316. lopo_df = validation.leave_one_protocol_out(
  317. data["cal"], data["cell_obs"], best_fit,
  318. n_starts=20, cell_weight=CONFIG["CELL_WEIGHT"], seed=CONFIG["SEED"] + 1)
  319. lopo_df.to_csv(TABLE_DIR / "table5_lopo.csv", index=False)
  320. lopo_valid = lopo_df[~lopo_df["is_stress_test"]]
  321. lopo_stress = lopo_df[lopo_df["is_stress_test"]]
  322. print(f" LOPO válido={lopo_valid['rmse'].mean():.4f} "
  323. f"stress={lopo_stress['rmse'].mean():.4f}")
  325. fitting.compare_mechanistic_vs_baselines_rmse(
  326. fits, data["cal_arr"], data["val_arr"], data["val"]).to_csv(
  327. TABLE_DIR / "table_rmse_vs_baselines.csv", index=False)
  328. validation.save_validation_results(
  329. val_out, val_metrics, group_metrics, pd.DataFrame(), cal_diag, cal_group)
  330. bci = TABLE_DIR / "bootstrap_ci.csv"
  331. if bci.exists() and bci.stat().st_size < 50: bci.unlink()
  333. # Decay comparison
  334. fits_d = {}
  335. for model in ("access", "mixed"):
  336. for dec in (True, False):
  337. key = f"{model}_fixR_{'decay' if dec else 'nodecay'}"
  338. fits_d[key] = fitting.fit_model(
  339. model, data["cal_arr"], data["cell_obs"],
  340. n_starts=CONFIG["N_STARTS"] // 2, cell_weight=1.0,
  341. decay=dec, fix_rescues=True, seed=CONFIG["SEED"] + 99, verbose=True)
  342. fitting.information_criteria(
  343. fits_d, len(data["cal_arr"]["y"]) + len(data["cell_obs"])).to_csv(
  344. TABLE_DIR / "tableS_decay_comparison.csv", index=False)
  346. # [7] Bootstrap
  347. print(f"\n[7] Bootstrap (n={CONFIG['N_BOOT']})...")
  348. boot_df = validation.run_bootstrap(
  349. cal_df=data["cal"], cell_df=data["cell_raw"], best_fit=best_fit,
  350. n_boot=CONFIG["N_BOOT"], n_starts_boot=3,
  351. cell_weight=CONFIG["CELL_WEIGHT"], seed=CONFIG["SEED"] + 2, verbose=True)
  352. ci_df = validation.bootstrap_ci(boot_df)
  353. identif = validation.parameter_identifiability(boot_df)
  354. boot_df.to_csv(TABLE_DIR / "bootstrap_raw.csv", index=False)
  355. ci_df.to_csv(TABLE_DIR / "table4_bootstrap_ci.csv", index=False)
  356. identif.to_csv(TABLE_DIR / "tableS2_identifiability_corr.csv", index=False)
  357. print(ci_df[["parameter", "mean", "ci_low", "ci_high"]].to_string(index=False))
  359. print(" Bootstrap individual (n=200)...")
  360. boot_ind = validation.run_bootstrap_individual(
  361. cal_df=data["cal"], cell_ind=data["cell_ind"],
  362. best_fit=best_fit_ind, n_boot=200, n_starts_boot=3,
  363. seed=CONFIG["SEED"] + 3, verbose=True)
  364. ci_ind = validation.bootstrap_ci(boot_ind)
  365. boot_ind.to_csv(TABLE_DIR / "bootstrap_individual_raw.csv", index=False)
  366. ci_ind.to_csv(TABLE_DIR / "table4b_bootstrap_ci_individual.csv", index=False)
  368. # Model averaging con varianza incondicional (usando bootstrap ya calculado)
  369. # Usa modelos agregados (access/mixed/storage fixR) — declarado como tal
  370. fits_fixR = {k: v for k, v in fits.items() if k in ("access_fixR", "mixed_fixR", "storage_fixR")}
  371. ic_fixR = ic[ic["key"].isin(fits_fixR.keys())].copy()
  372. ma = _model_averaging_unconditional(fits_fixR, ic_fixR, boot_df)
  373. (TABLE_DIR / "model_averaging.json").write_text(
  374. json.dumps({k: _s(v) for k, v in ma.items() if k != "details"}, indent=2))
  375. print(f" MA (agregado): d_C={ma['d_C_MA']:.3f} "
  376. f"CI_uncond=[{ma['ci_lo_unconditional']:.3f},{ma['ci_hi_unconditional']:.3f}]")
  378. # [8] Perfiles, recuperación, ΔRMSE correcto, omega-lambda
  379. print(f"\n[8] Perfiles, recuperación, ΔRMSE, omega-lambda...")
  380. model_b = best_fit["model"]; decay_b = best_fit["decay"]; fix_r_b = best_fit["fix_rescues"]
  381. names_b, lo_b, hi_b = parameter_spec(model_b, decay_b, fix_r_b)
  383. profiles = {}
  384. for pname in ["d_C", "theta_r"]:
  385. if pname not in names_b: continue
  386. idx_p = names_b.index(pname)
  387. p_range = np.linspace(lo_b[idx_p] + 0.01, hi_b[idx_p] - 0.02, 40)
  388. print(f" Perfil {pname}...")
  389. prof = fitting.likelihood_profile(pname, p_range, best_fit,
  390. data["cal_arr"], data["cell_obs"],
  391. cell_weight=CONFIG["CELL_WEIGHT"],
  392. n_starts_refit=5, seed=CONFIG["SEED"] + 3)
  393. profiles[pname] = prof
  394. prof.to_csv(TABLE_DIR / f"profile_{pname}.csv", index=False)
  395. if prof["in_ci95"].any():
  396. print(f" [{prof.loc[prof['in_ci95'], 'param_value'].min():.3f},"
  397. f"{prof.loc[prof['in_ci95'], 'param_value'].max():.3f}]")
  399. # Recuperación cuadrícula
  400. dc_grid = [0.20, 0.35, 0.50, 0.65, 0.80]
  401. frames = []; base_params = best_fit["parameters"].copy()
  402. for dc_true in dc_grid:
  403. rdf = fitting.synthetic_recovery(
  404. true_params={**base_params, "d_C": dc_true}, model=model_b,
  405. arrays_template=data["cal_arr"], cell_obs_template=data["cell_obs"],
  406. n_sim=CONFIG["N_SIM_REC"] // len(dc_grid), n_starts=5,
  407. cell_weight=1.0, decay=decay_b, fix_rescues=fix_r_b,
  408. seed=CONFIG["SEED"] + 4 + int(dc_true * 100))
  409. frames.append(rdf)
  410. rec_df = pd.concat(frames, ignore_index=True)
  411. rec_df.to_csv(TABLE_DIR / "synthetic_recovery_raw.csv", index=False)
  412. fitting.recovery_summary(rec_df, ["d_C"]).to_csv(TABLE_DIR / "tableS4_recovery_summary.csv", index=False)
  413. fitting.recovery_summary_detailed(rec_df).to_csv(TABLE_DIR / "tableS4b_recovery_detailed.csv", index=False)
  415. # ΔRMSE correcto: predicciones fijas, índices remuestreados (N=1000)
  416. print(" Bootstrap ΔRMSE (N=1000, predicciones fijas)...")
  417. delta_rmse = fitting.bootstrap_rmse_difference(
  418. data["val"], data["val_arr"], best_fit,
  419. cal_arr=data["cal_arr"], n_boot=1000, seed=55)
  420. (TABLE_DIR / "delta_rmse_bootstrap.json").write_text(json.dumps(delta_rmse, indent=2))
  421. print(f" ΔRMSE: [{delta_rmse['ci_low']:.4f},{delta_rmse['ci_high']:.4f}] "
  422. f"P(+)={delta_rmse['pct_positive']:.2f} incl.cero={delta_rmse['ci_low']<0}")
  424. # Omega-lambda
  425. omega_map = None
  426. p11 = data["struct"][data["struct"]["protocol_id"] == "P11_selective_opto_retrieval"]
  427. g1 = p11[p11["condition_label"] == "Group 1"]["spi_for_model"].dropna().values
  428. g2 = p11[p11["condition_label"] == "Group 2"]["spi_for_model"].dropna().values
  429. if len(g1) >= 3 and len(g2) >= 3:
  430. omega_map = omega_lambda_compatibility_map(
  431. best_fit["result"].x, model_b, g1, g2,
  432. decay=decay_b, fix_rescues=fix_r_b, n_omega=80, n_lam=80)
  433. print(f" omega-lambda: {omega_map['n_compatible']}/6400 compatible")
  434. (TABLE_DIR / "omega_lambda_map.json").write_text(json.dumps({
  435. "obs_mean": omega_map["obs_mean"], "obs_se": omega_map["obs_se"],
  436. "n_compatible": omega_map["n_compatible"], "note": omega_map["note"]}, indent=2))
  438. # [9] Figuras
  439. print("\n[9] Figuras...")
  440. plots.fig1_behavior(data["memory"])
  441. plots.fig2_cellular(data["cell_raw"])
  442. plots.fig3_identifiability(data["cell_obs"])
  443. plots.fig4_model_comparison(ic)
  444. plots.fig5_calibration(cal_out, cal_group)
  445. plots.fig6_validation(val_out, group_metrics, val_metrics, lopo_df)
  446. plots.fig7_latent_trajectories(best_fit)
  447. plots.fig8_two_memories(best_fit, omega_map)
  448. plots.figS1_bootstrap(boot_df, ci_df)
  449. if profiles: plots.figS2_likelihood_profiles(profiles)
  450. if not rec_df.empty: plots.figS3_synthetic_recovery(rec_df, ["d_C"])
  452. # [10] Summary y verificación
  453. elapsed = time.time() - t0
  454. dc_ind_ci = ci_ind[ci_ind["parameter"] == "d_C"]
  455. summary = {
  456. "version": "V5.1",
  457. "run_date": datetime.datetime.now().isoformat(),
  458. "data_sha256": data_hash, "code_sha256": code_hash,
  459. "config": {k: _s(v) for k, v in CONFIG.items()},
  460. "equivalence_test_passed": True,
  461. # Individual (PRIMARY)
  462. "individual_model": {
  463. "best_key": "access_IND",
  464. "d_C_MLE": float(best_fit_ind["parameters"].get("d_C", 0)),
  465. "d_C_boot_mean": float(dc_ind_ci["mean"].values[0]) if len(dc_ind_ci) else None,
  466. "d_C_boot_ci": [float(dc_ind_ci["ci_low"].values[0]),
  467. float(dc_ind_ci["ci_high"].values[0])] if len(dc_ind_ci) else None,
  468. "AICc_corrected": float(ic_ind_corr.loc[ic_ind_corr["key"] == "access_IND",
  469. "AICc_corrected"].values[0]),
  470. "delta_AICc_mixed": float(ic_ind_corr.loc[ic_ind_corr["key"] == "mixed_IND",
  471. "delta_AICc"].values[0]),
  472. "delta_AICc_storage": float(ic_ind_corr.loc[ic_ind_corr["key"] == "storage_IND",
  473. "delta_AICc"].values[0]),
  474. },
  475. # Aggregate (SENSITIVITY)
  476. "aggregate_model": {
  477. "best_key": best_key,
  478. "d_C_MLE": float(best_fit["parameters"].get("d_C", 0)),
  479. "d_C_boot_ci": [float(ci_df.loc[ci_df["parameter"] == "d_C", "ci_low"].values[0]),
  480. float(ci_df.loc[ci_df["parameter"] == "d_C", "ci_high"].values[0])],
  481. },
  482. # Model averaging (aggregate fixR models)
  483. "model_averaging_aggregate": {
  484. "d_E_MA": ma["d_E_MA"], "d_C_MA": ma["d_C_MA"],
  485. "ci_lo_unconditional": ma["ci_lo_unconditional"],
  486. "ci_hi_unconditional": ma["ci_hi_unconditional"],
  487. },
  488. # Pseudocount
  489. "zero_sensitivity": zero_sens,
  490. # Predictive check
  491. "val_metrics": {k: _s(v) for k, v in val_metrics.items()},
  492. "delta_rmse_bootstrap": {k: _s(v) for k, v in delta_rmse.items()},
  493. "lopo_rmse_valid": float(lopo_valid["rmse"].mean()),
  494. "omega_compatible": omega_map["n_compatible"] if omega_map else None,
  495. "elapsed_s": elapsed,
  496. }
  497. (MODEL_DIR / "analysis_summary.json").write_text(json.dumps(summary, indent=2))
  498. _verify_sync(CONFIG)
  500. print(f"\n{'='*65}")
  501. print(f" V5.1 completado en {elapsed:.0f}s")
  502. print(f" Individual: d_C MLE={summary['individual_model']['d_C_MLE']:.3f} "
  503. f"ΔAICc_storage={summary['individual_model']['delta_AICc_storage']:.2f}")
  504. print(f" Aggregate: d_C MLE={summary['aggregate_model']['d_C_MLE']:.3f}")
  505. print(f" MA: d_C={summary['model_averaging_aggregate']['d_C_MA']:.3f} "
  506. f"CI=[{summary['model_averaging_aggregate']['ci_lo_unconditional']:.3f},"
  507. f"{summary['model_averaging_aggregate']['ci_hi_unconditional']:.3f}]")
  508. print(f" ΔRMSE CI=[{delta_rmse['ci_low']:.4f},{delta_rmse['ci_high']:.4f}] incl.0={delta_rmse['ci_low']<0}")
  509. print(f" Val RMSE={val_metrics['rmse']:.4f}")
  510. print("=" * 65)
  511. return summary
  514. if __name__ == "__main__":
  515. main()

run_analysis.py, under MIT · at the source

Overview

Authors: Adithya Sarma1, Camilla Paraciani1, Junfei Cao1, Evgeniya Tyumeneva1, Caterina Stacchiola1, Elroy L. Meijer1, Nienke de Vries1, Soraya Smit1, Fleur Meijer1, Marit Bonne1, Jean-Christophe Billeter2, Peter Meerlo1, Robbert Havekes1
ORCID iDs: Adithya Sarma, Camilla Paraciani, Evgeniya Tyumeneva, Caterina Stacchiola, Elroy L. Meijer, Nienke de Vries, Soraya Smit, Jean-Christophe Billeter, Peter Meerlo, Robbert Havekes
  1. Neurobiology Expertise Group, Groningen Institute for Evolutionary Life Sciences (GELIFES), University of Groningen, Groningen, Netherlands
  2. Evolutionary Genetics, Development & Behaviour Group, Groningen Institute for Evolutionary Life Sciences (GELIFES), University of Groningen, Groningen, Netherlands
Institutions: University of Groningen (Netherlands)
Journal: Science advances, volume 12, issue 24, article eadu9805
Dates: received 29 November 2024; accepted 1 May 2026; published online 10 June 2026; in print June 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1126/sciadv.adu9805 · PMID 42268954 · PMCID PMC13251821 · OpenAlex W7164130601
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: mouse (organism)
Methods: Statistics, fMRI & imaging, Smoothing, state filtering, decompositions
MeSH: Memory, Long-Term*, Recognition, Psychology*, Sleep Deprivation*, Social Behavior*, Amnesia, Animals, Dentate Gyrus, Male, Mice, Optogenetics (* major topic)
Journal subjects: Neuroscience
Topic: Sleep and Wakefulness Research (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: University of Groningen (Adaptive Life Program Fellowship); Nederlandse Organisatie voor Wetenschappelijk Onderzoek (OCENW.M.23.219); United States Air Force Office of Scientific Research (FA9550-21-1-0310)
Citations: not cited yet (Europe PMC); 68 references in the paper

Abstract

Long-term social memories are vital for forming and maintaining relationships, and social amnesia can disrupt daily life. Using a new paradigm to study these hippocampus-dependent memories, we show that mice can distinguish between multiple social experiences even when they occur in the same context across days. Sleep deprivation immediately after socialization disrupts memory consolidation, leading to social amnesia. Treatment with the Food and Drug Administration–approved phosphodiesterase type 4 inhibitor roflumilast during sleep deprivation protects memory consolidation, while administration immediately before testing reverses social amnesia temporarily. Optogenetic reactivation of dentate gyrus engram cells restores social memory access and enables selective retrieval of individual memories. In addition, using two cFos-based engram-tagging strategies, we find that sleep deprivation selectively affects reactivation of experience-specific engrams while not affecting engram formation or overlap. These results suggest that impaired engram reactivation contributes to sleep deprivation–induced social amnesia and highlights the role of the hippocampal dentate gyrus in maintaining and distinguishing social experiences in a single context.

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

Repositories

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

Zenodo 20786379

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: DataCite
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (8 files), pandas (8 files), SciPy (6 files), Matplotlib (2 files), statsmodels (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
13 files
At the source: doi.org/10.5281/zenodo.20786379

Zenodo 20786380

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: DataCite
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (8 files), pandas (8 files), SciPy (6 files), Matplotlib (2 files), statsmodels (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
13 files
At the source: doi.org/10.5281/zenodo.20786380

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

Tracing map

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

What the map holds:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 22 scripts, each with its path and the digest of its content;
  • 6 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, 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. This study did not generate new materials. The datasets generated and codes used during this study are made publicly available on Dataverse-NL (permalink: https://doi.org/10.34894/0VQF4B) at the date of acceptance.

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, 13 authors, 10 MeSH terms, 3 funders, 65 references.

Cite

This paper

Sarma, A., Paraciani, C., Cao, J., Tyumeneva, E., Stacchiola, C., Meijer, E. L., de Vries, N., Smit, S., Meijer, F., Bonne, M., Billeter, J.-C., Meerlo, P., & Havekes, R. (2026). Restoring access to long-term social recognition memories disrupted by sleep deprivation. Science advances, 12(24), eadu9805. https://doi.org/10.1126/sciadv.adu9805

BibTeX

@article{sarma2026restoring,
author = {Sarma, Adithya and Paraciani, Camilla and Cao, Junfei and Tyumeneva, Evgeniya and Stacchiola, Caterina and Meijer, Elroy L. and de Vries, Nienke and Smit, Soraya and Meijer, Fleur and Bonne, Marit and Billeter, Jean-Christophe and Meerlo, Peter and Havekes, Robbert},
title = {{Restoring access to long-term social recognition memories disrupted by sleep deprivation}},
journal = {Science advances},
year = {2026},
month = jun,
volume = {12},
number = {24},
pages = {eadu9805},
publisher = {American Association for the Advancement of Science},
issn = {2375-2548},
doi = {10.1126/sciadv.adu9805},
url = {https://doi.org/10.1126/sciadv.adu9805},
pmid = {42268954},
pmcid = {PMC13251821}
}

RIS

TY - JOUR
AU - Sarma, Adithya
AU - Paraciani, Camilla
AU - Cao, Junfei
AU - Tyumeneva, Evgeniya
AU - Stacchiola, Caterina
AU - Meijer, Elroy L.
AU - de Vries, Nienke
AU - Smit, Soraya
AU - Meijer, Fleur
AU - Bonne, Marit
AU - Billeter, Jean-Christophe
AU - Meerlo, Peter
AU - Havekes, Robbert
TI - Restoring access to long-term social recognition memories disrupted by sleep deprivation
T2 - Science advances
J2 - Sci Adv
PY - 2026
DA - 2026/06/10
VL - 12
IS - 24
SP - eadu9805
SN - 2375-2548
PB - American Association for the Advancement of Science
DO - 10.1126/sciadv.adu9805
UR - https://doi.org/10.1126/sciadv.adu9805
LA - en
ER -

CSL-JSON

{
"id": "10.1126/sciadv.adu9805",
"type": "article-journal",
"title": "Restoring access to long-term social recognition memories disrupted by sleep deprivation",
"container-title": "Science advances",
"author": [
{
"family": "Sarma",
"given": "Adithya"
},
{
"family": "Paraciani",
"given": "Camilla"
},
{
"family": "Cao",
"given": "Junfei"
},
{
"family": "Tyumeneva",
"given": "Evgeniya"
},
{
"family": "Stacchiola",
"given": "Caterina"
},
{
"family": "Meijer",
"given": "Elroy L."
},
{
"family": "de Vries",
"given": "Nienke"
},
{
"family": "Smit",
"given": "Soraya"
},
{
"family": "Meijer",
"given": "Fleur"
},
{
"family": "Bonne",
"given": "Marit"
},
{
"family": "Billeter",
"given": "Jean-Christophe"
},
{
"family": "Meerlo",
"given": "Peter"
},
{
"family": "Havekes",
"given": "Robbert"
}
],
"container-title-short": "Sci Adv",
"volume": "12",
"issue": "24",
"page": "eadu9805",
"DOI": "10.1126/sciadv.adu9805",
"PMID": "42268954",
"PMCID": "PMC13251821",
"ISSN": "2375-2548",
"publisher": "American Association for the Advancement of Science",
"URL": "https://doi.org/10.1126/sciadv.adu9805",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
10
]
]
}
}

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.1186/s13041-026-01298-y
Repeated sleep deprivation selectively reactivates hippocampal CA1 pyramidal neurons.
Journal: Molecular brain
In common: mouse, 6 references
[2] doi:10.1038/s41467-026-75162-x [code]
Dyscoordination of thalamic reticular spindles is associated with social memory deficits in mice and humans with autism spectrum disorder.
Journal: Nature communications
In common: mouse, 6 references
[3] doi:10.3389/fnbeh.2026.1895371 [code]
Discrete and differentiated encoding of distinct components of spatial experience occurs within the proximodistal subfields of the dorsoventral axis of the hippocampus.
Journal: Frontiers in behavioral neuroscience
In common: SciPy, Matplotlib, NumPy, 3 references
[4] doi:10.7554/elife.108593 [code]
A stress-activated neuronal ensemble in the supramammillary nucleus produces anxiety-like behavior in male mice.
Journal: eLife
In common: mouse, 4 references
[5] doi:10.1038/s41386-026-02492-1 [code]
Proximity in mice induced by an auditory-conditioned stimulus.
Journal: Neuropsychopharmacology : official publication of the American College of Neuropsychopharmacology
In common: pandas, SciPy, Matplotlib, 1 other tool, mouse, 2 references
[6] doi:10.1002/advs.202518450 [code]
A Brain-Wide Atlas of Astrocytic Oxytocin Receptors Reveals a Glial Basis for Nucleus Accumbens Modulation of Affiliative Behavior.
Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany)
In common: statsmodels, pandas, SciPy, 2 other tools, mouse, 1 reference
[7] doi:10.1073/pnas.2517639123 [code]
Transformations of the spatial activity manifold convey aversive information in CA3.
Journal: Proceedings of the National Academy of Sciences of the United States of America
In common: statsmodels, pandas, SciPy, 2 other tools, mouse, 1 reference
[8] doi:10.1126/sciadv.adv3557 [code]
GRASPing experience-dependent protein expression signatures enriched for hippocampal engram cell synapses.
Journal: Science advances
In common: pandas, mouse, 3 references
[9] doi:10.1073/pnas.2530482123 [code]
The detection of episodic memory in others biases social choice.
Journal: Proceedings of the National Academy of Sciences of the United States of America
In common: pandas, 3 references
[10] doi:10.1038/s41593-026-02292-2
Individual differences in prosocial learning are represented in the hippocampal dorsal CA1.
Journal: Nature neuroscience
In common: mouse, 3 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.