Multilevel Monte Carlo for a class of partially observed processes in neuroscience

Mohamed Maama; Ajay Jasra; Kengo Kamatani

doi:10.3934/fods.2025012

Article Contents

2026, Volume 9: 61-87. Doi: 10.3934/fods.2025012

This volume Previous Article The persistent Laplacian of non-branching complexes Next Article Convergence rates of non-stationary and deep Gaussian process regression

Multilevel Monte Carlo for a class of partially observed processes in neuroscience

1.
Applied Mathematics and Computational Science Program, KAUST, Kingdom of Saudi Arabia
2.
School of Data Science, The Chinese University of Hong Kong, Shenzhen, China
3.
Institute of Statistical Mathematics, Tokyo 190-0014, Japan

^*Corresponding author: Mohamed Maama
^*Corresponding author: Mohamed Maama

Received: November 2024

Revised: July 2025

Early access: July 17, 2025

Published online: July 17, 2025

Abstract / Introduction Full Text(HTML) Figure(14) / Table(2) Related Papers Cited by

Abstract

In this paper, we consider Bayesian parameter inference associated with a class of partially observed stochastic differential equations (SDE) driven by jump processes. Such type of models can be routinely found in applications, for which we focus upon the case of neuroscience. The data are assumed to be observed regularly in time and driven by the SDE model with unknown parameters. In practice, the SDE may not have an analytically tractable solution, and this leads naturally to a time discretization. We adapt the multilevel Markov chain Monte Carlo method of [18], which works with a hierarchy of time discretizations, and show empirically and theoretically that this is preferable to using one single time discretization. The improvement is in terms of the computational cost needed to obtain a pre-specified numerical error. Our approach is illustrated on models applied to examples of Bayesian inference problems using both simulated and real observed data from neuroscience.

Keywords:

Mathematics Subject Classification: Primary: 65C05, 60H35; Secondary: 65C30.

Citation:

Full Text(HTML)

Figure 1. Excitatory-inhibitory neuronal networks scheme

Download: Full-size image PowerPoint slide

Figure 2. Raster plots of firing-activity for two different systems during $ 400 \; ms $, inhibitory (blue - bottom half) and excitatory (red - top half) integrate-and-fire neurons showing two different regimes. Top: A homogeneous property with minimal to no correlations among spike times, $ S^{\text{EE}} = S^{\text{II}} = S^{\text{EI}} = S^{\text{IE}} = 0.0028 $. Bottom: A more synchronous regime with the majority of the network firing in perfect synchronization, $ S^{\text{EE}} = S^{\text{II}} = 0.009 $, $ S^{\text{EI}} = S^{\text{IE}} = 0.007 $. Each regime has network parameters, $ N_{\text{E}} = N_{\text{I}} = 200 $, and the Poisson random inputs are the amplitude $ S^{\text{dr}} = 0.065 $, the frequency per ms $ \lambda = 0.55 $

Download: Full-size image PowerPoint slide

Figure 3. PMMH outputs for our first model. (A) ACF plot for $ S^{\text{dr}} $; (B) Trace plot for estimated $ S^{\text{dr}} $, with the average value indicated by a black line; (C) Posterior distribution histogram of the estimated parameter of equation (22)

Download: Full-size image PowerPoint slide

Figure 4. Mean-squared error versus cost function for the parameter of interest $ S^{\text{dr}} $ in equation (22)

Download: Full-size image PowerPoint slide

Figure 5. Simulation output of equation (25) : all plots are for $ S^{IE} $. (A) Trace plot of PMCMC chains for $ S^{IE} $; (B) Posterior function plot; (C) MSE vs. Cost plot

Download: Full-size image PowerPoint slide

Figure 6. Simulation outputs of equation (25) : all plots are for $ S^{EI} $. (A) Trace plot of PMCMC chains for $ S^{EI} $; (B) Posterior function plot; (C) Mean-squared error vs. Cost function plot

Download: Full-size image PowerPoint slide

Figure 7. Membrane Potential vs. Time and Power Spectral Density for $ S^{dr} = 0.1 $

Download: Full-size image PowerPoint slide

Figure 8. Membrane Potential vs. Time and Power Spectral Density for $S^{dr} = 0.5 $

Download: Full-size image PowerPoint slide

Figure 9. Membrane Potential vs. Time and Power Spectral Density for $ S^{dr} = 1 $

Download: Full-size image PowerPoint slide

Figure 10. Posterior and Prior Distributions for Parameter $a$ derived from Resting and Stimulated Real Data. The figure compares the posterior distributions (red) with the prior distributions (blue) for the parameter $a$

Download: Full-size image PowerPoint slide

Figure 11. Posterior and Prior Distributions for Parameter $b$ derived from Resting and Stimulated Real Data. The figure highlights the shift from the prior to the posterior for the parameter $b$

Download: Full-size image PowerPoint slide

Figure 12. Posterior and Prior Distributions for Parameter $c$ derived from Resting and Stimulated Real Data. The posterior distributions provide insights into the inferred reset membrane potential values for parameter $c$

