Inference on the Macroscopic Dynamics of Spiking Neurons Open Access

Neural computation involves the complex information processing performed by spiking neurons, where their interactions within neural circuits contribute to higher-level computations and cognitive functions. Consequently, the cornerstone of theoretical neuroscience lies in the use of mathematical models to understand the intricacies of neural computation (Dayan & Abbott, 2005; Hertz, 2018). This approach enables researchers to explore the complexities of neural processes, thereby revealing insights into the mechanisms that drive brain function and behavior. These mathematical models can range from simple idealized representations of single neurons (Hopfield, 1982; Izhikevich, 2003) to complex network models that simulate the interactions within neural circuits (Marder, 1998; Sussillo, 2014; O’Leary et al., 2015; Bittner et al., 2021). The biological neural computation forms the basis for the brain’s ability to execute various inference tasks, such as recognizing patterns, making decisions, and generating responses to stimuli. This capability relies on the collective behavior of spiking neurons, where the interactions and coordination among these neurons within neural circuits enable the brain to process, integrate, and interpret sensory information, ultimately leading to higher-level cognitive functions and adaptive behavior (Kandel et al., 2000; Gerstner et al., 2014; Friston & Kiebel, 2009; Friston et al., 2017).

Mean-field (MF) models serve as effective computational abstractions that represent the collective behavior of large populations of neurons while maintaining a degree of mathematical tractability (Amari, 1977; Wilson & Cowan, 1973; Jirsa & Haken, 1996; David & Friston, 2003; Deco et al., 2008; Hutt et al., 2015; Coombes & Byrne, 2018; Bandyopadhyay et al., 2021; Cook et al., 2022). Hence, MF models facilitate the study of information processing and computation within the brain, which underlie cognitive processes such as perception, memory, and learning. Nevertheless, deriving and validating MF models from spiking neural networks present many challenges, particularly when assessing how the interaction between individual neurons and parameters leads to macroscopic behavior that aligns with the averaged activity of neural populations.

Recently, an analytically driven MF model of spiking neurons has been formulated that can effectively describe all potential macroscopic dynamical states of the network, including states of synchronous spiking activity (Montbrió et al., 2015). However, the operation of such complex systems (governed by coupled nonlinear differential equations) is determined by the selection of (biological or phenomenological) parameters that, when set in a specific configuration, give rise to a measurable signature of a computation (Achard & De Schutter, 2006; Sussillo, 2014). Analyzing MF models and comparing their emergent dynamics against a network of neurons involves solving inverse problems to ascertain the optimal parameter setting. This process requires swift and robust outcomes to inform real-time decisions and also to deal with observation and dynamical noise. Yet even in the simplest models, there can be a degenerate relationship between the model parameters and its overall emergent function (Edelman & Gally, 2001; Prinz et al., 2004; Alonso & Marder, 2019), making the inverse problem more challenging.

Maintaining the interdependency between parameters by traversing across scales adds a layer of complexity to the validation process. It is crucial to distinguish between genuine, biologically relevant correlations and artificial correlations that may arise from the inference process or modeling assumptions. Furthermore, due to the intricate nature of computation within neural circuits, it becomes intractable to analytically derive MF models that include more biological realism (such as adaption, neuromodulation, extra-synaptic transmission, and E/I ratios). Statistical inference offers an efficient and adaptable approach to solving the inverse problem by identifying approximate parameter distributions that are responsible for generating computations in a biologically realistic model (Achard & De Schutter, 2006; Liepe et al., 2014; Lueckmann et al., 2017; Gonçalves et al., 2020; Bittner et al., 2021; Młynarski et al., 2021).

The performance of statistical inference algorithms depends on the task, and there is no universally best algorithm for different inverse problems. Therefore, we conducted a benchmarking analysis against state-of-the-art optimization and Bayesian estimation algorithms to discern their respective advantages and limitations. In practice, optimization methods are commonly used to quickly determine unknown quantities through a single point estimate (Mendes & Kell, 1998; Nocedal & Wright, 1999; Kelley, 1999; Floudas & Gounaris, 2009). These methods involve iteratively adjusting parameters to minimize or maximize an objective function, scoring the model’s performance against observed data—for example, through minimizing distance errors or maximizing correlation (Banga & Balsa-Canto, 2008; Tashkova et al., 2011; Svensson et al., 2012; Hashemi et al., 2018).

The Bayesian approach offers a principled method for making inferences and predictions, establishing relationships between parameters, and quantifying uncertainty in the decision-making process (Gelman et al., 1995, 2020; Bishop, 2006; van de Schoot et al., 2021). In Bayesian modeling, all model parameters are treated as random variables, and their values are subject to variation based on their underlying probability distributions. Such probabilistic techniques provide the full posterior distribution of unknown quantities hidden in the underlying data-generating process. The uncertainty and interdependency in Bayesian estimation are naturally quantified by assigning a probability distribution to each parameter (known as the prior distribution), which is then updated based on the information provided by the data (referred to as the likelihood function). To conduct a fully Bayesian procedure, the state-of-the-art MCMC method is adaptive Hamiltonian Monte Carlo (HMC; Duane et al., 1987; Neal, 2010; Hoffman & Gelman, 2014), which utilizes gradient information to avoid random walk behavior. This enables efficient sampling from high-dimensional distributions that may exhibit strong correlations (Betancourt, 2017).

