OSCR

A generative framework for estimating brain energy landscapes from continuous neuroimaging observations.

Code ↔ Paper

7 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 7 matches · 2 of them tie a paragraph to a whole file without given lines: weak matches, whose lines are not tinted
  1. [1] § Methods › EEG leadfield observation model ↔ pmem_leadfield.py, lines 1 to 36 · score 0.71 · sensor noise, sensor space, compound symmetry, forward model, maps, mixing · Report a mismatch
  2. [2] § Methods › EEG leadfield observation model ↔ pmem_leadfield.py, lines 116 to 171 · score 0.69 · signal subtracted residuals, source signal, compound symmetry, bulk, raw, fit · Report a mismatch
  3. [3] § Methods › Experimental datasets and target network ↔ __init__.py, the whole file · a weak match · score 0.66 · continuous neuroimaging, source space, forward model, pMEM, MEG, simulation · Report a mismatch
  4. [4] § Results › Parameter recovery under cross-channel residual correlation › Recovery evaluation at a representative correlation level ↔ pmem_obs.py, lines 133 to 158 · score 0.65 · cross channel noise, GMM thresholding, full model, correlation, demean, node · Report a mismatch
  5. [5] § Methods › The EM algorithm › Initialization and convergence ↔ pmem_leadfield.py, lines 116 to 171 · score 0.63 · subtracted residual, source signal, regressing, median, raw, leadfield · Report a mismatch
  6. [6] § Methods › EEG leadfield observation model ↔ pmem_leadfield.py, lines 1 to 36 · score 0.61 · cortical source, sensor space, binary states, leadfield, EEG, model · Report a mismatch
  7. [7] § Methods › EEG leadfield observation model ↔ __init__.py, the whole file · a weak match · score 0.53 · leadfield matrix, sensor space, observation model, MEG, EEG · Report a mismatch

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 409 lines · 18 KB · no license · 4 matches

  1. """
  2. Hierarchical pMEM for EEG with a lead-field forward model.
  4. Model
  5. -----
  6. z(t) ~ pMEM(h, J) binary state at each of N cortical sources
  7. s(t) = beta0 + beta1 ⊙ z(t) source amplitude
  8. y(t) = L s(t) + eps(t), eps ~ N(0, Sigma_eps)
  9. Sigma_eps = sigma2 * (I + rho*(J_M - I)) compound-symmetry sensor noise (background EEG)
  11. What this file fixes relative to `pmem.pMEM_Leadfield`
  12. ------------------------------------------------------
  13. 1. Mean-field reference for the M-step uses EXACT enumeration of the pMEM
  14. when N <= MAX_EXACT_N, else block-Gibbs marginals. No more chasing a
  15. biased single-mode mean-field fixed point.
  16. 2. Closed-form WLS update for (beta0, beta1) jointly given m, accounting for
  17. the lead-field mixing (i.e. L^T P L is dense). Original code used
  18. element-wise heuristics that ignored cross-source coupling through L.
  19. 3. Gauge / sign-flip protection: enforces the canonical convention "active
  20. state corresponds to higher signal in the sensor-projected source space"
  21. so that latent identifiability is resolved consistently across runs.
  22. 4. Sensor-space noise sigma2 estimated from data.
  23. 5. ELBO tracking, MAP latent extraction, source-projected warm start using a
  24. Tikhonov-regularized minimum-norm inverse.
  25. 6. Ridge prior on beta0 to handle the null space of L (M < N).
  26. """
  27. from __future__ import annotations
  29. import itertools
  30. from dataclasses import dataclass, field
  32. import numpy as np
  33. from scipy.special import expit
  36. MAX_EXACT_N = 14 # 2^14 = 16384 states; beyond this fall back to Gibbs
  39. # ------------------------------------------------------------
  40. # pMEM prior helpers
  41. # ------------------------------------------------------------
  42. def _enumerate_states(N: int) -> np.ndarray:
  43. return np.array(list(itertools.product([0, 1], repeat=N)), dtype=float)
  46. def exact_prior_moments(h, J, Zs=None):
  47. N = len(h)
  48. if Zs is None:
  49. Zs = _enumerate_states(N)
  50. J_up = np.triu(J, 1)
  51. E = Zs @ h + np.einsum('ki,ij,kj->k', Zs, J_up, Zs)
  52. E -= E.max()
  53. w = np.exp(E)
  54. w /= w.sum()
  55. m0 = w @ Zs
  56. m0_outer = (Zs.T * w) @ Zs
  57. return m0, m0_outer
  60. def gibbs_prior_moments(h, J, n_burn=200, n_samples=2000, rng=None):
  61. """Fallback prior moment estimator for N > MAX_EXACT_N."""
  62. if rng is None:
  63. rng = np.random.default_rng()
  64. N = len(h)
  65. z = rng.integers(0, 2, size=N).astype(float)
  66. m1 = np.zeros(N)
  67. m2 = np.zeros((N, N))
  68. for k in range(n_burn + n_samples):
  69. for i in range(N):
  70. field = h[i] + np.dot(J[i], z) - J[i, i] * z[i]
  71. z[i] = float(rng.random() < expit(field))
  72. if k >= n_burn:
  73. m1 += z
  74. m2 += np.outer(z, z)
  75. m1 /= n_samples
  76. m2 /= n_samples
  77. return m1, m2
  80. def make_sym_upper(M):
  81. M = np.triu(M, 1)
  82. M = M + M.T
  83. np.fill_diagonal(M, 0.0)
  84. return M
  87. # ------------------------------------------------------------
  88. # Lead-field pMEM
  89. # ------------------------------------------------------------
  90. @dataclass
  91. class pMEM_Leadfield_v2:
  92. Y: np.ndarray # [T, M]
  93. L: np.ndarray # [M, N]
  94. n_iter: int = 80
  95. lr_h: float = 0.2
  96. lr_J: float = 0.2
  98. sigma2: float | None = None # sensor noise variance; estimated if None
  99. rho: float | None = None # compound-symmetry correlation; estimated if None
  101. use_exact_prior: bool = True
  102. use_exact_posterior: bool = True # exact enumeration for E-step (N<=MAX_EXACT_N)
  103. sign_flip_protection: bool = True
  104. estimate_noise: bool = True
  105. beta0_ridge: float = 1e-2 # Tikhonov on beta0 (handles M<N null space)
  106. seed: int = 0
  108. # populated in fit
  109. h: np.ndarray = field(init=False)
  110. J: np.ndarray = field(init=False)
  111. beta0: np.ndarray = field(init=False)
  112. beta1: np.ndarray = field(init=False)
  113. m: np.ndarray = field(init=False)
  114. elbo_history: list = field(default_factory=list)
  116. def __post_init__(self):
  117. self.T, self.M = self.Y.shape
  118. _, self.N = self.L.shape
  119. rng = np.random.default_rng(self.seed)
  121. # ---- sensor noise: initialize rho from signal-subtracted residuals ----
  122. # Raw Y includes source signal (L s) which is NOT compound-symmetric and
  123. # would bias rho upward. A quick Tikhonov inverse removes the bulk of
  124. # the signal before fitting the compound-symmetry noise model.
  125. if self.sigma2 is None:
  126. Yc = self.Y - self.Y.mean(axis=0, keepdims=True)
  127. LtL = self.L.T @ self.L
  128. lam_init = 0.1 * np.trace(LtL) / max(self.N, 1)
  129. W = np.linalg.solve(LtL + lam_init * np.eye(self.N), self.L.T)
  130. R0 = Yc - (Yc @ W.T) @ self.L.T # [T, M] approx noise residual
  131. S = (R0.T @ R0) / self.T
  132. self.sigma2 = float(np.diag(S).mean())
  133. off = (S.sum() - np.trace(S)) / (self.M * (self.M - 1))
  134. self.rho = float(np.clip(off / (self.sigma2 + 1e-9), -0.3, 0.95))
  135. if self.rho is None:
  136. self.rho = 0.0
  138. self._update_P()
  140. # ---- pMEM init ------------------------------------------------
  141. self.h = np.zeros(self.N)
  142. J0 = 0.01 * rng.standard_normal((self.N, self.N))
  143. self.J = make_sym_upper(J0)
  145. # ---- warm start: MNE-style Tikhonov inverse -------------------
  146. Yc = self.Y - self.Y.mean(axis=0, keepdims=True)
  147. # noise-normalized leadfield gain ≈ scale of L^T y_t
  148. LtL = self.L.T @ self.L
  149. lam_w = 1e-2 * np.trace(LtL) / max(self.N, 1)
  150. W = np.linalg.solve(LtL + lam_w * np.eye(self.N), self.L.T)
  151. S_hat = Yc @ W.T # [T, N] source estimates
  152. med = np.median(S_hat, axis=0, keepdims=True)
  153. z_warm = (S_hat > med).astype(float)
  154. self.m = 0.05 + 0.90 * z_warm
  156. # beta1 init by sign-aware regression of S_hat on z_warm
  157. b1 = np.zeros(self.N)
  158. for i in range(self.N):
  159. zi = z_warm[:, i] - z_warm[:, i].mean()
  160. num = ((S_hat[:, i] - S_hat[:, i].mean()) * zi).sum()
  161. den = (zi ** 2).sum() + 1e-9
  162. b1[i] = num / den
  163. self.beta1 = np.maximum(np.abs(b1), 0.1)
  165. # beta0 closed form below; placeholder
  166. self.beta0 = S_hat.mean(axis=0) - self.beta1 * z_warm.mean(axis=0)
  168. if self.N <= MAX_EXACT_N:
  169. self._Zs = _enumerate_states(self.N)
  170. else:
  171. self._Zs = None
  173. # --------------------------------------------------------
  174. # Precision and projected leadfield quantities
  175. # --------------------------------------------------------
  176. def _update_P(self):
  177. a = max(1.0 - self.rho, 1e-4)
  178. b = self.rho
  179. denom = a * (a + self.M * b)
  180. I = np.eye(self.M)
  181. O = np.ones((self.M, self.M))
  182. self.P = (1.0 / max(self.sigma2, 1e-6)) * ((1.0 / a) * I - (b / denom) * O)
  183. self.LtPL = self.L.T @ self.P @ self.L # [N, N]
  184. self.LtP = self.L.T @ self.P # [N, M]
  185. self.LtPL_diag = np.diag(self.LtPL)
  187. # --------------------------------------------------------
  188. # Posterior natural parameters at time t
  189. # --------------------------------------------------------
  190. def _posterior_nat_params(self, y_t):
  191. # residual to baseline drive Lβ0
  192. r = y_t - self.L @ self.beta0
  193. # ∂/∂z_i of -0.5 (y - Lβ0 - L diag(β1) z)^T P (...) yields:
  194. # β1_i * [L^T P r]_i for the linear term
  195. # -0.5 * β1_i^2 * (L^T P L)_ii for the diagonal quadratic
  196. # -β1_i β1_j (L^T P L)_ij for the cross term in tilde_J
  197. Lp_r = self.LtP @ r # [N]
  198. tilde_h = self.h + self.beta1 * Lp_r \
  199. - 0.5 * (self.beta1 ** 2) * self.LtPL_diag
  200. B1B1 = np.outer(self.beta1, self.beta1)
  201. cross = B1B1 * self.LtPL
  202. np.fill_diagonal(cross, 0.0)
  203. tilde_J = self.J - cross
  204. return tilde_h, tilde_J
  206. # --------------------------------------------------------
  207. # E-step (mean-field OR exact enumeration)
  208. #
  209. # When exact posterior is used, we also store the pairwise posterior
  210. # <z_i z_j> per time step so the M-step J update is not biased by the
  211. # mean-field independence assumption.
  212. # --------------------------------------------------------
  213. def E_step(self, sweeps: int = 20):
  214. if self.use_exact_posterior and self.N <= MAX_EXACT_N:
  215. self._E_step_exact()
  216. else:
  217. self._E_step_meanfield(sweeps=sweeps)
  219. def _E_step_meanfield(self, sweeps: int = 20):
  220. # Sweep i = N-1, N-2, ..., 0 (downward) so that j > i nodes (already
  221. # updated in this sweep) feed into node i's field — consistent with the
  222. # manuscript's sequential update direction.
  223. for t in range(self.T):
  224. tilde_h, tilde_J = self._posterior_nat_params(self.Y[t])
  225. m_t = self.m[t].copy()
  226. for _ in range(sweeps):
  227. m_old = m_t.copy()
  228. for i in range(self.N - 1, -1, -1):
  229. field = tilde_h[i] + np.dot(tilde_J[i], m_t) # tilde_J diagonal is 0
  230. m_t[i] = expit(field)
  231. if np.max(np.abs(m_t - m_old)) < 1e-6:
  232. break
  233. self.m[t] = np.clip(m_t, 1e-4, 1.0 - 1e-4)
  234. # mean-field has no exact pairwise moments stored
  235. self._m_outer_post = None
  237. def _E_step_exact(self):
  238. """
  239. Vectorised exact E-step: no Python loop over T.
  241. tilde_J is constant across time (depends only on fixed params).
  242. tilde_h(t) varies linearly with Y[t]. We batch all T at once:
  243. E_all[k,t] = z_k · tilde_h(t) + z_k · tJ_up · z_k (K×T)
  244. then softmax over k, and accumulate m_post and m_outer_post via
  245. matrix multiplies — loop-free and GPU-ready.
  246. """
  247. Zs = self._Zs # (K, N), K=2^N
  249. # ---- constant part: tilde_J and its quadratic contribution ----
  250. cross = np.outer(self.beta1, self.beta1) * self.LtPL
  251. np.fill_diagonal(cross, 0.0)
  252. tJ_up = np.triu(self.J - cross, 1) # (N, N) constant
  253. E_quad = np.einsum('ki,ij,kj->k', Zs, tJ_up, Zs) # (K,) constant
  255. # ---- batch tilde_h over all T ----
  256. # r(t) = Y[t] - L @ beta0 → R = Y - L@beta0 shape (T, M)
  257. # LtP @ r(t) = LtP @ R[t] → batch: R @ LtP.T shape (T, N)
  258. R = self.Y - (self.L @ self.beta0) # (T, M)
  259. LtP_R = R @ self.LtP.T # (T, N)
  260. h_const = self.h - 0.5 * (self.beta1 ** 2) * self.LtPL_diag # (N,)
  261. tilde_H = h_const + self.beta1 * LtP_R # (T, N)
  263. # ---- energy matrix ----
  264. E_lin = Zs @ tilde_H.T # (K, T)
  265. E_all = E_lin + E_quad[:, None] # (K, T)
  266. E_all -= E_all.max(axis=0) # numerical stability
  268. # ---- softmax weights ----
  269. w_all = np.exp(E_all) # (K, T)
  270. w_all /= w_all.sum(axis=0) # normalise over K
  272. # ---- posterior statistics ----
  273. m_post = w_all.T @ Zs # (T, N)
  274. w_sum = w_all.sum(axis=1) # (K,) Σ_t w[k,t]
  275. m_outer_post = (Zs.T * w_sum) @ Zs / self.T # (N, N)
  277. self.m = np.clip(m_post, 1e-4, 1.0 - 1e-4)
  278. self._m_outer_post = m_outer_post # for M-step
  280. # --------------------------------------------------------
  281. # Gauge protection: canonical "active ⇒ stronger projected source"
  282. # --------------------------------------------------------
  283. def _fix_sign_flips(self):
  284. # use Tikhonov-projected source estimate as the reference signal per node
  285. Yc = self.Y - self.Y.mean(axis=0, keepdims=True)
  286. LtL = self.L.T @ self.L
  287. lam_w = 1e-2 * np.trace(LtL) / max(self.N, 1)
  288. S_hat = Yc @ np.linalg.solve(LtL + lam_w * np.eye(self.N), self.L.T).T
  289. for i in range(self.N):
  290. if self.m[:, i].std() < 1e-6 or S_hat[:, i].std() < 1e-6:
  291. continue
  292. c = np.corrcoef(self.m[:, i], S_hat[:, i])[0, 1]
  293. if np.isnan(c) or c >= -0.05:
  294. continue
  295. self.m[:, i] = 1.0 - self.m[:, i]
  296. self.beta0[i] = self.beta0[i] + self.beta1[i]
  297. self.beta1[i] = -self.beta1[i]
  298. neg = self.beta1 < 0
  299. if np.any(neg):
  300. self.beta1[neg] = -self.beta1[neg]
  301. self.m[:, neg] = 1.0 - self.m[:, neg]
  302. self.beta0[neg] = self.beta0[neg] - self.beta1[neg]
  303. self.beta1 = np.maximum(self.beta1, 0.05)
  305. # --------------------------------------------------------
  306. # M-step
  307. # --------------------------------------------------------
  308. def M_step(self):
  309. # ---- (beta0, beta1): two alternating closed-form sweeps ----------
  310. # Both updates are optimal given the other fixed; two sweeps of the
  311. # outer EM loop are enough for practical convergence.
  312. for _ in range(2):
  313. # β0 — Tikhonov-regularized source-space WLS (δI handles null(L))
  314. r_bar = self.Y.mean(axis=0) - self.L @ (self.beta1 * self.m.mean(axis=0))
  315. A0 = self.LtPL + self.beta0_ridge * np.eye(self.N)
  316. self.beta0 = np.linalg.solve(A0, self.LtP @ r_bar)
  318. # β1 — N×N linear system A β1 = b, A_{kk'} = K_{kk'} * sum_t C_{kk'}(t)
  319. # K = L^T P L. Diagonal exact: C_{kk}(t) = m_k(t) (binary identity).
  320. # Off-diagonal: exact posterior E_q[z_k z_k'] when available (enumeration),
  321. # mean-field m_k m_k' + diagonal V(t) correction otherwise.
  322. R0 = self.Y - (self.L @ self.beta0)[None, :] # [T, M]
  323. G = R0 @ self.LtP.T # [T, N]
  324. numer = (self.m * G).sum(axis=0) # [N]
  325. if self._m_outer_post is not None:
  326. sum_C = self._m_outer_post * self.T # [N, N]: sum_t C(t)
  327. A1 = self.LtPL * sum_C # diagonal already exact
  328. else:
  329. mm = self.m.T @ self.m # [N, N]: sum_t m_i m_j
  330. sum_V = (self.m * (1.0 - self.m)).sum(axis=0) # [N]
  331. A1 = self.LtPL * mm + np.diag(self.LtPL_diag * sum_V)
  332. A1 += 1e-9 * np.eye(self.N)
  333. self.beta1 = np.maximum(np.linalg.solve(A1, numer), 0.05)
  335. # ---- h and J moment matching ---------------------------------
  336. m_bar = self.m.mean(axis=0)
  337. # Use exact posterior pairwise <z_i z_j> if available; otherwise fall
  338. # back to the (biased) mean-field outer product.
  339. if getattr(self, "_m_outer_post", None) is not None:
  340. m_outer = self._m_outer_post
  341. else:
  342. m_outer = (self.m.T @ self.m) / self.T
  344. if self.use_exact_prior and self.N <= MAX_EXACT_N:
  345. m0, m0_outer = exact_prior_moments(self.h, self.J, Zs=self._Zs)
  346. else:
  347. m0, m0_outer = gibbs_prior_moments(self.h, self.J)
  348. self.h += self.lr_h * (m_bar - m0)
  349. grad_J = np.triu(m_outer - m0_outer, 1)
  350. self.J += self.lr_J * grad_J
  351. self.J = make_sym_upper(self.J)
  353. # ---- sensor noise update (compound-symmetry) ------------------
  354. # E[r̄ r̄^T] = Σ_ε + L diag(β₁) Cov_q[z] diag(β₁) L^T
  355. # Subtract the sensor-space posterior-variance contamination to recover Σ_ε.
  356. if self.estimate_noise:
  357. S_drive = self.beta0[None, :] + self.beta1[None, :] * self.m
  358. mu = S_drive @ self.L.T # [T, M]
  359. R = self.Y - mu # [T, M]
  360. S = (R.T @ R) / self.T # [M, M]
  361. if self._m_outer_post is not None:
  362. avg_cov_q = self._m_outer_post - (self.m.T @ self.m) / self.T # [N, N]
  363. else:
  364. avg_cov_q = np.diag((self.m * (1.0 - self.m)).mean(axis=0)) # diagonal
  365. B = self.L * self.beta1[None, :] # [M, N]: L diag(β₁)
  366. S = S - B @ avg_cov_q @ B.T # [M, M]
  367. diag_var = max(float(np.diag(S).mean()), 1e-4)
  368. off_cov = float((S.sum() - np.trace(S)) / (self.M * (self.M - 1)))
  369. sigma2_new = diag_var
  370. rho_new = float(np.clip(off_cov / (sigma2_new + 1e-9), -0.3, 0.95))
  371. self.sigma2 = 0.5 * self.sigma2 + 0.5 * sigma2_new
  372. self.rho = 0.5 * self.rho + 0.5 * rho_new
  373. self._update_P()
  375. # --------------------------------------------------------
  376. # ELBO (exact prior log-Z when N small)
  377. # --------------------------------------------------------
  378. def _elbo(self):
  379. m_bar = self.m.mean(axis=0)
  380. m_outer = (self.m.T @ self.m) / self.T
  381. if self._Zs is not None:
  382. E = self._Zs @ self.h + np.einsum('ki,ij,kj->k', self._Zs, np.triu(self.J, 1), self._Zs)
  383. logZ = E.max() + np.log(np.exp(E - E.max()).sum())
  384. else:
  385. logZ = 0.0
  386. prior_term = self.T * (m_bar @ self.h
  387. + np.einsum('ij,ij->', np.triu(self.J, 1), m_outer)
  388. - logZ)
  389. mu = (self.beta0 + self.beta1 * self.m) @ self.L.T
  390. R = self.Y - mu
  391. obs = -0.5 * np.einsum('ti,ij,tj->', R, self.P, R)
  392. eps = 1e-12
  393. mc = np.clip(self.m, eps, 1 - eps)
  394. H = -(mc * np.log(mc) + (1 - mc) * np.log(1 - mc)).sum()
  395. return float(prior_term + obs + H)
  397. # --------------------------------------------------------
  398. def fit(self, verbose=False):
  399. for it in range(self.n_iter):
  400. self.E_step()
  401. if self.sign_flip_protection and (it < 5 or it % 20 == 0):
  402. self._fix_sign_flips()
  403. self.M_step()
  404. if verbose and ((it + 1) % 10 == 0 or it == self.n_iter - 1):
  405. elbo = self._elbo()
  406. self.elbo_history.append(elbo)
  407. print(f"[pMEM_Leadfield_v2] it {it+1}/{self.n_iter} ELBO={elbo:.2f} "
  408. f"σ²={self.sigma2:.3f} ρ={self.rho:.3f}")
  409. return self.h, self.J, self.beta0, self.beta1, self.m

