Pith. sign in

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 →

arxiv 1908.03926 v1 pith:2LRYHMWU submitted 2019-08-11 stat.AP stat.ME

classification stat.APstat.ME MSC 62F1562M0562P10
keywords MEGinverseproblemdiscreteposteriordistributionExpectation-Maximizationspatio-temporalmodelsourcelocalizationhiddenMarkovforward-backwardalgorithmmagnetoencephalography
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper aims to quantify where a brain source is located over time in magnetoencephalography (MEG), a problem that is ill-posed and nonlinear in the location parameter. Its central claim is that discretizing the region of interest into voxels and treating the source location as a hidden Markov state yields a closed-form discrete posterior distribution at every time point, computed from the entire sequence of measurements rather than only the past. This avoids the slow convergence of sequential sampling methods while still describing the time evolution of the source. If the claim holds, practitioners get a tractable, probabilistic map of source activity across the whole scan.

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.

Watch

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

Editorial extensions of the paper, not claims the author makes directly.

  • 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.
Share X Bluesky LinkedIn Reddit HN

Signed reviews

No signed human review yet.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

5 major / 6 minor

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)
  1. [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.
  2. [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.
  3. [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.
  4. [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.
  5. [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)
  1. [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.
  2. [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.
  3. [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.
  4. [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.
  5. [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.
  6. [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

0 steps flagged · score 0.0 of 10

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 4 free parameters · 7 assumptions · 0 invented entities

The central claim rests on a standard HMM framework plus several domain-specific modeling choices. The most consequential inputs pulled from outside the derivation are the fixed dipole moment, the pre-specified ROI and mesh, the hand-set noise and initial-state parameters in the real data, and the heuristic shrink factor. The model itself also assumes Gaussian AR dynamics and uncorrelated Gaussian sensor noise, both of which are known simplifications. No new physical or mathematical entities are introduced.

free parameters (4)
  • Dipole moment vector q = q = (3, 3, 3) in simulations; unit moment in real data
    The moment (orientation and strength) is fixed by hand rather than estimated or updated; the location posterior strongly depends on this value (Eqs. 2.1, 2.3).
  • Initial ROI and mesh grid sizes = K_i = 10 in simulations; ROI around MNE for real data
    The discrete posterior and the EM estimates depend on the chosen initial ROI and mesh resolution; the paper provides no sensitivity analysis.
  • Manually set mu_0, Sigma_0, Sigma in real data = Chosen from MNE and Fourier noise estimates
    In Section 4.1, the source model parameters mu_0, Sigma_0, and Sigma are set by hand, which constrains the source dynamics before EM.
  • ROI shrinking factor (3-sigma) = 3
    The dynamic procedure shrinks the ROI to the posterior mean plus/minus three posterior standard deviations; this constant is chosen without justification.
assumptions (7)
  • domain assumption AR(1) Gaussian process for the source state
    Eq. (2.4) assumes the six-dimensional source state (location and moment) follows a first-order auto-regression with Gaussian innovations; the posterior is misspecified if the true dynamics are nonlinear or non-Gaussian.
  • domain assumption Gaussian, homogeneous, uncorrelated sensor noise
    Section 2.1 assumes Ut ~ N(0, sigma^2 I), uncorrelated across sensors; real MEG sensor noise is spatially correlated.
  • domain assumption Quasi-static approximation of Maxwell's equations and the Biot-Savart forward model
    Eq. (2.1) is the standard quasi-static MEG forward model (citing [18]); deviations from this model are not considered.
  • ad hoc to paper Discretization approximates the continuous posterior
    Eqs. (2.6)-(2.8) replace continuous transition and emission densities with voxel-level probabilities evaluated at voxel centers; no error bound is provided for this approximation.
  • standard math Regularity conditions (C1)-(C7)
    Theorem 2.1 relies on these conditions for EM convergence; they are standard in the EM literature but not verified for this model.
  • domain assumption A true ROI exists and contains the source for all time points
    Section 2.5 assumes the source is restricted to ROI0; if a source leaves the ROI, the posterior cannot represent it.
  • ad hoc to paper Switch procedure fixed point approximates the joint posterior marginals
    Eq. (2.17) and Theorem 2.2 claim the switch marginal approximates the true joint marginal, but the theorem's bound is vacuous ('there exists epsilon > 0' with no practical size).

how reviews work

0 comments
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.

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

26 extracted references · 26 canonical work pages

  1. [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

  2. [2]

    Baillet and L

    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

  3. [3]

    Baillet, J

    S. Baillet, J. C. Mosher, and R. M. Leahy. Electromagnetic brain mapping. IEEE Signal Processing Magazine, 18(6):14–30, 2001

  4. [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

  5. [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

  6. [6]

    Fukushima, O

    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

  7. [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

  8. [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

Show all 26 references
  1. [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

  2. [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

  3. [11]

    J. S. Liu and R. Chen. Sequential monte carlo methods for dynamic systems. Journal of the American Statistical Association , 93:1032–1044, 1998

  4. [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

  5. [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

  6. [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

  7. [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

  8. [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

  9. [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

  10. [18]

    J. Sarvas. Basic mathematical and electromagnetic concepts of the biomagnetic in- verse problem. Physics in Medicine and Biology , 32:11–12, 1984

  11. [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

  12. [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

  13. [21]

    F. Vaida. Parameter convergence for em and mm algorithms. Statistica Sinica , 15(3):831–840, 2005

  14. [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

  15. [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

  16. [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

  17. [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

  18. [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, ...

Pith tools

Reviewed August 14, 2026 · model on record in the stance chip above.