Simulation-based inference (SBI; Cranmer et al., 2020; Brehmer, 2021) or likelihood-free inference (Papamakarios & Murray, 2016; Brehmer et al., 2020) leverages deep generative models to conduct approximate Bayesian estimation, using low-dimensional data features that are generated by random simulations (Gonçalves et al., 2020; Lueckmann et al., 2021; Boelts et al., 2022). In this efficient approach, a simple base probability distribution (prior) is transformed into a more complex distribution (posterior) through a sequence of invertible transformations (i.e., normalizing flows; Rezende & Mohamed, 2015; Papamakarios, Nalisnick, et al., 2019). Notably, it allows for direct estimation of joint posterior distributions, bypassing the need for MCMC sampling (Greenberg et al., 2019; Papamakarios, Sterratt, et al., 2019). Moreover, expressive deep generative models have the potential to capture parameter nonlinear relationships between parameters and multimodalities in the distributions (Hashemi et al., 2023).

Data-driven methods for learning dynamical models from time-series data have been extensively researched for several decades (Juang, 1994; Ljung, 1998; Brunton et al., 2016; Linderman et al., 2017; Duncker et al., 2019; Koppe et al., 2019; Sip et al., 2023). Instead of relying on discretized maps, neural ordinary differential equations (neural ODEs; Chen et al., 2018) form a new family of deep neural network models for modeling continuous-time dynamics. Neural ODEs define the vector fields and ODE solution as a black-box differential equation solver, allowing for uncovering the dynamics of a system even when the governing equations are unknown (Dupont et al., 2019; Biloš et al., 2021). This data-driven approach involves parameterizing system dynamics as continuous functions, enabling smooth and uninterrupted modeling of temporal evolution (Yan et al., 2019; Kim et al., 2021). Neural ODEs naturally adapt to varying time intervals and can accommodate fluctuations in the frequency of data observations (Zhu et al., 2022; Goyal & Benner, 2023).

Through an exploration of the aforementioned methods, we demonstrate that global optimization algorithms, such as differential evolution algorithm, offer fast and accurate point estimation of the true generative parameters when the dynamical noise is absent. However, when dealing with dynamic evolution subject to noise, SBI, using deep neural density estimators, emerges as the superior approach, outperforming other algorithms, such as adaptive HMC sampling. Additionally, when dealing with missing data (such as population firing rate) in state-space modeling, HMC fails to capture the dynamics of bistable switching behavior. Instead, SBI is able to accurately recover the diverse dynamics in the phase-space representation. Nevertheless, this approach may lead to an overestimation of uncertainty and correlation between parameters, which can be mitigated by using a time-delay embedding technique to improve the results. We validate this by employing an MF model of quadratic integrate-and-fire (QIF) neurons, which demonstrates that the interdependencies of parameters are maintained when traversing across scales. Finally, we demonstrate that training deep neural ODEs on spiking neurons enables the inference of vector fields at macroscopic level. This allows for the prediction of emergent behaviors and system dynamics based on the microscopic state of the spiking neurons. Note that we use the following abbreviations throughout this article: MF, mean-field; QIF, quadratic Integrate-and-fire; SBI, simulation-based inference; SNPE, sequential neural posterior estimation; MCMC, Markov chain Monte Carlo; HMC, Hamiltonian Monte Carlo; MAP, maximum a posteriori; DE, differential evolution; PSO, particle swarm optimization; BO, Bayesian optimization; ODEs, ordinary differential equations; SDEs, stochastic differential equations; neural ODEs, neural ordinary differential equations.

The quadratic integrate-and-fire (QIF) neurons are a class of simplified computational models that are extensively used to study the dynamics of spiking neurons (Gerstner & Kistler, 2002; Izhikevich, 2007). In the QIF model, the membrane potential of each neuron evolves according to a quadratic differential equation until it reaches a threshold, at which point the neuron emits a spike and the potential is reset.

We validate the MF approximation by comparing it against detailed spiking neurons using various parameter estimation inference methods. To evaluate these methods, we use synthetic data generated from the MF model and a network of QIF neurons. Inference methods investigated in this work can be broadly divided into two classes, optimization and Bayesian methods.

From the optimization methods, we use the global search algorithms that incorporate a bio-inspired random search principle: differential evolution (DE; Storn & Price, 1997; Price, 1999), and particle swarm optimization (PSO; Kennedy & Eberhart, 1995; Eberhart & Kennedy, 1995). These algorithms do not require an initial guess for the parameters or the gradient information of the objective function. We also consider Bayesian optimization (BO; Snoek et al., 2012; Shahriari et al., 2015), an algorithm that constructs a probabilistic model for the objective function using gaussian processes. This approach allows for the integration of uncertainty in the optimization process.

From the Bayesian methods, we compare the results of two state-of-the-art Bayesian computation algorithms: Hamiltonian Monte Carlo (HMC; Duane et al., 1987; Neal, 2010) which is unbiased and exact in infinite runs, and a simulation-based inference (SBI; Cranmer et al., 2020; Brehmer, 2021), which approximates (or parameterizes) the posterior using deep generative models. Here, generative modeling is an unsupervised machine learning method to model a probability distribution based on the samples drawn from that distribution. From the Bayesian methods, we also report the results of maximum a posteriori (MAP) estimation.

Markov chain Monte Carlo (MCMC) is a powerful class of computational algorithms used for sampling from a distribution, in which the sampling process does not require knowledge of the entire distribution, making it a versatile tool (Andrieu et al., 2003; Murphy, 2022; McElreath, 2020). MCMC is unbiased and asymptotically exact in the limit of infinite runs. Hamiltonian Monte Carlo (HMC; Duane et al., 1987; Neal, 2010) is a gradient-based MCMC designed to avoid random walk behavior, and it can efficiently sample from high-dimensional distributions that may exhibit strong correlations (Betancourt, 2017). However, the efficiency of HMC is sensitive to the algorithm parameters.