pmem_leadfield.py at commit b1ca5dd, no license · at the source

Overview

Authors: Seok-Oh Jeong1, Euisun Kim2,3, Jiyoung Kang4, Jinseok Eo2,3, Dongmyeong Lee3, Hae-Jeong Park2,3,5,6
ORCID iDs: Euisun Kim, Hae-Jeong Park
  1. Department of Statistics, Hankuk University of Foreign Studies, Yong-In, Republic of Korea
  2. Graduate School of Medical Science, Brain Korea 21 Project, Department of Nuclear Medicine, Department of Psychiatry, Yonsei University College of Medicine, Seoul, Republic of Korea
  3. Center for Systems and Translational Brain Science, Institute of Human Complexity and Systems Science, Yonsei University, Seoul, Republic of Korea
  4. Department of Scientific Computing, Pukyong National University, Busan, Republic of Korea
  5. Department of Cognitive Science, Yonsei University, Seoul, Republic of Korea
  6. Brain Research Institute, Epilepsy Research Institute, Yonsei University College of Medicine, Seoul, Republic of Korea
Institutions: Hankuk University of Foreign Studies (South Korea); Yonsei University (South Korea); Yonsei University Health System (South Korea); Pukyong National University (South Korea)
Journal: Imaging neuroscience (Cambridge, Mass.), volume 4, article IMAG.a.1415
Dates: received 13 June 2026; accepted 22 September 2026; published online 27 October 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1162/imag.a.1415 · PMID 42953883 · PMCID PMC13716286 · OpenAlex W7220429905
Open access: diamond, a free copy (OpenAlex)
Status: code verified (how the code is found and verified)
Categories: EEG (modality), MEG (modality), human (organism), methods / tools (subfield)
Methods: Spectral & time-frequency, Connectivity, fMRI & imaging, Smoothing, state filtering, decompositions
Keywords: pairwise maximum entropy model, nonlinear state dynamics, latent binary states, correlated noise, EEG/MEG source, energy landscape
MeSH: Brain*, Models, Neurological*, Neuroimaging*, Electroencephalography, Entropy, Humans, Magnetoencephalography (* major topic)
Topic: Neural dynamics and brain function (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Citations: not cited yet (Europe PMC); 27 references in the paper
Research resources: RRID:SCR_007345

Abstract

Pairwise maximum entropy models (pMEMs) are gaining increasing attention as a principled way to study nonlinear brain-state transitions. By representing each brain region as a binary active/inactive state and estimating the least-biased probability distribution over all such network configurations, pMEMs reconstruct an energy landscape in which local minima correspond to stable or metastable brain states and transitions reflect movement across attractor basins. However, neuroimaging recordings are continuous and spatially mixed, forcing conventional pMEM analyses to demean and hard-threshold each channel. This ad hoc binarization discards amplitude information and can conflate spatially correlated noise or source mixing with genuine neural coupling. We introduce a hierarchical generative pMEM that links latent binary network states to continuous observations. Node-specific baseline and gain parameters model amplitude modulation, compound-symmetry noise captures shared residual covariance, and, for electroencephalography/magnetoencephalography (EEG/MEG), the electromagnetic leadfield is included as a generative forward model rather than a preprocessing inverse. Model parameters and latent states are jointly estimated using exact posterior enumeration for the small networks analyzed here and a mean-field extension for larger networks. Controlled simulations demonstrated two complementary advantages of the proposed framework. Residual-covariance modeling improved recovery of pMEM fields, couplings, and latent states by accounting for shared fluctuations, with the clearest gains when residual spatial correlation was sufficiently strong. The generative leadfield model further improved sensor-space inference by avoiding inverse-first source reconstruction, including settings in which covariance modeling alone offered limited gains. Resting-state EEG validation showed higher residual correlation during eyes-closed than eyes-open rest, consistent with elevated and spatially coherent alpha activity. Applied further to propofol sedation EEG (alpha-band amplitude envelopes, 8–13 Hz), the fitted 6-node default mode network (DMN) model showed pipeline-dependent condition effects in residual spatial correlation and basin entropy, not apparent under the conventional demean-based pMEM.

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

psalmhim/pMEM-EEGMEG

License: none: the authors keep all their rights
State: the link answers, verified on 30 October 2026
Evidence: files inventoried
Commit: b1ca5dd2adb28f102a69d9dd0b660b8c2ae2f795, 13 June 2026
Languages: Python (5)
Size: 8 files, 5 scripts
Software Heritage: not archived
Found in: “Data and Code Availability”
Holds: README, environment (requirements.txt)
Not found: license file, CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (3 files), SciPy (2 files), scikit-learn (1 file)
Availability: 1 check, the latest on 30 October 2026: the link answers
  • 30 October 2026: the link answers
6 files
At the source: github.com/psalmhim/pMEM-EEGMEG

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

The PhysioNet EEG Motor Movement/Imagery dataset is freely available at https://physionet.org/content/eegmmidb/1.0.0/ (Goldberger et al., 2000; Schalk, 2009). The MNE sample MEG dataset is distributed with MNE-Python and is available at https://mne.tools/stable/index.html (Gramfort et al., 2014). The propofol sedation EEG dataset (Chennu et al., 2016) is openly available from the University of Cambridge Data Repository at https://doi.org/10.17863/CAM.68959, under a Creative Commons Attribution 2.0 UK: England & Wales license.

The simulation and analysis code is available at https://github.com/psalmhim/pMEM-EEGMEG, implemented in Python (≥3.9) using NumPy (≥1.23), SciPy (≥1.9), scikit-learn (≥1.1), MNE-Python (≥1.3), Nilearn (≥0.10), and Matplotlib (≥3.6); exact package versions are pinned in requirements.txt in the repository.

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, 30 October 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 6 authors, 6 keywords, 7 MeSH terms, 27 references, 1 RRID.

Cite

This paper

Jeong, S.-O., Kim, E., Kang, J., Eo, J., Lee, D., & Park, H.-J. (2026). A generative framework for estimating brain energy landscapes from continuous neuroimaging observations. Imaging neuroscience (Cambridge, Mass.), 4, IMAG.a.1415. https://doi.org/10.1162/imag.a.1415

BibTeX

@article{jeong2026generative,
author = {Jeong, Seok-Oh and Kim, Euisun and Kang, Jiyoung and Eo, Jinseok and Lee, Dongmyeong and Park, Hae-Jeong},
title = {{A generative framework for estimating brain energy landscapes from continuous neuroimaging observations}},
journal = {Imaging neuroscience (Cambridge, Mass.)},
year = {2026},
month = oct,
volume = {4},
pages = {IMAG.a.1415},
publisher = {MIT Press},
issn = {2837-6056},
doi = {10.1162/imag.a.1415},
url = {https://doi.org/10.1162/imag.a.1415},
pmid = {42953883},
pmcid = {PMC13716286}
}

RIS

TY - JOUR
AU - Jeong, Seok-Oh
AU - Kim, Euisun
AU - Kang, Jiyoung
AU - Eo, Jinseok
AU - Lee, Dongmyeong
AU - Park, Hae-Jeong
TI - A generative framework for estimating brain energy landscapes from continuous neuroimaging observations
T2 - Imaging neuroscience (Cambridge, Mass.)
J2 - Imaging Neurosci (Camb)
PY - 2026
DA - 2026/10/27
VL - 4
SP - IMAG.a.1415
SN - 2837-6056
PB - MIT Press
DO - 10.1162/imag.a.1415
UR - https://doi.org/10.1162/imag.a.1415
LA - en
ER -

CSL-JSON

{
"id": "10.1162/imag.a.1415",
"type": "article-journal",
"title": "A generative framework for estimating brain energy landscapes from continuous neuroimaging observations",
"container-title": "Imaging neuroscience (Cambridge, Mass.)",
"author": [
{
"family": "Jeong",
"given": "Seok-Oh"
},
{
"family": "Kim",
"given": "Euisun"
},
{
"family": "Kang",
"given": "Jiyoung"
},
{
"family": "Eo",
"given": "Jinseok"
},
{
"family": "Lee",
"given": "Dongmyeong"
},
{
"family": "Park",
"given": "Hae-Jeong"
}
],
"container-title-short": "Imaging Neurosci (Camb)",
"volume": "4",
"page": "IMAG.a.1415",
"DOI": "10.1162/imag.a.1415",
"PMID": "42953883",
"PMCID": "PMC13716286",
"ISSN": "2837-6056",
"publisher": "MIT Press",
"URL": "https://doi.org/10.1162/imag.a.1415",
"language": "en",
"issued": {
"date-parts": [
[
2026,
10,
27
]
]
}
}

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

Similar papers

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

[1] doi:10.1016/j.patter.2026.101560 [code]
Automating region selection with genetic algorithms for energy landscape analyses of brain dynamics.
Journal: Patterns (New York, N.Y.)
In common: SciPy, NumPy, methods / tools, 6 references
[2] doi:10.1002/advs.202414016 [code]
Comparison of Large Language Model with Aphasia.
Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany)
In common: NumPy, 5 references
[3] doi:10.1038/s41562-025-02121-9 [code]
General anaesthesia decreases the uniqueness of brain functional connectivity across individuals and species.
Journal: Nature human behaviour
In common: scikit-learn, SciPy, NumPy, 5 references
[4] doi:10.1038/s41598-025-32207-3 [code]
Multi-scale EEG feature decoding with Swin Transformers for subject independent motor imagery BCIs.
Journal: Scientific reports
In common: scikit-learn, SciPy, NumPy, physionet.org/content/eegmmidb, methods / tools, EEG, 1 reference
[5] doi:10.1002/hbm.70029 [code]
The alterations of repetitive transcranial magnetic stimulation on the energy landscape of resting-state networks differ across the human cortex.
Journal: Human brain mapping
In common: 5 references
[6] doi:10.1371/journal.pcbi.1010811 [code]
LSD-induced increase of Ising temperature and algorithmic complexity of brain dynamics.
Journal: PLoS computational biology
In common: SciPy, NumPy, 5 references
[7] doi:10.7554/elife.97107 [code]
A dynamic generative model can extract interpretable oscillatory components from multichannel neurophysiological recordings.
Journal: eLife
In common: scikit-learn, SciPy, NumPy, MEG, EEG, 3 references
[8] doi:10.1038/s41598-025-07703-1
Scalable geometric learning with correlation-based functional brain networks.
Journal: Scientific reports
In common: EEG, 2 references, author Hae-Jeong Park
[9] doi:10.3389/fnhum.2022.992649 [code]
Paradoxical markers of conscious levels: Effects of propofol on patients in disorders of consciousness.
Journal: Frontiers in human neuroscience
In common: SciPy, NumPy, EEG, 4 references
[10] doi:10.1186/s12868-024-00854-3 [code]
State-transition dynamics of resting-state functional magnetic resonance imaging data: model comparison and test-to-retest analysis.
Journal: BMC neuroscience
In common: scikit-learn, SciPy, NumPy, MEG, EEG, 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.

Add the badge to the code's README

OSCR: the badge of this paper

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown
<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/doi/10.1162/imag.a.1415.svg)](https://openscicode.org/badge/doi/10.1162/imag.a.1415)
<!-- oscr-badge:end -->
reStructuredText
.. oscr-badge:start

.. image:: https://openscicode.org/badge/doi/10.1162/imag.a.1415.svg
   :target: https://openscicode.org/badge/doi/10.1162/imag.a.1415
   :alt: OSCR

.. oscr-badge:end
HTML
<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/doi/10.1162/imag.a.1415"><img src="https://openscicode.org/badge/doi/10.1162/imag.a.1415.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

On GitHub, open the README, choose Edit (the pencil), paste the snippet under the title, then "Commit changes".

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge 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: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

Versions of its code: the releases of its code tied to a version of this paper (the preprint, the accepted manuscript, the version of record…), each with the tracing map versioned with it.

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