REVIEW 5 major objections 6 minor 26 references
Quantifying Time-Varying Sources in Magnetoencephalography -- A Discrete Approach
T0 review · 5 major / 6 minor · reviewed 2026-08-14 · deepseek-v4-flash
Pith's one-line read This paper turns the time-varying MEG source-localization problem into a finite hidden Markov model, computing a closed-form discrete posterior for the source location at every time point from all measurements.
desk verdict A genuinely new discrete HMM formulation for MEG source localization with exact forward-backward posteriors, but the fixed moment assumption and an unquantified discretization error keep the central claim from being fully supported. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
The machinery is the discretization of the brain region of interest into $K$ voxels $\{V_k\}$ with binary indicator variables $v_{tk}$, which turns the continuous state-space model into a finite hidden Markov model. This discretization makes the posterior $P(J_t^p \in V_k \mid Y_T, \Theta)$ computable in closed form by the forward-backward algorithm, and it gives the EM algorithm closed-form parameter updates. The switch and dynamic procedures are auxiliary devices that make the same machinery tractable when there are multiple sources or when the region of interest is unknown; the switch handles state-space explosion by iterated marginal posteriors, and the dynamic procedure shrinks the ROI and increases mesh density during iterations.
What would settle it
Simulate a source moving along a known trajectory while its dipole moment direction rotates; run the discrete method with a fixed unit moment and check whether the posterior means track the true path or shift to wrong voxels. A cleaner check: with known parameters, compare the discrete posterior to the exact continuous posterior obtained by fine-grid numerical integration on a small ROI; disagreement beyond the chosen voxel resolution would refute the closed-form claim.
Extended reading notes
Core claim
The paper claims that the continuous nonlinear MEG inverse problem can be replaced by a finite hidden Markov model over voxels: for each time point $t$ and voxel $V_k$, the posterior probability $P(J_t^p \in V_k \mid Y_T, \Theta)$ is computed exactly (up to discretization) as the product $\alpha_{tk}(\Theta)\beta_{tk}(\Theta)$ from the forward-backward algorithm. The forward pass uses measurements up to $t$, the backward pass uses measurements from $t+1$ to $T$, so the posterior reflects the whole experiment. Parameters $\Theta$ are estimated by an EM algorithm whose M-step updates have closed forms because the state space is finite. For multiple sources, a switch procedure computes each source's marginal posterior conditional on the others, and a dynamic procedure shrinks the region of interest while refining the mesh across EM iterations. The result is a quantification of the source trajectory as a posterior probability map rather than a point estimate or a sample cloud.
Load-bearing premise
The load-bearing premise is that the dipole moment and strength vector $q_t$ is fixed for all time points, and in the real-data analysis is set to a unit vector; if the true moment changes over time or is misspecified, the computed location posterior can be systematically biased.
Editorial extensions
If this is right
- At every time point the method yields a full voxel-level probability map, so source uncertainty is reported directly rather than as a point estimate.
- Because the backward pass uses future measurements, the posterior can localize a source more accurately after the fact than any filtering-only scheme.
- The closed-form EM updates mean the computational bottleneck is the forward-backward sweep, so the method scales to long recordings without particle sampling.
- For multiple sources, the switch procedure offers a tractable approximation with a proven bound relative to the full joint posterior.
- When no region of interest is given, the dynamic procedure starts from a coarse whole-head grid and automatically concentrates the mesh around the active source.
Reading between the lines
- The same discretization could be extended to a joint location-orientation state by including the moment orientation in the voxel state, at the cost of a larger state space; the forward-backward and EM machinery would carry over unchanged.
- Comparing posteriors from two mesh resolutions on the same data would give a practical convergence diagnostic, since the paper does not quantify discretization error.
- If the number of sources is uncertain, the switch procedure could be embedded in a model-selection scheme, treating the source count as a hyperparameter and comparing marginal likelihoods.
- A misspecified moment should show up as a systematic shift of the posterior; this is a testable empirical prediction that the paper does not address.
Signed reviews
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The paper proposes a discrete hidden Markov model approach to the MEG inverse problem. The source location space is discretized into K voxels, and the hidden state at each time point is the voxel containing the source. A forward-backward algorithm computes the discrete posterior P(J_t^p in V_k | Y_T, Theta) using the entire set of measurements, and an EM algorithm estimates the model parameters. A dynamic procedure iteratively shrinks the ROI and increases mesh resolution, and a switch procedure approximates the multi-source posterior by iteratively conditioning on the other sources. The method is demonstrated on simulated single- and two-source data, on a BCI MEG data set, and on an EEG working-memory data set.
Significance. If the proposed framework were fully validated, it would offer an attractive alternative to sequential Monte Carlo for time-varying MEG source localization, providing closed-form discrete posteriors and smoothing over all measurements. The forward-backward recursions and EM updates are standard under the stated discrete model, and the dynamic and switch procedures address important computational scalability issues. The use of simulated ground truth and two real data applications is a strength. However, several load-bearing issues currently prevent the central claim from being supported: the discrete HMM is not completely specified, the moment parameter is fixed without sensitivity analysis, the discretization error is unquantified, and the switch-procedure theorem is vacuous as stated. These issues are fixable but require substantial additional work.
major comments (5)
- [Section 2.2, Eqs (2.6)-(2.8)] The discrete HMM is not fully specified. The paper states approximations of the continuous initial, transition, and emission probabilities, but never defines how P(v_1k=1|Theta) and P(v_tk=1|v_{t-1,l}=1,Theta) are computed from the Gaussian AR model of Eq (2.4). The forward-backward algorithm in Tables 1-2 and the EM updates in Table 3 require these quantities. Please provide explicit formulas, presumably integrals of the Gaussian densities over the voxels, and state how the categorical initial and transition probabilities are derived from the continuous model.
- [Section 2.2 and Section 4, Eq (2.1)] The moment/strength parameter q_t is fixed for all time points (Section 2.2) and set to a unit vector in the real data (Section 4). Because the emission B_l(J_t^p) depends on q_t, a time-varying or misspecified moment can produce time-varying field patterns that the fixed-q likelihood will attribute to source motion. The simulations in Table 5 keep q_t essentially constant (moment variances 1e-4), so this failure mode is never exercised. Add a sensitivity analysis with rotating or misspecified moments, or explicitly restrict the claim to sources with a known, constant moment.
- [Section 2.2 and Section 2.5] The paper claims that the discrete posterior 'approximates' the continuous posterior, but no discretization error bound or convergence analysis is provided. The dynamic procedure's ROI shrinkage also has no coverage guarantee. Provide either a theoretical bound (e.g., as K increases) or an empirical convergence study (e.g., showing stability of the posterior as mesh grids increase) to support the central approximation claim.
- [Appendix B / Theorem 2.2] Theorem 2.2 as stated is vacuous: 'there exists epsilon > 0 such that |...| <= c epsilon' is trivially true for arbitrary epsilon, and the proof relies on assumption (B.3), which essentially assumes the difference between the switch and non-switch marginals is small. The theorem therefore does not provide a meaningful bound on the switch procedure's approximation error. Either state the theorem with a concrete quantity (e.g., the distance ||Theta_s - Theta_ns||) or present the switch procedure as a heuristic supported by the two-source simulation.
- [Section 4] In the real-data analysis, the parameters mu_0, Sigma_0, and Sigma are manually set, and the reported posterior distributions are not accompanied by any sensitivity analysis over these choices. Since the method's output likely depends on these manual settings, please report robustness results, for example by varying the manual parameters and showing that the posterior means and distributions remain stable.
minor comments (6)
- [Section 2.2, Eq (2.9)] The product alpha_tk(Theta) beta_tk(Theta) equals the posterior P(v_tk=1|Y_T,Theta) only up to a normalizing constant that depends on t. The equality in (2.9) should be replaced by proportionality, or the normalization over k should be stated explicitly.
- [Section 2.1, bullet 1] The first Markov property is printed tautologically as p(J_t^p|J_{t-1}^p,Theta) = p(J_t^p|J_{t-1}^p,Theta); the intended statement is that J_t^p is first-order Markov in time.
- [Section 3, Table 6] The standard deviations in Table 6 are based on only four repetitions. This is a very small number of replications; please either increase the number of simulations or interpret the reported variability with appropriate caution.
- [Abstract and Introduction] The claim that sequential sampling methods 'may suffer from slow convergence due to the large volume of measurement' is not supported by any comparison in the paper. Adding a small benchmark against SIS or a particle filter would substantiate this motivation.
- [Notation] The symbol J_p_T for the collection of sources is easily confused with the single-source notation J_t^p; consider using clearer notation, such as a bold or uppercase calligraphic symbol.
- [Section 2.3 and Section 6] There are minor typos: in Section 2.3, 'P(J_p_T|YT, Theta(J-1)' is missing a closing bracket, and in Section 6, 'mutiple-source' should be 'multiple-source'.
Circularity Check
No significant circularity: the discrete posterior is derived from an explicit HMM and validated against independent simulations; the fixed-moment assumption and self-citations are modeling assumptions and prior work, not circular reductions.
full rationale
The derivation of the discrete posterior is self-contained. Section 2.2 discretizes the continuous spatio-temporal model (2.3)-(2.4) into a finite HMM over voxels, defines P(J_t^p in V_k | Y_T, Theta) through Bayes' rule in (2.5) and (2.9), and computes it with the standard forward-backward recursions in Tables 1-2. The closed-form EM updates in Table 3 maximize the same model's likelihood; the fitted parameters are not defined in terms of the posterior they are used to compute. The fixed-moment assumption (Section 2.2: 'the moment and strength parameter q_t is fixed for all the time points throughout this paper'; Section 4: 'we assume a unit moment for all possible sources for simplicity') conditions the posterior on q rather than inferring q, which is a model limitation that could bias location estimates if q is misspecified, but it is not a circular reduction because q is not fitted from the location posterior. The simulations are validated against independently generated ground truth (Tables 5-7 and Figures 2-4), so the central claim does not collapse into its inputs. The self-citations [23] and [24] are prior published work used for comparison and for the real-data number-of-sources input; the location-posterior derivation does not depend on their validity, so this is normal citation rather than load-bearing circularity.
Assumptions & free parameters
free parameters (4)
- Dipole moment vector q =
q = (3, 3, 3) in simulations; unit moment in real data
- Initial ROI and mesh grid sizes =
K_i = 10 in simulations; ROI around MNE for real data
- Manually set mu_0, Sigma_0, Sigma in real data =
Chosen from MNE and Fourier noise estimates
- ROI shrinking factor (3-sigma) =
3
assumptions (7)
- domain assumption AR(1) Gaussian process for the source state
- domain assumption Gaussian, homogeneous, uncorrelated sensor noise
- domain assumption Quasi-static approximation of Maxwell's equations and the Biot-Savart forward model
- ad hoc to paper Discretization approximates the continuous posterior
- standard math Regularity conditions (C1)-(C7)
- domain assumption A true ROI exists and contains the source for all time points
- ad hoc to paper Switch procedure fixed point approximates the joint posterior marginals
Cite this review
Pith. "Pith review of Quantifying Time-Varying Sources in Magnetoencephalography -- A Discrete Approach." pith.science (2026). https://pith.science/paper/2LRYHMWU
@misc{pith2026190803926,
author = {Pith},
title = {Pith review of: Quantifying Time-Varying Sources in Magnetoencephalography -- A Discrete Approach},
year = {2026},
howpublished = {\url{https://pith.science/paper/2LRYHMWU}},
note = {Machine review of arXiv:1908.03926}
}
read the original abstract
We study the distribution of brain source from the most advanced brain imaging technique, Magnetoencephalography (MEG), which measures the magnetic fields outside the human head produced by the electrical activity inside the brain. Common time-varying source localization methods assume the source current with a time-varying structure and solve the MEG inverse problem by mainly estimating the source moment parameters. These methods use the fact that the magnetic fields linearly depend on the moment parameters of the source, and work well under the linear dynamic system. However, magnetic fields are known to be non-linearly related to the location parameters of the source. The existing work on estimating the time-varying unknown location parameters is limited. We are motivated to investigate the source distribution for the location parameters based on a dynamic framework, where the posterior distribution of the source is computed in a closed form discretely. The new framework allows us not only to directly approximate the posterior distribution of the source current, where sequential sampling methods may suffer from slow convergence due to the large volume of measurement, but also to quantify the source distribution at any time point from the entire set of measurements reflecting the distribution of the source, rather than using only the measurements up to the time point of interest. Both a dynamic procedure and a switch procedure are proposed for the new discrete approach, balancing estimation accuracy and computational efficiency when multiple sources are present. In both simulation and real data, we illustrate that the new method is able to provide comprehensive insight into the time evolution of the sources at different stages of the MEG and EEG experiment.
Reference graph
Works this paper leans on
-
[1]
M. S. Arulampalam, S. Maskell, N. Gordon, and T. Clapp. A tutorial on particle filters for online nonlinear/non-gaussian bayesian tracking. IEEE Transactions on Signal Processing, 50(2):174–188, 2002
work page 2002
-
[2]
S. Baillet and L. Garnero. A bayesian approach to introducing anatomo-functional priors in the eeg/meg inverse problem. IEEE Transactions on Biomedical Engineer- ing, 44(5):374–385, 1997. QUANTIFYING TIME-VARYING SOURCES IN MEG 37
work page 1997
-
[3]
S. Baillet, J. C. Mosher, and R. M. Leahy. Electromagnetic brain mapping. IEEE Signal Processing Magazine, 18(6):14–30, 2001
work page 2001
-
[4]
E. Boto, N. Holmes, J. Leggett, G. Roberts, V. Shah, S. S. Meyer, L. D. Mu˜ noz, K. J. Mullinger, T. M. Tierney, S. Bestmann, R. G. Barnes, R. Bowtell, and J. M. Brookes. Moving magnetoencephalography towards real-world applications with a wearable system. Nature, 555(7698):657, 2018
work page 2018
-
[5]
A. P. Dempster, N. M. Laird, and D. B. Rubin. Maximum likelihood from incom- plete data via the em algorithm. Journal of the Royal Statistical Society. Series B (Methodological), 39(1):1–38, 1977
work page 1977
-
[6]
M. Fukushima, O. Yamashita, T. R. Kn¨ osche, and M.-a. Sato. Meg source reconstruc- tion based on identification of directed source interactions on whole-brain anatomical networks. NeuroImage, 105:408–427, 2015
work page 2015
-
[7]
M. S. H¨ am¨ al¨ ainen, R. Hari, R. J. Ilmoniemi, J. Knuutila, and O. V. Lounasmaa. Magnetoencephalography theory, instrumentation, and applications to noninvasive studies of the working human brain. Reviews of Modern Physics , 65:413–497, 1993
work page 1993
-
[8]
M. S. H¨ am¨ al¨ ainen and R. J. Ilmoniemi. Interpreting magnetic fields of the brain: min- imum norm estimates. Medical and Biological Engineering and Computing, 32(1):35– 42, 1994
work page 1994
Show all 26 references
-
[9]
Lamus, M
C. Lamus, M. S. H¨ am¨ al¨ ainen, S. Temereanca, E. N. Brown, and P. L. Purdon. A spa- tiotemporal dynamic distributed solution to the meg inverse problem. NeuroImage, 63(2):894–909, 2012
2012
-
[10]
F.-H. Lin, T. Witzel, S. P. Ahlfors, S. M. Stufflebeam, J. W. Belliveau, and M. S. H¨ am¨ al¨ ainen. Assessing and improving the spatial accuracy in meg source localization by depth-weighted minimum-norm estimates. NeuroImage, 31(1):160–171, 2006
2006
-
[11]
J. S. Liu and R. Chen. Sequential monte carlo methods for dynamic systems. Journal of the American Statistical Association , 93:1032–1044, 1998
1998
-
[12]
Liu, J.-H
S. Liu, J.-H. Poh, H. L. Koh, K. K. Ng, Y. M. Loke, J. K. W. Lim, J. S. X. Chong, and J. Zhou. Carrying the past to the future: Distinct brain networks underlie individual differences in human spatial working memory capacity. NeuroImage, 176:1–10, 2018
2018
-
[13]
C. J. Long, P. L. Purdon, S. Temereanca, N. U. Desai, M. S. H¨ am¨ al¨ ainen, and E. N. Brown. State-space solutions to the dynamic magnetoencephalography inverse prob- lem using high performance computing. The Annals of Applied Statistics, 5(2B):1207– 1228, 2011
2011
-
[14]
J. C. Mosher, P. S. Lewis, and R. M. Leahy. Multiple dipole modeling and localiza- tion from spatio-temporal meg data. IEEE Transactions on Biomedical Engineering, 39(6):541–557, 1992
1992
-
[15]
W. Ou, M. S. H¨ am¨ al¨ ainen, and P. Golland. A distributed spatio-temporal eeg/meg inverse solver. NeuroImage, 44(3):932–946, 2009
2009
-
[16]
R. D. Pascual-Marqui, C. M. Michel, and D. Lehmann. Low resolution electro- magnetic tomography: a new method for localizing electrical activity in the brain. International Journal of Psychophysiology , 18(1):49–65, 1994
1994
-
[17]
L. R. Rabiner. A tutorial on hidden markov models and selected applications in speech recognition. Proceedings of the IEEE, 77(2):257–286, 1989
1989
-
[18]
J. Sarvas. Basic mathematical and electromagnetic concepts of the biomagnetic in- verse problem. Physics in Medicine and Biology , 32:11–12, 1984
1984
-
[19]
N. J. Trujillo-Barreto, E. Aubert-V´ azquez, and W. D. Penny. Bayesian m/eeg source reconstruction with spatio-temporal priors. NeuroImage, 39(1):318–335, 2008. 38
2008
-
[20]
Uutela, M
K. Uutela, M. H¨ am¨ al¨ ainen, and E. Somersalo. Visualization of magnetoencephalo- graphic data using minimum current estimates. NeuroImage, 10(2):173–180, 1999
1999
-
[21]
F. Vaida. Parameter convergence for em and mm algorithms. Statistica Sinica , 15(3):831–840, 2005
2005
-
[22]
B. Veen, J. Joseph, and K. Hecox. Localization of intra-cerebral sources of electrical activity via linearly constrained minimum variance spatial filtering. In Proceedings of IEEE Workshop on Statistical Signal and Array Processing , 1:526–529, 1992
1992
-
[23]
Yao and W
Z. Yao and W. Eddy. A statistical approach to the inverse problem in magnetoen- cephalography. Annals of Applied Statistics , 8:1119–1144, 2014
2014
-
[24]
Z. Yao, Y. Zhang, Z. Bai, and W. F. Eddy. Estimating the number of sources in mag- netoencephalography using spiked population eigenvalues. Journal of the American Statistical Association, 113(522):505–518, 2018
2018
-
[25]
Zhang and C
J. Zhang and C. Liu. On linearly constrained minimum variance beamforming. The Journal of Machine Learning Research , 16(1):2099–2145, 2015
2015
-
[26]
Zhang and L
J. Zhang and L. Su. Temporal autocorrelation-based beamforming with meg neu- roimaging data. Journal of the American Statistical Association, 110(512):1375–1388, 2015. Department of Statistics and Applied Probability National University of Singapore, 21 Lower Kent Ridge Road, ...
2015
Reviewed August 14, 2026 · model on record in the stance chip above.
Discussion (0). Continue with ORCID to comment.