Download: Full-size image PowerPoint slide

Figure 13. Posterior and Prior Distributions for Parameter $d$ derived from Resting and Stimulated Real Data. The comparison indicates the adjustments in the reset recovery variable $d$

Download: Full-size image PowerPoint slide

Figure 14. Posterior and Prior Distributions for Parameter $S^{dr}$ derived from Resting and Stimulated Real Data. The posterior distribution reflects the estimated noise strength $S^{dr}$

Download: Full-size image PowerPoint slide

Table 1. Estimated rates of convergence of mean-squared error with respect to the cost function for the three key parameters $ \big( S^{\text{dr}}, \; S^{\text{EI}}, \; S^{\text{IE}} \big) $, adapted to the curves simulated above. MLPMCMC is a multilevel PMCMC

Model	Parameter	PMCMC	MLPMCMC
Case 1: Integrate-and-Fire Model with Poisson Input	$ S^{\text{dr}} $	-1.54	-1.02
Case 2: Small Network of E-I Coupled Neurons	$ S^{\text{EI}} $	-1.47	-1.05
Case 2: Small Network of E-I Coupled Neurons	$ S^{\text{IE}} $	-1.52	-1.02

| Show Table

DownLoad: CSV

Table 2. Summary of priors for parameter inference under resting and stimulated cases

Parameter	Resting (1551, 1608)	Stimulated (1554, 1612)
$a$ (Recovery time)	$\log\mathcal{N}(-3, 0.5^2)$	$\log\mathcal{N}(-2.5, 0.6^2)$
$b$ (Recovery sensitivity)	$\mathcal{U}(0.1, 0.5)$	$\mathcal{U}(0.2, 0.6)$
$c$ (Reset membrane potential)	$\mathcal{N}(-65, 5^2)$	$\mathcal{N}(-60, 5^2)$
$d$ (Reset recovery variable)	$\mathcal{N}(8, 1^2)$	$\mathcal{N}(9, 1^2)$
$S^{\text{dr}}$ (Poisson noise strength)	$\Gamma(2, 0.5)$	$\Gamma(3, 0.7)$

| Show Table

DownLoad: CSV

Related Papers

Cited by

References