In this study we use a self-tuning variant of HMC, known as the No-U-Turn Sampler (NUTS; Hoffman & Gelman, 2014) from a high-level statistical modeling tool called Stan (Carpenter et al., 2017). In particular, the NUTS calibrates the number of steps and step size of the leapfrog integrator (in solving the Hamiltonian equations of motion) during a warm-up phase to achieve a target Metropolis acceptance rate. For more details, see Betancourt (2013) and Baldy et al. (2023). Moreover, Stan offers alternative methods such as MAP estimation using gradient-based optimization, automatic differentiation for efficient gradient computation, and various diagnostics to assess the convergence of the inference process (see https://mc-stan.org).

Simulation-based inference (SBI) conducts efficient Bayesian inference for complex models when the calculation of the likelihood function is either analytically or computationally intractable (Cranmer et al., 2020; Brehmer, 2021). In computational models, where the data can be generated through stochastic simulations, SBI leverages repeated simulations from the generative model and employs probabilistic machine learning to estimate a target probability distribution. Instead of directly sampling from distributions using MCMC or explicitly evaluating the likelihood function, SBI overcomes these challenges by using deep neural density estimators, such as masked autoregressive flows (MAF; Papamakarios et al., 2017). These density estimators learn an invertible transformation between the distributions of parameters and low-dimensional data features at a very low computational cost to efficiently sample from distributions.

Taking the prior distribution $p(θ)$ over the parameters $θ$⁠, a limited number of $N$ simulations are generated for a training step from the generative model: ${(θi,xi)}i=1N∼p(θ,x)$⁠. After the training step, we are able to quickly estimate the approximated posterior $qφ(θ∣x)$ with learnable parameters $φ$⁠, so that for the observed data, $xobs$⁠: $qφ(θ∣xobs)≃p(θ∣xobs)$⁠. (For details, see Gonçalves et al., 2020; Hashemi et al., 2023.) Minimizing the Kullback-Leibler divergence between the parameterized (approximate) posterior and the true posterior is also the objective in dynamical causal modeling (Friston et al., 2003; Blei et al., 2017), which is based on maximizing a lower bound on the marginal likelihood of the data, or equivalently, minimizing the free energy.

The methods for SBI often include a sequential training procedure, which adaptively guides simulations to yield more informative estimates (Papamakarios, Sterratt, et al., 2019; Lueckmann et al., 2019; Durkan et al., 2020; Wiqvist et al., 2021; Deistler et al., 2022). In particular, sequential neural posterior estimation (SNPE; Greenberg et al., 2019; Gonçalves et al., 2020) dynamically refines the proposals, network weights, and posterior estimates to learn the relationships between model parameters and the observed summary statistics of the data. In this study, we used SNPE with a single round to take advantage of an amortized strategy; After incurring an initial computational cost for the simulation and training steps to learn all the joint posterior distributions, the posterior can be quickly estimated from any new observations (by a forward pass through neural networks) without any additional computational overhead or further simulations.

Time-delay embedding is a commonly used technique for characterizing dynamical systems based on limited measurements, time-series analysis, and prediction (Takens, 1981). In time-delay embedding, the reconstruction of a latent high-dimensional system relies on incorporating incomplete measurements, along with a temporal history of preceding measurements to create a comprehensive representation (Kennel et al., 1992; Hirsh et al., 2021). In a subsequent analysis, we challenge the inference process by assuming that the firing rate activity $r$ is not directly observed (missing data problem). Instead, to improve the inference, we recovered the latent time series $rrec$ from the observed activity $v$⁠, which is coupled to $r$ according to equation 2.1. By expanding on a method introduced by Abarbanel et al. (1994), which primarily focuses on predicting physical variables in time-delay embedding, we removed the assumption that we have the access to training data points from the hidden time series. Instead, we leverage our understanding of the generator to simulate data pairs $(r,v)$ that can be used for training purposes.

Time-delay embedding requires the setting of hyperparameters, such as the delay (time lag), $T$⁠, and the dimension of the embedding space, $d$⁠. These hyperparameters are typically set prior to training, often based on the minimization of mutual information to determine the appropriate delay and the false nearest neighbors method to determine the number of embedding dimensions (Kennel et al., 1992; Tan et al., 2023). However, in the our application here, we have found that the hyperparameters suggested by these methods were not the most effective in achieving an accurate fit. Instead, we select hyperparameters that minimize the mean square error of the fit to simulated data, with $T=160$ points (i.e., $0.16sec$⁠) and $d=12$⁠. To do this, a set of 100 pairs of coupled time series $(r,v)$ was simulated, with an observation noise intensity of 0.1 and varying parameters (⁠$Δ,η,J$⁠). These pairs were then used to infer regression coefficients that closely match $r$ to the delay embedding space representation of $v$⁠.

We use Python implementations of the optimization and Bayesian algorithms: DE from SciPy (Virtanen et al., 2020), PSO from the toolkit PySwarms (Miranda, 2018), and BO from the package Bayesian optimization (Nogueira et al., 2014). DE is set with a population size of 10 and the maximum iterations of 500. PSO is set with 10 particles and 500 iterations. BO is set with 50 random initialization, 500 iterations, and kappa $=$ 100 (higher values increase exploration).

The MAP estimation is obtained using a quasi-Newton L-BFGS optimizer from Python interface to Stan (Carpenter et al., 2017), an open-source probabilistic programming language for Bayesian modeling and statistical inference. We also use Stan’s implementation of the NUTS for fully Bayesian estimation by HMC, with a Metropolis acceptance rate of 0.8, maximum tree depth of 10, warm-up iterations of 1000, and 500 sampling phase iterations.

Finally, for SBI, we used SNPE from the PyTorch–based SBI toolkit (Tejero-Cantero et al., 2020). SNPE was run for a single round with 100,000 simulations, unless specified otherwise. For training step, we used an MAF, with five autoregressive layers, each with two hidden layers of 50 units. The set of data features includes the statistical moments of time series up to the fourth order (mean, standard deviation, skewness, and kurtosis) and the peak properties (number and location of first peak, as shown in Figure S2).

Parameter bounds for optimization and uniform prior for Bayesian inference were set as $Δ∈[0.1,5]$⁠, $η∈[-10,-3]$⁠, $J∈[5,20]$⁠. For the simulation of equation 2.1, we used Euler-Maruyama integration for $t=100sec$ with a time step of $dt=1msec$⁠. The ground truth parameters were set as $η=-4.6$⁠, $J=14.5$⁠, and $Δ=0.7$⁠, ensuring that the system is in a bistable regime, unless specified otherwise, with a one-step current with amplitude of $3v$ applied from time $30sec$ to time $60sec$⁠. The input current was considered known and not included in the inferred quantities unless otherwise specified. Each simulation took around $0.01sec$ to run using a just-in-time (JIT) compiler. To generate stochastic dynamics, a zero-mean gaussian noise with $σ=0.1$ was added to the both the $r$ and $v$ variables. The prior on dynamical noise was set as $σ∈[0,1]$⁠. For comparison with the MF model, we ran a network of $104$ all-to-all connected QIF neurons using the Brian simulator (Stimberg et al., 2019).

The model simulation and parameter estimation were performed on a Linux machine with 3.60 GHz Intel Core i7-7700 and 8 GB of memory.

We present the results of the optimization and Bayesian algorithms for three cases: (1) inferring from exact deterministic synthetic data, (2) inferring from stochastic synthetic data driven by dynamical (a.k.a. state) noise, and (3) inferring with stochastic synthetic data when the activity for $r$ is unknown (missing data). We report and compare the goodness-of-fit using the root mean square error (RMSE) to the true (observed or hidden) time series and parameters, the variance in the estimation process, and the computational cost.

Before running inference process, we conduct a sensitivity analysis to assess the structural identifiability of the parameters. This analysis involves calculating local sensitivity coefficients, which measure the effect of small adjustments in a parameter on the model’s output while keeping all other parameters constant. A small change in a sensitive parameter leads to a significant alteration in the model output, indicating its identifiability. Conversely, when there are no alterations in the model output despite adjustments in a parameter, it suggests the parameter’s nonidentifiability. To assess the identifiability of the parameters, we used the profile likelihood (Raue et al., 2009; Wieland et al., 2021) and Hessian matrix, a metric describing the local curvature of a function based on its second partial derivatives (Hashemi et al., 2018, 2023). Our results indicate that all three parameters $Δ$⁠, $η$⁠, and $J$ can be identifiable; however, the profile likelihood obtained using equation 2.2 reveals the presence of numerous local minima (see Figure S3). This highlights the difficulty in attaining the global minimum during the inference process.

Here, we compare the inference results of different algorithms when we have complete observations of both state variables $(r,v)$⁠, and the system operates without any noise (see Figure 1). The observed time series and trajectories in phase-plane, which exhibit a bistability between a stable node and a stable focus, driven by the input current, are illustrated in Figure 1A. From the results demonstrated in Figure 1B, it can be seen that all the algorithms qualitatively reproduce the bistability behavior in the phase-plane. However, when evaluating their performance by estimating the true generative parameters, all optimization methods except DE get stuck in a local minima (see Figure 1C). The results indicate that DE, SBI, and HMC algorithms correctly recover the ground-truth parameters, the former as an almost exact point estimate (see Figure 1C, top panels), while the latter two, SBI and HMC, yield full posterior distributions (see Figure 1C, bottom panels). In terms of uncertainty quantification, HMC offers a slightly more informative posterior distribution compared to SBI (see Figure S4).

When it comes to matching the observed time series, the optimization algorithms such as PSO and BO deviate significantly from true values, while other approaches equally retrieve an almost perfect fit to the observed time series (see Figure 1D, top panel). Nevertheless, MAP largely fails to accurately estimate model parameters due to overfitting (see Figure 1D, bottom panel). This type of overfitting emphasizes the importance of quantifying uncertainty to verify the reliability of inference, going beyond a point estimation. Note that for HMC and SBI, we report the results using posterior predictive check that is, re-generating data using random parameters drawn from the estimated posterior and then comparing simulations with the observed data.

In terms of computational cost for inference, DE has a clear advantage with its rapid performance, typically taking less than a minute to complete (see Figure 1E). Despite a high precision, the computational cost of HMC in this example is prohibitively expensive, taking almost 100 hours to complete the inference process. On the other hand, SBI was terminated in nearly one and a half hours (including 100,000 random simulations, training, and sampling), making it approximately 60 orders of magnitude faster than HMC. Additionally, due to the amortized approach adopted by SBI, each sampling process takes less than 1 minute to estimate the joint posterior distributions from new data.

In summary, DE demonstrates superior speed and accuracy, but it only provides a point estimate. Among Bayesian methods, HMC is exact and provides certain estimates, but it is computationally prohibitive. SBI effectively provides accurate estimates along with associated uncertainties at a reasonable computational cost (see Table 1).

The inherent randomness in stochastic systems introduces uncertainty into the behavior of the system under study, which makes it challenging to accurately identify its underlying dynamics. In particular, the presence of dynamical noise significantly increases the complexity of the inference process and necessitates the use of robust probabilistic methodologies that can effectively account for and handle the stochastic nature of the data.

Here, we report the results of inference on synthetic data with zero-centered gaussian dynamical noise, where the intensity is $σ=0.1$ (see Figure 2). The observed noisy time series and trajectories in phase-plane, exhibiting a bistable behavior between a stable node and a stable focus, are illustrated in Figure 2A. From the results demonstrated in Figure 2B, we can see that all the algorithms qualitatively replicate the bistable behavior in the phase-plane. Yet when assessing how well they recover the true parameters, only DE yields an accurate point estimate among the optimization algorithms (see Figure 2C, top panels).

In terms of uncertainty quantification, SBI generates posteriors that are tightly centered on the ground-truth parameters, in comparison to the HMC sampling (see Figure 2C, bottom panels and Figure S5). A detailed comparison of convergence diagnostics indicates that both SBI and HMC methods provide ideal Bayesian estimation (see Figure S6). Regarding the proximity to the observed time series and true parameters, both HMC and SBI offer a closer match compared to the other algorithms (see Figure 2D). This highlights the challenges inherent in inferring stochastic systems via optimization algorithms, as the error metrics such as RMSE used to define the objective function may not reliably gauge accuracy.

In terms of computational efficiency for inference, the running time of all algorithms increases when the noise is present, except for HMC (see Figure 2E). Nevertheless, when considering the entire process (including random simulations, training, and sampling), SBI remains approximately eight orders of magnitude faster than HMC. Interestingly, SBI is also able to accurately and efficiently estimate the dynamical noise in the system (see Figure S7).

In summary, these results demonstrate that when considering overall accuracy in the presence of noise and bistability, uncertainty quantification, and computational cost, the SBI outperforms all other algorithms, including HMC (see Table 2).

Next, we investigate how varying the intensity of dynamical noise affects the inference on stochastic data. Our results indicate that the overall quality of fit to the noisy data (see Figure 3A) and the accuracy of recovered parameters (see Figure 3B) remain stable up to $σ≤10-1$ for all algorithms. However, both metrics exhibit a significant increase beyond this threshold, particularly for optimization methods.

The MAP estimation is robust with regard to the amount of dynamical noise, but it tends to overfit. This is because the MAP estimation consistently provides the least accurate inference on parameters, even though its fit to the time-series data is almost perfect. In contrast, Bayesian inference methods such as HMC and SBI are overall more accurate and significantly more robust compared to optimization methods. Both HMC and SBI exhibit significantly lower error than others, even at high noise levels (e.g., $σ=1$⁠).

Interestingly, SBI appears to be more resilient to high levels of noise compared to HMC. Given a sufficient number of simulations for training (e.g., 100,000), SBI demonstrates the most accurate and robust fit to the data, consistently staying close to a perfect fit across various noise values. (See Figure S8 for a systematic investigation on the impact of the number of simulations on the performance of SBI.) Overall, these results indicate that SBI is the most accurate and robust algorithm in our benchmark.

Here, we explore the performance of inference methods in situations where data are available for only one of the variables in the state-space modeling. We consider the mean membrane potential $v$ as the observed data, which is simulated with dynamic noise of intensity $σ=0.1$⁠, while the firing rate $r$ remains latent (see Figure 4A). The noise intensity is fixed for inference process.

From Figure 4B, we observe that only DE and SBI were capable of retrieving the true dynamics in the phase-plane when $r$ was missing. As shown in Figure 4C (top panel), among optimization algorithms, only DE correctly estimates the true parameters. Notably, HMC fails considerably in this problem by proposing overly confident posterior distributions that are far from the ground truth (see Figure 4C, bottom panel). In contrast, the estimated posteriors using SBI are centered on the ground-truth parameters, but they exhibit a more diffuse uncertainty compared to the previous results (see Figure S9). Note that HMC fails to accurately reconstruct the latent variable $r$ from observed $v$ (see Figure S10), even though its error on the fitted trajectories is comparable to that of SBI (see Figure 4D, top panel). Nevertheless, the unreliability in the results produced by HMC becomes evident when observing the RMSE values for model parameters (see Figure 4D, bottom panel).

In terms of computational cost for inference (see Figure 1E), optimization by DE still has a clear advantage with its rapid performance. Interestingly, SBI remains efficient, approximately 68 orders of magnitude faster than HMC (see Table 3). Overall, SBI outperforms HMC in the recovery of the bistable dynamics, including the hidden firing rate. However, this approach can lead to an overestimation in the associated uncertainty, when compared to the full observed state-space dynamics (see Figure 2 versus Figure 4).

As demonstrated in the previous section, the inference process posed a challenge when only the membrane potential $v$ was observed and the firing rate $r$ was missing. To improve the inference in such cases, we approximately reconstruct the hidden $r$ from the observed counterpart $v$ using the time-delay embedding technique (see Figure 5).

As it can be seen from Figure 5A, the reconstructed firing rates closely follow the original time series (RMSE $=$ 0.151). This leads to accurately capturing the bistable switching behavior in the phase-plane, as shown in Figure 5B. We subsequently explored how the access to the trajectories of state variables influences the statistical relationships between parameters. When both the firing rate (⁠$r$⁠) and membrane potential (⁠$v$⁠) are observed, both HMC and SBI algorithms reveal a strong correlation between the average excitability $η$ and the synaptic weight $J$ (⁠$ρη,J≈-0.78$⁠), while the other parameters exhibit no codependency (see Figure 5C).

When the firing rate is latent and only the membrane potential is observed, the correlation between parameters is more pronounced, as estimated by HMC, while SBI tends to overestimate this codependency as fully degenerate (see Figure 5D). This highlights the challenges in the inference or parameter estimation process when lacking access to complete knowledge of the system dynamics. This issue can be improved by using time-delay embedding to reconstruct the full phase-space dynamics and inform the inference process, as shown by the reduced correlation between parameters in the joint Bayesian posterior in panel 5E.

In the previous sections, we performed inference against the system dynamics derived from the MF model given by equation 2.1. We now aim to investigate how the inference of the posterior distribution and the relationships between parameters remain consistent by traversing across scales. This can be achieved by conducting inference using observed data generated by QIF neurons (at the microscopic level) versus the data generated by the MF model (at the macroscopic level).

First, we compared the simulated membrane potential $v$ and firing rate $r$ using an MF model given by equation 2.1 to the averaged activities of a network of $104$ all-to-all connected QIF neurons (see the raster plot in Figure 6A). Figure 6B demonstrates that the MF model can accurately generate the (smoothed) transient dynamics that emerge from an ensemble of spiking neurons. However, the first spike emitted by the MF model after stimulation may exhibit a short lag compared to the averaged QIF neurons. Then we used the data generated by the MF and QIF models as the observation and conducted inference using the MF model through Bayesian estimation algorithms, such as HMC and SBI. Specifically, we compared the correlations between the estimated joint posteriors of model parameters in each case, as shown in Figure 6(C).

We observed consistent correlations between parameters across the two types of data sets (macroscopic MF and microscopic QIF) when fitted using the MF model. Our results indicate that the strong negative linear correlation between excitability and synaptic weight (⁠$ρη,J≈-0.78$⁠) persists when using HMC and SBI against both data sets. Considering the agreement between HMC and SBI across two data sets, we can conclude that the consistent high correlation between parameters $η$ and $J$ is intrinsic and not induced by the inference process or model assumptions.

Here, we challenge the SBI approach in inferring system dynamics while also accounting for an unknown input current, which plays a crucial role in emerging the bistable behavior. Given only the position or waveform of input currents, we estimate the intensity or angular velocity across various ground-truth values: $I0$ in a step current as $I(t)=I0·1stim(t)$ and $ω0$ in a sinusoidal current as $I(t)=sin(ω0t)·1stim(t)$⁠, as shown in Figures 7A, and 7B, respectively. Our results demonstrate that in both cases, the SBI approach accurately recovers the unknown parameters in the input currents. The posterior credibility intervals, visualized as error bars, indicate certain estimates that are close to a perfect fit (⁠$y=x$⁠, in red) across all different values. See Figure S11 for the observed and predicted time series. This validates the capability of SBI in accurately estimating the system dynamics, even when the characteristics of the input current are unknown.

Here, we show that SBI can be used to investigate the stability of system dynamics from low-dimensional summary statics of observed time-series. Figure 8 shows a phase diagram of the system as a function of the mean $η$ and synaptic weight $J$⁠, both normalized by the width of the input distribution $Δ$⁠. Using linear stability analysis, there are three qualitatively distinct regions of the phase diagram: (1) a single stable node corresponding to a low-activity state (shown in blue), (2) a single stable focus generally corresponding to a high-activity state (shown in red), and (3) a region of bistability between low and high firing rate (shown in cyan).

Interestingly, a similar basin of bistability in phase diagram can be readily reproduced using deep neural density estimators (such as MAF model) in the SBI approach. By training on the low-dimensional data features extracted from the time series (such as the presence or absence of damped oscillations before, during, and after stimulation), the generated posterior samples display a very close agreement with the results obtained from linear stability analysis. This demonstrates the capability of SBI in accurately estimating system dynamics from summary statics, including the presence of bistability in the phase diagram.

In this section, our aim is to infer the collective dynamics of QIF neurons without making any assumptions about the underlying generating dynamics. To achieve this, we used neural ODEs as a powerful tool for modeling continuous-time dynamics without assuming any prior knowledge of the underlying equations governing the system. Unlike traditional discrete-time models, neural ODEs parameterize the continuous-depth formulation, allowing for seamless interpolation between observed data points.

We first validate neural ODEs using the data garnered by MF model described by equation 2.1. We generated three different data sets using the MF model given by equation 2.1, with a set of parameters corresponding to a bistable regime: ${Δ=1,J=15,η=-5}$⁠. For each data set, we sampled the phase-space by varying initial conditions according to a regular grid so that $r0∈[0.1,3]$ and $v0∈[-2,2]$⁠. We then solved the system for 1000 time points (with an integration time step of $dt=0.01sec$⁠). To speed up training, each trajectory was downsampled by a factor of 10 and divided into 10 segments, each consisting of 10 time points. The training and test data were randomly split, with $75%$ of the data used for training and $25%$ used for testing.

The first training set, comprising $N=100$ deterministic trajectories, was generated (see Figure 9A). The results indicate that the neural ODE almost perfectly reconstructed the phase-space (see Figure 9B). In the second training set, which had the same size, a dynamical noise with a standard deviation of $σ=0.1$ was added during the integration process (see Figure 9C). In this case, the estimated nullclines suffered from overfitting although without affecting the overall reconstructed dynamics, as the predicted trajectories were still very similar to the original data (as shown in Figure 9D). Increasing the size of the training data set (see Figure 9E) significantly reduces overfitting. As a result, the reconstructed nullclines show very close agreement with those obtained from the deterministic data (see Figure 9F). To illustrate overfitting, Figure S12 shows the loss functions for the different scenarios, as well as snapshots of phase-space reconstruction during training. Overall, the neural ODE successfully reconstructs the underlying deterministic system, even in the presence of noise.

We then trained neural ODEs using data generated by $104$ QIF neurons with a uniform stimulus (see Figure 10). The data were partitioned using the first 400 points for training and predicting the remaining 1600 points. The results indicate that using the dynamics of derivatives, we can achieve a reliable understanding and prediction of the complex behavior of a network of spiking neurons, as illustrated in Figure 10. The emergent dynamics vary based on different parameter settings, resulting the stable node, stable focus, and bistable regime. (See Figure S13 for the loss function in the training and test sets.)

In the ever-evolving field of computational neuroscience, accurately estimating parameters that consistently govern the collective behavior of neural networks is a crucial endeavor, especially within the framework of recurrently coupled spiking neurons. In this study, we emphasize the use of mean-field (MF) theory to streamline the inference process for networks of spiking neurons. This choice is driven by the computational challenges in the calculation of the likelihood function—an essential ingredient for both frequentist and Bayesian inference methods—that becomes computationally prohibitive when attempting inference using these spiking networks in forward modeling.

The needs and objectives of researchers in the field of computational neuroscience are as diverse as the neural systems they aim to understand. Some studies may prioritize rapid point estimation through optimization, seeking quick outcomes to inform real-time decisions (Vattikonda et al., 2021; Penas et al., 2023). Other studies find value in exploring the full distribution of parameters, which provides a nuanced understanding of parameter uncertainty for reliable decision making (Hashemi et al., 2020; Jha et al., 2022). This article offers a comprehensive but not exhaustive benchmarking of the state-of-the-art inference methods applied to an MF model of spiking neurons (see Figures 1, 2, and 4). Our comparative analyses offer practical guidance, assisting researchers in selecting the most suitable method for their specific data sets and research inquiries. Note that the comparison of computational cost and accuracy in parameter estimation can heavily depend on the selection of hyperparameters, such as population size in DE, warm-up phase in HMC, or the number of simulations and features in SBI. To make an unbiased assessment, we used optimal values for these hyperparameters in each algorithm, tailored to the dynamical model used in this study. Nevertheless, these optimal values may vary for different inverse problems, depending on parameter space and data dimensions, nonlinearity, sparsity, and the mapping function to measurements.

While optimization and approximate Bayesian computing methods can construct a confidence interval for the estimation based on bootstraping or on a threshold for accepting or rejecting the estimates, the results are highly dependent on the chosen threshold value (Beaumont et al., 2002; Cranmer et al., 2020). In contrast, Bayesian inference naturally provides uncertainty quantification by placing a distribution over parameters and treating them as random variable. Following Bayes’s rule, this distribution is updated with evidence from observed data to form the posterior distribution, which furnishes comprehensive information for inference and prediction.

Our results indicated that evolutionary algorithms for solving global optimization problems, such as DE, provide rapid and accurate point estimation of the true generative parameters when there is no dynamical noise present. Challenges arise when using optimization methods in the presence of noisy data, mainly due to the lack of an efficient form of the objective function. The selection of the objective function plays a crucial role in determining the estimation through optimization methods (Svensson et al., 2012; Hashemi et al., 2018). The error explanation with distance metrics such as RMSE is limited in accurately capturing the underlying data generation process (Baldy et al., 2023). This is because when generative parameters remain unchanged but when dynamic noise is introduced, the time series can show large fluctuations, resulting in deviations from the observed data. In particular, the presence of noise can easily lead to unreliable estimations, increasing the risk of overfitting where the model fits to the noise rather than capturing the true underlying relationship. Consequently, it becomes more conspicuous to conduct inference using distributions in the Bayesian framework, particularly when dealing with dynamical noise. Diagnosing of overfitting, as shown using MAP estimation (see Figures 2 and 4), can be better understood through uncertainty quantification.

Furthermore, Bayesian inference can reveal the true relationships between parameters, capturing degeneracies in parameter space (Edelman & Gally, 2001; Hashemi et al., 2023). For example, when we assess the agreement between HMC and SBI across different scales (as shown in Figures 5, S4, S5, and S9), it can be concluded that the persistent strong correlation between parameters $η$ and $J$ is inherent and not influenced by the inference procedure or model assumptions. This is in line with previous findings that have reported a strong and robust correlation between firing rates and synaptic weights across different brain states, environments, and situations (Buzsáki & Mizuseki, 2014).

One of the main findings of this study is the effectiveness of deep neural networks in generating probability distributions for parameters of networks of spiking neurons. This approach outperforms other computational algorithms, such as MCMC, particularly in real-world applications involving missing data in bistable systems (see Figure 4). The effectiveness of deep generative models for inference from the mechanistic model of networks of spiking neurons is confirmed by the robustness of the estimation under significant dynamic noise (see Figures 3 and S7), as well as the precise estimation of input current (see Figure 7) and the consistency with linear stability analysis (see Figure 8). However, it is important to acknowledge that the use of deep generative models can result in an overestimation of uncertainty and correlations between parameters. To address this challenge, incorporating time-delay embedding is an effective remedy (see Figure 5).

Using high-performance computing, model simulations can be run independently, creating a large training data set for training deep neural density estimators in SBI approach (Hashemi et al., 2023). In contrast, HMC is limited to embarrassingly parallel execution with only independent chains on computational nodes (Hashemi et al., 2021). Moreover, when dealing with bistability in the state-space representation, HMC methods require significant computational time to detect state transitions in the latent space (see Tables 1, 2, and 3) or need to be augmented with generative models such as normalizing flows (Hoffman et al., 2019; Gabrié et al., 2022). On the other hand, SBI offers efficient Bayesian estimation, even without detailed knowledge of the system’s state-space representation. This aligns with findings from recent studies that highlight the efficiency of SBI across various challenging inverse problems (Gonçalves et al., 2020; Deistler et al., 2022; Boelts et al., 2022, 2023; Hashemi et al., 2023; Lavanga et al., 2023; Yalccinkaya et al., 2023; Rabuffo et al., 2023; Sorrentino et al., 2023). By benchmarking and addressing questions such as computational cost, uncertainty quantification, interdependency exploration, and data availability, we conclude that SBI is more efficient than alternatives in making informed choices from microscopic states to emergent dynamics at the macroscale.

Note that our results should not rule out the use of HMC for dynamical MF models of spiking neurons. The implementation of HMC with parameterization tricks, adaptive integrators to solve ODEs/SDEs, and an effective initialization strategy in other tools such as Numpyro (Phan et al., 2019) could significantly improve the computational cost of state-space estimation. Given that the main challenge lies in the complex posterior geometries (Betancourt, 2017) and integration process, substantial time savings can be achieved through reparameterization techniques and optimizing the integrator. For instance, a careful manifold reparameterization to reorganize the model configuration space or replacing carry-over for-loops with a more efficient evaluation method (Bradbury et al., 2018b) will facilitate exploration of the posterior distribution in terms of computational time and convergence diagnostics.

Moreover, our analyses were focused on a single neural population, with the assumption of heterogeneous all-to-all interaction between QIF neurons in the thermodynamic limit (Montbrió et al., 2015). Considering coupled excitatory and inhibitory populations and incorporating anatomical features of cortico-cortical connections would bring us closer to traversing scales, up to the whole-brain level (Hashemi et al., 2020, 2021). In essence, this presents a more challenging and sensitive inversion problem, and it is necessary to ascertain whether the competitive performance of SBI still holds in this context, especially considering the exploding number of parameters.

Considering the aforementioned points, extending the generality of our conclusions necessitates further investigation. The first part of this work has focused on machine learning procedures for state-space modeling, some with a Bayesian twist to approximate the posterior distribution in one way or another. The alternative approach for the latter, which predominates in the identification of MF approximations to neuronal networks, is based on variational procedures, as used in dynamic causal modeling (DCM; Friston et al., 2003). DCM rests on minimizing the free energy, that is, the Kullback-Leibler divergence between the true and approximate posterior, minus the log evidence (or marginal likelihood). This approach enables analytic solutions to the approximated posterior, eliminating the need for sampling, and provides a variational bound on model evidence for Bayesian model comparison (Penny, 2012; Zeidman et al., 2023). Hence, when the computed posterior is sufficiently accurate, free energy serves as a reliable approximation to the negative log evidence. From the perspective of system identification (Linderman et al., 2017), DCM can be regarded as solving a triple estimation problem, namely, inferring the latent states, parameters, and precisions of a state-space model, with a known functional form (Friston et al., 2010; Schiff & Sauer, 2007). Consequently, this approach to network inference can be interpreted as predictive coding (Millidge et al., 2020), which itself serves as a model or description of MF approximations to neural networks in the brain (Friston & Kiebel, 2009).

The main body of this work relies on the modeling assumption of mean-field approximation to spiking neurons. This assumption prompts questions regarding the suitability of the chosen MF approximation. Such inquiries are relevant to system identification, which is often framed in terms of model selection, structure learning, or network discovery (Friedman & Koller, 2004; Gershman & Niv, 2010; Seghier & Friston, 2013; Wipf & Rao, 2007). For instance, one might explore the possibility of replacing equation 2.1 with a chaotic system (augmented with suitable exogenous inputs) and evaluate whether the MF model given by equation 2.1 offers a superior explanation for the data compared to the chaotic system. Assessing model evidence can be achieved through Bayesian model comparison by computing information criteria from a Bayesian perspective, such as the widely applicable information criterion (Watanabe, 2010) and expected log predictive density using leave-one-out cross-validation (Vehtari et al., 2016; Gelman et al., 2013). In this context, a thorough investigation of approaches such as SBI is necessary in future studies. Nevertheless, HMC may be indispensable for achieving greater accuracy in Bayesian model comparison.

Our study highlights the use of deep neural ODEs in inferring vector fields at a macroscopic level, enabling the prediction of system dynamics from microscopic states. This approach has the potential to make interpretable predictions at larger scales from simulations at a detailed level, aiding in the prognosis and diagnosis of brain diseases. Recently, Sip et al. (2023) introduced a method using variational autoencoders for nonlinear dynamical system identification at the whole-brain level to infer both the neural mass model and the region- and subject-specific parameters from the functional data while respecting the known network structure. The scalability of neural ODEs at the whole-brain network level remains to be investigated in future studies. In addition, symbolic regression applied to the outcomes from neural ODEs may unveil closed-form equations for neural mass models, offering promising avenues for future research. This approach may lead to the discovery of concise data-driven MF representations of complex neural dynamics, contributing to our understanding of brain function.

In conclusion, this work highlights the improved accuracy and efficiency that deep learning techniques bring to the inference from networks of spiking neurons. It opens up exciting possibilities for future research in neural computation, where the trade-off between accuracy and uncertainty needs to be carefully considered. As we continue to explore the capabilities of SBI and neural ODEs at larger scales, this study serves as a valuable step forward in our quest to unravel the complexities of neural networks and their computational mechanisms.

All code is available on GitHub (https://github.com/ins-amu/Inference_MFM).

This research has received funding from EU’s Horizon 2020 Framework Programme for Research and Innovation under the Specific Grant Agreements 101147319 (EBRAINS 2.0 Project), 101137289 (Virtual Brain Twin Project), and government grant managed by the Agence Nationale de la Recherche reference ANR-22-PESN- 0012 (France 2030 program). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. We thank Lionel Kusch, Matthieu Gilson, and Spase Petkoski for fruitful discussions.

Conceptualization: M.W., V.J., and M.H. Methodology: M.W. and M.H. Software: N.B, M.B., M.W., and M.H. Investigation: N.B., M.B., and M.H. Visualization: N.B., M.B., and M.H. Supervision: V.J. and M.H. Funding acquisition: V.J. Writing of the original draft: N.B., and M.H. Writing review and editing: N.B, M.B, M.W., V.J, and M.H.