[1]	C. Andrieu, A. Doucet and R. Holenstein, Particle Markov chain Monte Carlo methods (with discussion), J. R. Statist. Soc. Ser. B, 72 (2010), 269-342. doi: 10.1111/j.1467-9868.2009.00736.x.
[2]	J.-M. Bernardo and A. F. M. Smith, Bayesian Theory, Wiley, London, 1994. doi: 10.1002/9780470316870.
[3]	N. Bruti-Liberati and E. Platen, Strong approximations of stochastic differential equations with jumps, J. Comput. Appl. Math., 205 (2007), 982-1001. doi: 10.1016/j.cam.2006.03.040.
[4]	A. N. Burkitt, A review of the integrate-and-fire neuron model: I. Homogeneous synaptic input, Biol. Cybern., 95 (2006), 1-19. doi: 10.1007/s00422-006-0068-6.
[5]	N. K. Chada, J. Franks, A. Jasra, K. J. Law and M. Vihola, Unbiased inference for discretely observed hidden Markov model diffusions, SIAM/ASA J. Uncertain. Quantif., 9 (2021), 763-787. doi: 10.1137/20M131549X.
[6]	L. Chariker, R. Shapley and L.-S. Young, Rhythm and synchrony in a cortical network model, J. Neurosci., 38 (2018), 8621-8634. doi: 10.1523/JNEUROSCI.0675-18.2018.
[7]	L. Chariker and L.-S. Young, Emergent spike patterns in neuronal populations, J. Comput. Neurosci., 38 (2015), 203-220. doi: 10.1007/s10827-014-0534-4.
[8]	D. J. Daley and D. Vere-Jones, An Introduction to the Theory of Point Processes, Volume I: Elementary Theory and Methods, Springer, New York, 2003.
[9]	P. Del Moral, Feynman-Kac Formulae: Genealogical and Interacting Particle Systems with Applications, Springer, New York, 2004. doi: 10.1007/978-1-4684-9393-1.
[10]	S. Ditlevsen and A. Samson, Estimation in the partially observed Morris-Lecar neuronal model with particle filter and stochastic approximation methods, Ann. Appl. Stat., 8 (2014), 674-702. doi: 10.1214/14-AOAS729.
[11]	T. J. Dodwell, C. Ketelsen, R. Scheichl and A. L. Teckentrup, A hierarchical multilevel Markov chain Monte Carlo algorithm with applications to uncertainty quantification in subsurface flow, SIAM/ASA J. Uncertain. Quantif., 3 (2015), 1075-1108. doi: 10.1137/130915005.
[12]	W. Gerstner and W. M. Kistler, Spiking Neuron Models, Cambridge University Press, Cambridge, 2002. doi: 10.1017/CBO9780511815706.
[13]	M. B. Giles, Multilevel Monte Carlo path simulation, Oper. Res., 56 (2008), 607-617. doi: 10.1287/opre.1070.0496.
[14]	S. Heinrich, Multilevel Monte Carlo methods, in Large-Scale Scientific Computing, Springer, Berlin, 2001, 58-67. doi: 10.1007/3-540-45346-6_5.
[15]	V. H. Hoang, C. Schwab and A. M. Stuart, Complexity analysis of accelerated MCMC methods for Bayesian inversion, Inverse Probl., 29 (2013), 085010, 37 pp. doi: 10.1088/0266-5611/29/8/085010.
[16]	E. M. Izhikevich, Simple model of spiking neurons, IEEE Trans. Neural Netw., 14 (2003), 1569-1572. doi: 10.1109/TNN.2003.820440.
[17]	A. Jasra, J. Heng, Y. Xu and A. N. Bishop, A multilevel approach for stochastic nonlinear optimal control, Int. J. Control, 95 (2022), 1290-1304. doi: 10.1080/00207179.2020.1849805.
[18]	A. Jasra, K. Kamatani, K. Law and Y. Zhou, Bayesian static parameter estimation for partially observed diffusions via multilevel Monte Carlo, SIAM J. Sci. Comput., 40 (2018), A887-A902. doi: 10.1137/17M1112595.
[19]	A. Jasra, K. Law and C. Suciu, Advanced multilevel Monte Carlob methods, Int. Stat. Rev., 88 (2020), 548-579. doi: 10.1111/insr.12365.
[20]	A. Jasra, M. Maama and A. Mijatović, Modeling of measurement error in financial returns data, preprint, 2024, arXiv: 2408.07405.
[21]	C. Koch, Biophysics of Computation: Information Processing in Single Neurons, Oxford University Press, Oxford, 2004.
[22]	S. Koyama and L. Paninski, Efficient computation of the maximum a posteriori path and parameter estimation in integrate-and-fire and more general state-space models, J. Comput. Neurosci., 29 (2010), 89-105. doi: 10.1007/s10827-009-0150-x.
[23]	K. W. Latimer, E. J. Chichilnisky, F. Rieke and J. W. Pillow, Inferring synaptic conductances from spike trains with a biophysically inspired point process model, Adv. Neural Inf. Process. Syst., 27 (2014).
[24]	M. Maama, B. Ambrosio, M. A. Aziz-Alaoui and S. M. Mintchev, Emergent properties in a V1-inspired network of Hodgkin-Huxley neurons, Math. Model. Nat. Phenom., 19 (2024), Paper No. 3, 26 pp. doi: 10.1051/mmnp/2024001.
[25]	M. Maama, A. Jasra and H. Ombao, Bayesian parameter inference for partially observed SDEs driven by fractional Brownian motion, Stat. Comput., 33 (2023), Paper No. 19, 9 pp. doi: 10.1007/s11222-022-10193-0.
[26]	L. Meng, M. A. Kramer and U. T. Eden, A sequential Monte Carlo approach to estimate biophysical neural models from spikes, J. Neural Eng., 8 (2011), 065006. doi: 10.1088/1741-2560/8/6/065006.
[27]	B. Metcalfe, A. Hunter, J. Graham-Harper-Cater and J. Taylor, A dataset of action potentials recorded from the L5 dorsal rootlet of rat using a multiple electrode array, Data Brief, 33 (2020), 106561. doi: 10.1016/j.dib.2020.106561.
[28]	M. J. Moye and C. O. Diekman, Data assimilation methods for neuronal state and parameter estimation, J. Math. Neurosci., 8 (2018), 1-38. doi: 10.1186/s13408-018-0066-8.
[29]	L. Paninski, J. Pillow and E. Simoncelli, Comparing integrate-and-fire models estimated using intracellular and extracellular data, Neurocomputing, 65 (2005), 379-385. doi: 10.1016/j.neucom.2004.10.032.
[30]	A. V. Rangan and D. Cai, Fast numerical methods for simulating large-scale integrate-and-fire neuronal networks, J. Comput. Neurosci., 22 (2007), 81-100. doi: 10.1007/s10827-006-8526-7.
[31]	A. V. Rangan and L.-S. Young, Emergent dynamics in a model of visual cortex, J. Comput. Neurosci., 35 (2013), 155-167. doi: 10.1007/s10827-013-0445-9.
[32]	A. V. Rangan and L.-S. Young, Dynamics of spiking neurons: Between homogeneity and synchrony, J. Comput. Neurosci., 34 (2013), 433-460. doi: 10.1007/s10827-012-0429-1.
[33]	C. P. Robert and G. Casella, Monte Carlo Statistical Methods, Springer, New York, 2004. doi: 10.1007/978-1-4757-4145-2.
[34]	J. W. Zhang and A. V. Rangan, A reduction for spiking integrate-and-fire network dynamics ranging from homogeneity to synchrony, J. Comput. Neurosci., 38 (2015), 355-404. doi: 10.1007/s10827-014-0543-3.