1. Introduction

Understanding the complex dynamics of the brain and their neurobiological underpinnings, with the potential to advance precision medicine [1, 2, 3, 4], is a central goal in neuroscience. Modeling these dynamics provides crucial insights into causality and mechanisms underlying both normal brain function and various neurological disorders [5, 6, 7]. By integrating the average activity of large populations of neurons (e.g., neural mass models; [8, 9, 10, 11, 12, 13]) with information provided by structural imaging modalities (i.e., connectome; [14, 15, 16, 17]), whole-brain network modeling has proven to be a powerful tractable approach for simulating brain activities and emergent dynamics as recorded by functional imaging modalities (such as (s)EEG/MEG/fMRI; [18, 19, 20, 21, 22, 23].

The whole-brain models have been well-established in network neuroscience [24, 25] for understanding the brain structure and function [26, 27, 28, 29, 30, 31, 32] and investigating the mechanisms underlying brain dynamics at rest [33, 34, 35, 36], normal aging [37, 38], and also altered states such as anaesthesia and loss of consciousness [39, 40, 41, 42]. This class of computational models, also known as virtual brain models [43, 44, 18, 19, 45, 46], have shown remarkable capability in delineating the pathophysiological causes of a wide range of brain diseases, such as epilepsy [47, 48, 6], multiple sclerosis [49, 50], Alzheimer’s disease [51, 52], Parkinson’s disease [53, 54], neuropsychiatric disorders [55, 56], stroke [57, 58] and focal lesions [59]. In particular, they enable the personalized simulation of both normal and abnormal brain activities, along with their associated imaging recordings, thereby stratifying between healthy and diseased states [60, 61, 52] and potentially informing targeted interventions and treatment strategies [47, 62, 6, 45, 23]. Although there are only a few tools available for forward simulations at the whole-brain level, e.g., the brain network simulator The Virtual Brain (TVB; [44]), there is a lack of tools for addressing the inverse problem, i.e., finding the set of control (generative) parameters that best explains the observed data. This study aims to bridge this gap by addressing the inverse problem in large-scale brain networks, a crucial step toward making these models operable for clinical applications.

Accurately and reliably estimating the parameters of whole-brain models remains a formidable challenge, mainly due to the high-dimensionality and non-linearity inherent in brain activity data, as well as the non-trivial effects of noise and network inputs. A large number of previous studies in whole-brain modeling have relied on optimization techniques to identify a single optimal value from an objective function, scoring the model’s performance against observed data [34, 36, 63, 64]. This approach often involves minimizing metrics such as the Kolmogorov-Smirnov distance or maximizing the Pearson correlation between observed and generated data features such as functional connectivity (FC), functional connectivity dynamics (FCD), and/or power spectral density (PSD). Although fast, such a parametric approach results in only point estimates and fails to capture the relationship between parameters and their associated uncertainty. This limits the generalizability of findings, and hinders identifiability analysis, which explores the uniqueness of solutions. Furthermore, optimization algorithms can easily get stuck in local extrema, requiring multi-start strategies to address potential parameter degeneracies. These additional steps, while necessary, ultimately increase the computational cost. Critically, the estimation heavily depends on the form of the objective function defined for optimization [65, 66]. These limitations can be overcome by employing Bayesian inference, which naturally quantifies the uncertainty in the estimation and statistical dependencies between parameters, leading to more robust and generalizable models. Bayesian inference is a principal method for updating prior beliefs with information provided by data through the likelihood function, resulting in a posterior probability distribution that encodes all the information necessary for inferences and predictions. This approach has proven essential for understanding the intricate relationships between brain structure and function [67, 59, 37], as well as for revealing the pathophysiological causes underlying brain disorders [68, 51, 49, 46, 69].

In this context, simulation-based inference (SBI; [70, 71, 68, 69]) has gained prominence as an efficient methodology for conducting Bayesian inference in complex models where traditional inference techniques become inapplicable. SBI leverages computational simulations to generate synthetic data and employs advanced probabilistic machine learning methods to infer the joint distribution over parameters that best explain the observed data, along with associated uncertainty. This approach is particularly well-suited for Bayesian inference on whole-brain models, which often exhibit complex dynamics that are difficult to retrieve from neuroimaging data with conventional estimation techniques. Crucially, SBI circumvents the need for explicit likelihood evaluation and the Markovian (sequential) property in sampling. Markov chain Monte Carlo (MCMC) is a standard non-parametric technique and asymptotically exact for sampling from a probability distribution [72]. However, for Bayesian inference on whole-brain models given high-dimensional data, the likelihood function becomes intractable, rendering MCMC sampling computationally prohibitive. SBI offers significant advantages, such as parallel simulation while leveraging amortized learning, making it effective for personalized inference from large datasets [69]. Amortization strategies based on artificial neural networks allow for immediate application to arbitrary inference on new recordings without the need for repeated training [73]. Following an initial computational cost during simulation and training to learn all posterior distributions, subsequent evaluation of new hypotheses can be conducted efficiently, without additional computational overhead for further simulations [68]. Importantly, SBI sidesteps the convergence issues caused by complex geometries that are often encountered when using gradient-based MCMC methods [74, 75, 76]. It also substantially outperforms ABC methods, which rely on a threshold to accept or reject samples [77, 78, 71]. Such a likelihood-free approach provides us with generic inference on complex systems as long as we can provide three modules:

1. A prior distribution, describing the possible range of parameters from which random samples can be easily drawn, i.e., .

2. A simulator in computer code, that takes parameters as input and generates data as output, i.e., .

3. A set of low-dimensional data features, which are informative of the parameters that we aim to infer.

These elements prepare us with a training data set with a budget of N_sim simulations. Then, using a class of deep neural density estimators (such as masked autoregressive flows [79] or neural spline flows [80]), we can approximate the posterior distribution of parameters given a set of observed data, i.e., . Therefore, a versatile toolkit should be flexible and integrative, adeptly incorporating these modules to enable efficient Bayesian inference over complex models.

To address the need for widely applicable, reliable, and efficient parameter estimation from different neuroimaging modalities, we introduce Virtual Brain Inference (VBI), a flexible and integrative toolkit for probabilistic inference at whole-brain level. This open-source toolkit offers fast simulation through just-in-time (JIT) compilation of various brain models in different programming languages (Python/C++) and devices (CPUs/GPUs). It supports space-efficient storage of simulated data (HDF5/NPZ/PT), provides a memory-efficient loader for batched data, and facilitates the extraction of low-dimensional data features (FC/FCD/PSD). Additionally, it enables the training of deep neural density estimators (MAFs/NSFs), making it a versatile tool for inference on various neuroimaging data types ((s)EEG/MEG/fMRI). VBI leverages high-performance computing, significantly enhancing computational efficiency through parallel processing of large-scale datasets, which would be impractical with current alternative methods. Although SBI has been used for low-dimensional parameter spaces [71, 49, 81], we demonstrate that it can scale to whole-brain models with high-dimensional known parameters across neuroimaging modalities, as long as informative data features are provided. VBI is now accessible on the cloud platform EBRAINS (https://ebrains.eu), enabling users to explore more realistic brain dynamics underlying brain (dys)functioning using Bayesian inference. In the following sections, we will provide an overview of the theoretical foundations of whole-brain models and simulation-based inference, describe the architecture and features of the VBI toolkit, and demonstrate the validation through a series of case studies using in-silico data. We explore various whole-brain models corresponding to different types of brain recordings: a whole-brain network model of Wilson-Cowan [8], Jansen-Rit [82, 83], and Stuart-Landau [84] for simulating neural activity associated with EEG/MEG signals, the Epileptor [11] related to stereoelectro-EEG (sEEG) recordings, and Montbrió [12], and Wong-Wang [85, 86] mapped to fMRI BOLD signals. Although these models represent source signals and could be applied to other modalities (e.g., Stuart-Landau representing generic oscillatory dynamics), we focused on their capabilities to perform optimally in specific contexts. For instance, some are better suited for encephalographic signals (e.g., EEG/MEG) due to their ability to preserve spectral properties, while others have been used for fMRI data, emphasizing their ability to capture dynamic features such as bistability and time-varying functional connectivity.

2. Materials and methods

2.1. The virtual brain models

To build a virtual brain model (see Figure 1), the process begins with parcellating the brain into regions using anatomical data, typically derived from T1-MRI scans. Each region, represented as nodes in the network, is then equipped with a neural mass model to simulate the collective behavior of neurons within that area. These nodes are interconnected using a structural connectivity (SC) matrix, typically obtained from diffusion-weighted magnetic resonance imaging (DW-MRI). The entire network of interconnected nodes is then simulated using neuroinformatic tools, such as The Virtual Brain (TVB; [18]), replicating the intricate dynamics of brain activity and its associated brain imaging signals ((s)EEG/MEG/fMRI). This approach offers insights into both normal brain function and neurological disorders. In the following, we describe commonly used whole-brain network models corresponding to different types of neuroimaging recordings.

Fig. 1.

2.1.1. Wilson-Cowan model

The Wilson-Cowan model [8] is a seminal neural mass model that describes the dynamics of connected excitatory and inhibitory neural populations, at cortical microcolumn level. It has been widely used to understand the collective behavior of neurons and simulate neural activities recorded by methods such as local field potentials (LFPs) and EEG. The model effectively captures phenomena such as oscillations, wave propagation, pattern formation in neural tissue, and responses to external stimuli, offering insights into various brain (dys)functions, particularly in Parkinson’s disease [87, 88].

We focused on a simplified model for generation of beta oscillation within the cortex-subthalamic nucleus-globus pallidus network [89]. The model incorporates a closed-loop connection from the STN back to the cortex, represented by a single inhibitory connection with a time delay. However, it does not include feedback via the indirect pathway (cortex-striatum-GPe), as experimental evidence suggests this pathway is not essential for generating beta oscillations [90]. Instead, the GPe receives a constant inhibitory input from the striatum, consistent with observations from Parkinson’s disease models:

where the functions S, G, E, and I represent the firing rates of the STN, GPe, and the excitatory and inhibitory populations, respectively. The parameters T_ij denote the synaptic connection time delays from population i to population j, while T_ii represents the time delay of self-connections. The synaptic weights, w_ij, follow the same subscript conventions as the time delays, indicating the influence of the presynaptic neuron’s firing rate on the postsynaptic neuron. The membrane time constants are denoted by τ_i. A constant input, C, is provided to the excitatory population in the cortex to account for a constant component of both extrinsic and intrinsic excitatory inputs, while Str represents the constant inhibitory input from the striatum to the GPe. Lastly, F_i are the activation functions. The nominal parameter values and the prior range for the target parameters are summarized in Table 1.

Table 1.

2.1.2. Jansen-Rit model

The Jansen-Rit neural mass model [91] has been widely used to simulate physiological signals from various recording methods like intracranial LFPs, and scalp MEG/EEG recordings. For example, it has been shown to recreate responses similar to evoked-related potentials after a series of impulse stimulations [83, 92], generating high-alpha and low-beta oscillations (with added recurrent inhibitory connections and spike-rate modulation) [93], and also seizure patterns similar to those seen in temporal lobe epilepsy [94]. This biologically motivated model comprises of three main populations of neurons: excitatory pyramidal neurons, inhibitory interneurons and excitatory interneurons. These populations interact with each other through synaptic connections, forming a feedback loop that produces oscillatory activity governed by a set of nonlinear ordinary differential equations [82, 83, 95]:

where P(t) represents an external input current. The model’s output at i_th region, which corresponds to the membrane potential of pyramidal cells, is given by y_1i −y_2i. The nominal parameter values and the prior range for the target parameters are summarized in Table 2.

Table 2.

2.1.3. Stuart-Landau oscillator

The Stuart-Landau oscillator [84] is a generic mathematical model used to describe oscillatory phenomena, particularly those near a Hopf bifurcation, which is often employed to study the nonlinear dynamics of neural activity [96, 97, 63, 49]. One approach uses this model to capture slow hemodynamic changes in BOLD signal [96], while others apply it to model fast neuronal dynamics, which can be linked directly to EEG/MEG data [97, 63, 49]. Note that this is a phenomenological frame-work, and both applications operate on completely different time scales.

In the network, each brain region, characterized by an autonomous Stuart-Landau oscillator, can exhibit either damped or limit-cycle oscillations depending on the bifurcation parameter a. If a < 0, the system shows damped oscillations, similar to a pendulum under friction. In this regime, the system, when subjected to perturbation, relaxes back to its stable fixed point through damped oscillations with an angular frequency ω₀. The rate of amplitude damping is determined by |a|. Conversely, if a > 0, the system supports limit cycle solutions, allowing for self-sustained oscillations even in the absence of external noise. At a = 0, the system undergoes a Hopf bifurcation, i.e., small changes in parameters can lead to large variations in the system’s behavior.

Using whole-brain network modeling of EEG/MEG data [49, 63], the oscillators are interconnected via white-matter pathways, with coupling strengths specified by subject-specific DTI fiber counts, i.e., elements of SC matrix. This adjacency matrix is then scaled by a global coupling parameter G. Note that coupling between regions accounts for finite conduction times, which are often estimated by dividing the Euclidean distances between nodes by an average conduction velocity T_jk = d_jk/v. Knowing the personalized time-delays [98, 99], we can use the distance as a proxy, assuming a constant propagation velocity. The distance itself can be defined as either the length of the tracts or the Euclidean distance. Taking this into account, the activity of each regions is given by a set of complex differential equations:

where Z is a complex variable, and Re[Z(t)] is the corresponding time series. In this particular realization, each region has a natural frequency of 40 Hz (ω_j= ω₀ = 2π·40 rad/s), motivated by empirical studies demonstrating the emergence of gamma oscillations from the balance of excitation and inhibition, playing a role in local circuit computations [100].

In this study, for the sake of simplicity, a common cortico-cortical conduction velocity is estimated, i.e., the distance-dependent average velocities v_jk = v. We also consider a = −5, capturing the highly variable amplitude envelope of gamma oscillations as reported in experimental recordings [101, 63]. This choice also best reflects the slowest decay time constants of GABA_B inhibitory receptors–approximately 1 second [102]. A Gaussian noise with an intensity of σ = 10⁻⁴ is added to each oscillator to mimic stochastic fluctuations. The nominal parameter values and the prior range for the target parameters are summarized in Table 3.

Table 3.

2.1.4. Epileptor model

In personalized whole-brain network modeling of epilepsy spread [47], the dynamics of each brain region are governed by the Epileptor model [11]. The Epileptor model provides a comprehensive description of epileptic seizures, encompassing the complete taxonomy of system bifurcations to simultaneously reproduce the dynamics of seizure onset, progression, and termination [103]. The full Epileptor model comprises five state variables that couple two oscillatory dynamical systems operating on three different time scales [11]. Then motivated by Synergetic theory [104, 105] and under time-scale separation [106], the fast variables rapidly collapse on the slow manifold [107], whose dynamics is governed by the slow variable. This adiabatic approximation yields the 2D reduction of whole-brain model of epilepsy spread, also known as the Virtual Epileptic Patient (VEP) as follows:

where x_i and z_i indicate the fast and slow variables corresponding to i_th brain region, respectively, and the set of unknown η_i is the spatial map of epileptogenicity to be estimated. In real-world epilepsy applications [67, 68, 6], we compute the envelope function from sEEG data to perform inference. The nominal parameter values and the prior range for the target parameters are summarized in Table 4.

Table 4.

2.1.5. Montbrió model

The exact macroscopic dynamics of a specific brain region (represented as a node in the network) can be analytically derived in thermodynamic limit of infinitely all-to-all coupled spiking neurons [12] or Θ neuron representation [108]. By assuming a Lorentzian distribution on excitabilities in large ensembles of quadratic integrate- and- fire neurons with synaptic weights J and a half-width Δ centered at η, the macroscopic dynamics has been derived in terms of the collective firing activity and mean membrane potential [12]. Then, by coupling the brain regions via an additive current (e.g., in the average membrane potential equations), the dynamics of the whole-brain network can be described as follow [109, 110]:

where v_i and r_i are the average membrane potential and firing rate, respectively, at the i_th brain region, and parameter G is the network scaling parameter that modulates the overall impact of brain connectivity on the state dynamics. The SC_ij denotes the connection weight between i_th and j_th regions, and the dynamical noise ξ(t) ∼ N(0, σ²) follows a Gaussian distribution with mean zero and variance σ².

The model parameters are tuned so that each decoupled node is in a bistable regime, exhibiting a down-state stable fixed point (low-firing rate) and an up-state stable focus (high-firing rate) in the phase-space [12, 81]. The bistability is a fundamental property of regional brain dynamics to ensure a switching behavior in the data (e.g., to generate FCD), that has been recognized as representative of realistic dynamics observed empirically [109, 111, 112].

The solution of the coupled system yields a neuroelectric dataset that describes the evolution of the variables (r_i(t), v_i(t)) in each brain region i, providing measures of macroscopic activity. The surrogate BOLD activity for each region is then derived by filtering this activity through the Balloon-Windkessel model [113]. The input current I_stim represents the stimulation to selected brain regions, which increase the basin of attraction of the up-state in comparison to the down-state, while the fixed points move farther apart [109, 111, 112].

The nominal parameter values and the prior range for the target parameters are summarized in Table 5.

Table 5.

2.1.6. Wong-Wang model

Another commonly used whole-brain model for simulation of neural activity is the so-called parameterized dynamics mean-field (pDMF) model [114, 36, 86]. At each region, it comprises a simplified system of two non-linear coupled differential equations, motivated by the attractor network model, which integrates sensory information over time to make perceptual decisions, known as Wong-Wang model [85]. This biophysically realistic cortical network model of decision making then has been simplified further into a single-population model [86], which has been widely used to understand the mechanisms underpinning brain resting state dynamics [36, 115, 38]. The pDMF model has been also used to study whole-brain dynamics in various brain disorders, including Alzheimer’s disease [116], schizophrenia [117], and stroke [118]. The pDMF model equations are given as:

where H(x_i) and S_i, and x_i denote the population firing rate, the average synaptic gating variable, and the total input current at the i_th brain region, respectively. ξ_i(t) is uncorrelated standard Gaussian noise and the noise amplitude is controlled by σ. The nominal parameter values and the prior range for the target parameters are summarized in Table 6.

Table 6.

According to recent studies [36, 38], we can parameterize the set of w, I and σ as linear combinations of group-level T1w/T2w myelin maps [119] and the first principal gradient of functional connectivity:

where Mye_i and Grad_i are the average values of the T1w/T2w myelin map and the first FC principal gradient, respectively, within the i_th brain region. Therefore, the set of unknown parameters to estimate includes G and linear cofficients (a_w, b_w, c_w, a_I, b_I, c_I, a_σ, b_σ, c_σ) ∈ ℝ⁹, hence 10 parameters in total.

2.1.7. The Balloon-Windkessel model

The Balloon-Windkessel model is a biophysical framework that links neural activity to the BOLD signals detected in fMRI. This is not a neuronal model but rather a representation of neurovascular coupling, describing how neural activity influences hemodynamic responses. The model is characterized by two state variables: venous blood volume (v) and deoxyhemoglobin content (q). The system’s input is blood flow (f_in), and the output is the BOLD signal (y):

where V₀ represents the resting blood volume fraction, E₀ is the oxygen extraction fraction at rest, ɛ is the ratio of intra- to extravascular signals, r₀ is the slope of the relationship between the intravascular relaxation rate and oxygen saturation, ϑ₀ is the frequency offset at the surface of a fully deoxygenated vessel at 1.5T, and TE is the echo time. The dynamics of venous blood volume v and deoxyhemoglobin content q are governed by the Balloon model’s hemodynamic state equations:

where τ₀ is the transit time of blood flow, α reflects the resistance of the venous vessel (stiffness), and f(t) denotes blood inflow at time t, given by

where s is an exponentially decaying vasodilatory signal defined by

where, ɛ represents the efficacy of neuronal activity x(t) (i.e., integrated synaptic activity) in generating a signal increase, τ_s is the time constant for signal decay, and τ_f is the time constant for autoregulatory blood flow feedback [113]. For parameter values, see Table 7, taken from [113, 120, 121]. The resulting time series is downsampled to match the TR value in seconds.

Table 7.

2.2. Simulation-based inference

In the Bayesian framework [122], parameter estimation involves quantifying and propagating uncertainty through probability distributions placed on the parameters (prior information before seeing data), which are updated with information provided by the data (likelihood function). The formidable challenge to conducting efficient Bayesian inference is evaluating the likelihood function p(x | θ). This typically involves intractable integrating over all possible trajectories in the latent space: p(x | θ) = ∫ p(x, z | θ)dz, where p(x, z | θ) is the joint probability density of the data x and latent variables z, given parameters θ. For whole-brain network models with high-dimensional and nonlinear latent spaces, the computational cost can be prohibitive, making likelihood-based inference with MCMC sampling challenging to converge [76, 67, 123].

Simulation-based inference (SBI; [70, 71, 68]), or likelihood-free inference [124, 125, 126], addresses issues with inference on complex systems by using deep networks to learn an invertible transformation between parameters and data, avoiding direct sampling and likelihood evaluation. Given a prior distribution ) placed on the parameters , N random simulations are generated (with samples from prior), resulting in pairs , where and is the simulated data given . By training a deep neural density estimator F (such as normalizing flows [79, 80]), we can approximate the posterior with by minimizing the loss function

over network parameters ϕ. Once the parameters of the neural network F are optimized, for observed data we can readily estimate the target posterior . This allows for rapid approximation and sampling from the posterior distribution for any new observed data through a forward pass in the trained network [68, 69].

This approach uses a class of generative machine learning models known as normalizing flows ([127, 128]) to transform a simple based distribution into any complex target through a sequence of invertible mappings. Here, generative modeling is an unsupervised machine learning method for modeling a probability distribution given samples drawn from that distribution. The state-of-the-art normalizing flows, such as masked autoregressive flows (MAF; [79]) and neural spline flows (NSF; [80]), enable fast and exact density estimation and sampling from high-dimensional distributions. These models learn mappings between input data and probability densities, capturing complex dependencies and multi-modal distributions [128, 129].

In our work, we integrate the implementation of these models from the open-source SBI tool, leveraging both MAF and NSF architectures. The MAF model comprises 5 flow transforms, each with two blocks and 50 hidden units, tanh nonlinearity and batch normalization after each layer. The NSF model consists of 5 flow transforms, two residual blocks of 50 hidden units each, ReLU nonlinearity, and 10 spline bins. We apply these generative models to virtual brain simulations conducted with random parameters, to approximate the full posterior distribution of parameters from low-dimensional data features. Note that we employ a single round of SBI to benefit from amortization strategy rather than using a sequential approach that is designed to achieve a better fit but only for a specific dataset [68, 69].

2.3. Evaluation of posterior fit

To assess the reliability of Bayesian inference using synthetic data, we evaluate the posterior z-scores (denoted by z) against the posterior shrinkage (denoted by s), as defined by [75]:

where and θ^∗ are the posterior mean and the true values, respectively, σ_prior is the standard deviation of the prior, and σ_post is the standard deviation of the posterior. The posterior z-score measures how well the posterior distribution of a parameter captures its true value, while the posterior shrinkage quantifies the contraction of the posterior distribution relative to the prior. A concentration of estimations towards large shrinkages indicates that the posteriors are well-identified, whereas a concentration towards small z-scores signifies that the true values are accurately captured within the posteriors.

2.4. Flexible simulator and model building

A key feature of the VBI pipeline is its modularity and flexibility in integrating various simulators (see Figure 2). The Simulation module of the VBI pipeline is designed to be easily interchangeable, allowing researchers to replace it with other simulators, such as TVB [18] (see Figure S10), Neurolib [130], Brian [131], Brainpy [132]. This adaptability supports a wide range of simulation needs and computational environments, making the VBI a versatile tool for inference in system neuroscience. In particular, Simulation module offers a comprehensive implementation of commonly used whole-brain models. This is a customized version of implementation from open-source TVB simulator. While VBI does not encompass all the features of the original TVB, it is mainly designed to leverage the computational power of GPU devices and significantly reduce RAM requirements (see Figure S1A). This optimization ensures that high-performance clusters can be fully utilized, enabling parallel and scalable simulations, as often is required to perform scalable SBI.

Figure 2.

2.5. Comprehensive feature extraction

VBI offers a comprehensive toolbox for feature extraction across various datasets. The Features module includes but not limited to: (1) Statistical features, including mean (average of elements), variance (spread of the elements around mean), kurtosis (tailedness of the distribution of elements), and skewness (the asymmetry of the distribution of elements), that can be applied to any matrix. (2) Spectral features, such as low-dimensional summary statistics of power spectrum density (PSD). (3) Temporal features, including zero crossings, area under the curve, average power, and envelope. (4) Connectivity features, including functional connectivity (FC), which represents the statistical dependencies or correlations between activity patterns of different brain regions, and functional connectivity dynamics (FCD), which captures the temporal variations and transitions in these connectivity patterns over time. These calculations are performed for the whole-brain and/or subnetwork (e.g., limbic system, resting state networks). However, since these matrices are still high-dimensional, we use standard dimensional reduction techniques, such as principal component analysis (PCA) on FC/FCD matrices, to extract their associated low-dimensional summary statistics. (5) Information theory features, such as mutual information and transfer entropy.

In this study, we use the term spatio-temporal data features to refer to the statistical features and temporal features derived from time series, while we refer to the connectivity features extracted from FC/FCD matrices as functional data features.

The Features module uses parallel multiprocessing to speed up feature calculation. Additionally, it provides flexibility for users to add their own custom feature calculations with minimal effort and expertise, or to adjust the parameters of existing features based on the type of input time series. The feature extraction module is designed to be interchangeable with existing feature extraction libraries such as tsfel [133], pyspi [134], hctsa, [135], and scikit-learn [136]. Note that, some lightweight libraries such as catch22 [137] are directly accessible from the VBI feature extraction module.

2.6. VBI workflow

Figure 1 illustrates an overview of our approach in VBI, which combines virtual brain models and simulation-based inference to make probabilistic predictions on brain dynamics from neuroimaging recordings. The inputs to the pipeline include the structural imaging data (for building the connectome), functional imaging data such as (s)EEG/MEG, and fMRI as the target for fitting, and prior information as a plausible range over control parameters for generating random simulations. The main computational costs involve model simulations and data feature extraction. The output of the pipeline is the joint posterior distribution of control parameters (such as excitability, synaptic weights, or effective external input) that best explains the observed data. Since the approach is amortized (i.e., it learns across all combinations in the parameter space), it can be readily applied to any new data from a specific subject.

In the first step, non-invasive brain imaging data, such as T1-weighted MRI and Diffusion-weighted MRI (DW-MRI), are collected for a specific subject (Figure 1A). T1-weighted MRI images are processed to obtain brain parcellation, while DW-MRI images are used for tractography. Using the estimated fiber tracts and the defined brain regions from the parcellation, the connectome (i.e., the complete set of links between brain regions) is constructed by counting the fibers connecting all regions. The SC matrix, with entries representing the connection strength between brain regions, forms the structural component of the virtual brain which constrains the generation of brain dynamics and functional data at arbitrary brain locations (e.g., cortical and subcortical structures).

Subsequently, each brain network node is equipped with a computational model of average neuronal activity, known as neural mass models (see Figure 1B and subsection 2.1). They can be represented in the generic form of a dynamical model as , with the system variables (such as membrane potential and firing rate), the control parameters (such as excitability), and the input current I_input(such as stimulation). This integration of mathematical mean-field modeling (neural mass models) with anatomical information (connectome) allows us to efficiently analyze functional neuroimaging modalities at the whole-brain level.

To quantify the posterior distribution of control parameters given a set of observations, , we first need to define a plausible range for the control parameters based on background knowledge , i.e., a simple base distribution known as a prior. We draw random samples from the prior and provide them as input to the VBI simulator (implemented by Simulation module) to generate simulated time series associated with neuroimaging recordings, as shown in Figure 1C. Subsequently, we extract low-dimensional data features (implemented by Features module), as shown in Figure 1D for FC/FCD/PSD, to prepare the training dataset , with a budget of N_sim simulations. Then, we use a class of deep neural density estimators, such as MAF or NSF models, as schematically shown in Figure 1E, to learn all the posterior . Finally, we can readily sample from , which determines the probability distribution in parameter space that best explains the observed data. Figure 2 depicts the structure of the VBI toolkit, which consists of three main modules. The first module, referred to as Simulation module, is designed for fast simulation of whole-brain models, such as Wilson-Cowan (subsubsection 2.1.1), Jansen-Rit (subsubsection 2.1.2), Stuart-Landau (subsubsection 2.1.3), Epileptor (subsubsection 2.1.4), Montbrió (subsubsection 2.1.5), and Wong-Wang (subsubsection 2.1.6). These whole-brain models are implemented across various numerical computing libraries such as Cupy (GPU-accelerated computing with Python), C++ (a high-performance systems programming language), Numba (a just-in-time compiler for accelerating Python code), and PyTorch (an open-source machine learning library for creating deep neural network).

The second module, Features, provides a versatile tool for extracting low-dimensional features from simulated time series (see subsection 2.5). The features include, but are not limited to, spectral, temporal, connectivity, statistical, and information theory related features, and the associated summary statistics. The third module focuses on Inference, i.e., training the deep neural density estimators, such as MAF and NSF (see subsection 2.2), to learn the joint posterior distribution of control parameters.

3. Results

In the following, we demonstrate the capability of VBI for inference on the state-of-the-art whole-brain network models using in-silico testing, where the ground truth is known. We apply this approach to simulate neural activity and associated measurements, including (s)EEG/MEG, and fMRI, while also providing diagnostics for the accuracy and reliability of the estimation. Note that for (s)EEG/MEG neuroimaging, we perform inference at the regional level rather than at the sensor level, whereas for fMRI, it is mapped using the Balloon-Windkessel model (see subsubsection 2.1.7).

First, we demonstrate the inference on the whole-brain network of Wilson-Cowan model which generates beta oscillations characteristic of Parkinson’s disease [87, 88]. Figure 3 demonstrates the inference on the synaptic weights in an adapted Wilson-Cowan model (see subsubsection 2.1.1). Figure 3A and B show the observed and predicted neural activities (e.g., LFP recordings), and the corresponding power spectrum density (PSD) up to 20 Hz, respectively. Here, to estimate the set of unknown parameters , we conducted 10k random simulations for training of MAF density estimator. The non-informative prior was placed on all parameters (see Table 1), and the spatio-temporal features (such as statistical moments of time series) was used for training. Figure 3C demonstrates an accurate posterior estimation (in red), given the true values (in green). The model training was completed within a few minutes, and sampling from the posterior was nearly instantaneous. Figure 3D shows the evaluation of the inference by measuring the shrinkage and z-scores of the estimated posterior distributions. This result demonstrates an ideal Bayesian inference and prediction for this class of whole-brain network model. Note that removing the spatio-temporal features of time series and considering only the PSD as the data features leads to larger uncertainty in estimation (see Figure S2). Nevertheless, the observed and predicted PSD show close agreement.

Figure 3.

Then, we demonstrate the inference on heterogeneous control parameters in the whole-brain network of Jansen-Rit (see subsubsection 2.1.2), commonly used for modeling EEG/MEG data, e.g., in dementia and Alzheimer’s disease [138, 139]. Figure 4A and B show the observed and predicted EEG signals, given by (y_1i − y_2i) at each region, while Figure 4C and D illustrate the observed and predicted features such as PSD, respectively. Figure 4E and F show the estimated posterior distributions of synaptic connections C_i, and the global coupling parameter G, respectively, given the set of unknown parameters . Here we conducted 50k random simulations with samples drawn from uniform priors and (shown in blue, see Table 2). After approximately 45 min of training (MAF density estimator), the posterior sampling took only a few seconds. With such a sufficient number of simulations and informative data features, VBI shows accurate estimation of high-dimensional heterogeneous parameters (given the ground truth, shown in green), leading to a strong correspondence between the observed and predicted PSD of EEG/MEG data. Figure 4G displays the shrinkage and z-score as the evaluation metrics, indicating an ideal Bayesian estimation for C_i parameters, compared to the coupling parameter G. This occurred because the network input did not induce a significant change in the intrinsic frequency of activities at regional level, resulting in greater uncertainty in its estimation.

Figure 4.

Note that relying on only alpha-peak while excluding other summary statistics, such as total power (i.e., area under the curve), leads to poor estimation of synaptic connections across brain regions (see Figure S3). This results in less accurate predictions of the PSD, with more dispersion in their amplitudes. This example demonstrates that VBI provides a valuable tool for hypothesis evaluation and improved insight into data features by uncertainty quantification, and their impact on predictions.

To demonstrate efficient inference on the whole-brain time delay from EEG/MEG data, we used a whole-brain network model of coupled generic oscillators (StuartLandau model (see subsubsection 2.1.3). This model could establish a causal link between empirical spectral changes and the slower conduction velocities observed in multiple sclerosis patients, resulting from immune system attacks on the myelin sheath [49, 50]. The parameter set to estimate is , consisting of the global scaling parameter G and the averaged velocity of signal transmission V. The training was performed using a budget of only 2k simulations, which was sufficient due to the low dimensionality of the parameter space. Figure 5A illustrates the comparison between observed (in green) and predicted neural activities (in red). Figure 5B shows a close agreement between observed and predicted PSD signals, as the data feature used for training. Figure 5C and D provide visualizations of the posterior distributions for the averaged velocity V and the global coupling G. In these panels, we can see a large shrinkage in the posterior (in red) the uniform prior (in blue) centered around the true values (vertical green lines). Importantly, Figure 5E presenting the joint posterior distribution, indicates a high correlation of ρ = 0.7 between parameters G and V. This illustrates the advantage of Bayesian estimation in identifying statistical relationships between parameters, which helps to detect degeneracy among them. This is crucial for causal hypothesis evaluation and guiding conclusions in clinical settings. Finally, Figure 5F illustrates the sensitivity analysis (based on the eigenvalues of the posterior distribution), revealing that the posterior is more sensitive to changes in V compared to G. This highlights the relative impact of these parameters on the model’s posterior estimates.

Figure 5.

Next, we demonstrate the inference on a whole-brain model of epilepsy spread, known as the Virtual Epileptic Patient (VEP; [47, 76]), used to delineate the epileptogenic and propagation zone networks from the invasive sEEG recordings (subsubsection 2.1.4). Here, we used a large value for system time constant τ = 90 ms (see Table 4) to generate slow-fast dynamics in pathological areas, corresponding to seizure envelope at each brain region. Figure 6 demonstrates the inference the set of inferred parameters , with the global scaling parameter G and spatial map of epileptogenicity η_i, given i ∈ {1, 2,…, N_n = 88} parcelled regions. Figure 6A and B show the observed and predicted envelope, respectively, at each brain region. Here, the whole brain regions are classified into two epileptogenic zones (in red, corresponding to high excitability), three propagation zones (in yellow, corresponding to excitability close to bifurcation), and the rest as healthy regions (in green, corresponding to low excitability). Figure 6C and D illustrate the observed and predicted data features as the total power energy per region, calculated as the area under the curve. Additionally the seizure onset at each region was used as a data feature for training the MAF density estimator. From these panels, we observe accurate recovery of seizure envelopes in pathological regions. Figure 6E and F show that the posterior distribution of heterogeneous η_i, and global coupling parameter G, respectively, indicating 100% accurate recovery of the true values (in green). Figure 6G confirms the reliability and accuracy of the estimates through shrinkage and z-score diagnostics. With our efficient implementation, generating 10k whole-brain simulations took less than a minute (using 10 CPU cores). The training took approximately 13 minutes to converge, while posterior sampling required only a few seconds. See Figure S4 for a similar analysis with a faster time separation (τ = 10 ms). These results demonstrate an ideal and fast Bayesian estimation, despite the stiffness of equations in each region and the high dimensionality of the parameters.

Figure 6.

Targeting the fMRI data, we demonstrate the inference on the whole-brain dynamics using Montbrió model (see subsubsection 2.1.5). Figure 7 demonstrates the inference on heterogeneous control parameters of the Montbrió model, operating in a bistable regime (Table 5). Figure 7A and B show the observed and predicted BOLD time series, respectively, while Figure 7C and D illustrate the observed and predicted data features, such as the static and dynamical functional connectivity matrices (FC and FCD, respectively). Figure 7E and F show the estimated posterior distributions of excitability η_i per brain region, and the global coupling parameter G. Figure 7G displays the reliability and accuracy of estimation through the evaluation of posterior shrinkage and z-score (see Equation 12 and 13).

Figure 7.

Due to the large number of simulations for training and the informativeness of the data features (both spatio-temporal and functional data features), the results indicate that we achieved accurate parameter estimation, consequently, a close agreement between the observed and predicted features of BOLD data. This required 500k simulations for training (see Figure S6), given the uniform priors and . After approximately 10 h of training (of MAF density estimator), posterior sampling took only 1 min. Our results indicate that training the MAF model was 2 to 4 times faster than the NSF model (see Figure S1B). Note that removing the spatio-temporal features and considering only FC/FCD as the data features (see Figure S7) leads to poor estimation of the excitability parameter across brain regions (see Figure S5). Interestingly, accurate estimation of the only global coupling parameter, , from only FC/FCD requires around 100 simulations (see Figure S8). This showcase demonstrates the capability of VBI in inferring heterogeneous excitability, given the bistable brain dynamics, for fMRI studies.

Finally, in Figure 8, we show the inference on the so-called pDMF model, i.e., a whole-brain netwrok model of the reduced Wong-Wang equation (see subsubsection 2.1.6), comprising 10 control parameters: the global scaling of connections G and the linear coefficients (a_w, b_w, c_w, a_I, b_I, c_I, a_σ, b_σ, c_σ) ∈ ℝ⁹. These parameters are introduced to reduce the dimension of whole-brain parameters as recurrent connection strength w_i, external input current I_i, and noise amplitude σ_i for each region (in total, 264 parameters were reduced to 9 dimensions; see Equation 7). Here, we used summary statistics of both spatio-temporal and functional data features extracted from simulated BOLD data to train the MAF density estimator, with a budget of 50k simulations. The training took around 160 min to converge, whereas posterior sampling took only a few seconds.

Figure 8.

The diagonal panels in Figure 8A show estimated posterior distributions (in red), along with the prior (in blue) and true values (in green). The upper diagonal panels illustrate the joint posterior distributions between parameters (i.e., statistical dependency between parameters). Figure 8B illustrates the observed and predicted BOLD time series, generated by true and estimated parameters (in blue and red, respectively). From Figure 8C, we can see a close agreement between the observed and predicted data features (FC/FCD matrices). Note that due to the stochastic nature of the generative process, we do not expect an exact element-wise correspondence between these features, but rather a match in their summary statistics, such as the mean, variance, and higher-order moments (see Figure S9). Figure 8D shows the posterior z-score versus shrinkage, indicating less accurate estimation for the coefficients c_w, c_I, and c_σ compared to others, as they are not informed by anatomical data such as the T1w/T2w myelin map and the first FC principal gradient (see Equation 7). This showcase demonstrates the advantage of Bayesian inference over optimization in assessing the accuracy and reliability of parameter estimation, whether informed by anatomical data.

Note that in the whole-brain network of Wong-Wang model (6), the global scaling parameter G and synaptic coupling J exhibit structural non-identifiability, meaning their combined effects on the system cannot be uniquely disentangled. This is evident in the parameter estimations corresponding to selected observations, where the posterior distributions appear diffuse. The joint posterior plots reveal a nonlinear dependency (banana shape) between G and J, arising from their product in the neural mass equation (see Equation 6). Such a nonlinear relationship between parameters poses challenges for deriving causal conclusions, as often occurs in other neural mass models. This is a demonstration of how Bayesian inference facilitates causal hypothesis testing without requiring additional non-identifiability analysis.

4. Discussion

This study introduces the Virtual Brain Inference (VBI), a flexible and integrative toolkit designed to facilitate probabilistic inference on complex whole-brain dynamics using connectome-based models (forward problem) and simulation-based inference (inverse problem). The toolkit leverages high-performance programming languages (C++) and dynamic compilers (such as Python’s JIT compiler), alongside the computational power of parallel processors (GPUs), to significantly enhance the speed and efficiency of simulations. Additionally, VBI integrates popular feature extraction libraries with parallel multiprocessing to efficiently convert simulated time series into low-dimensional summary statistics. Moreover, VBI incorporates state-of-the-art deep neural density estimators (such as MAF and NSF generative models) to estimate the posterior density of control parameters within whole-brain models given low-dimensional data features. Our results demonstrated the versatility and efficacy of the VBI toolkit across commonly used whole-brain network models, such as Wilson-Cowan, Jansen-Rit, Stuart-Landau, Epileptor, Montbrió, and Wong-Wang equations placed at each region. The ability to perform parallel and rapid simulations, coupled with a taxonomy of feature extraction, allows for detailed and accurate parameter estimation from associated neuroimaging modalities such as (s)EEG/MEG/fMRI. This is crucial for advancing our understanding of brain dynamics and the underlying mechanisms of various brain disorders. Overall, VBI represents a substantial improvement over alternative methods, offering a robust framework for both simulation and parameter estimation, and contributing to the advancement of network neuroscience, potentially, to precision medicine.

The alternatives for parameter estimation include optimization techniques [66, 34, 36, 63, 64], ABC method, and MCMC sampling. Optimization techniques are sensitive to the choice of the objective function (e.g., minimizing distance error or maximizing correlation) and do not provide estimates of uncertainty. Although multiple runs and thresholding can be used to address these issues, such methods often fall short in revealing relationships between parameters, such as identifying degeneracy, which is crucial for reliable causal inference. Alternatively, a technique known as ABC compares observed and simulated data using a distance measure based on summary statistics [140, 141, 142]. It is known that, ABC methods suffer from the curse of dimensionality, and their performance also depends critically on the tolerance level in the accepted/rejected setting [71, 70]. The self-tuning variants of MCMC sampling have also been used for model inversion at the whole-brain level [76, 67]. Although MCMC is unbiased and exact with infinite runs, it can be computationally prohibitive, and sophisticated reparameterization methods are often required to facilitate convergence at whole-brain level [123, 143]. This becomes more challenging for gradient-based MCMC algorithms, due to the bistability and stiffness of neural mass models. Tailored to Bayes’ rule, SBI sidesteps these issues by relying on expressive deep neural density estimators (such as MAF and NSF) on low-dimensional data features to efficiently approximate the posterior distribution of model parameters. Taking spiking neurons as generative models, this approach has demonstrated superior performance compared to alternative methods, as it does not require model or data features to be differentiable [71, 81].

In previous studies, we demonstrated the effectiveness of SBI on virtual brain models of neurological [68, 49, 50, 54], and neurodegenerative diseases [69, 51, 23] as well as focal intervention [59] and healthy aging [37]. In this work, we extended this probabilistic methodology to encompass a broader range of whole-brain network models, highlighting its flexibility and scalability in leveraging diverse computational resources, from CPUs/GPUs to high-performance computing facilities. Our results indicated that the VBI toolkit effectively estimates posterior distributions of control parameters in whole-brain network modeling, offering a deeper understanding of the mechanisms underlying brain activity. For example, using the Montbrio and Wong-Wang models, we achieved a close match between observed and predicted FC/FCD matrices derived from BOLD time series (Figure 7 and Figure 8). Additionally, the Jansen-Rit and Stuart-Landau models provided accurate inferences of PSD from neural activity (Figure 4, and Figure 5), while the Epileptor model precisely captured the spread of of seizure envelopes (Figure 6). These results underscore the toolkit’s capability to manage complex, high-dimensional data with precision. Uncertainty quantification using VBI can illuminate and combine the informativeness of data features (e.g., FC/FCD) and reveal the causal drivers behind interventions [69, 59]. This adaptability ensures that VBI can be applied across various neuroimaging modalities, accommodating different computational capabilities and research needs.

Since data features are key in SBI, some factors can lead to the collapse of this method. For instance, high noise levels in observations or dynamical noise can compromise the accurate calculation of data features, undermining the inference process. Identifying the set of low-dimensional data features that are relevant to the control parameters for each case study is another challenge in effectively applying SBI. Nevertheless, the uncertainty of the posterior informs us about the predictive power of these features. Statistical moments of time series could be effective candidates for most models. However, this poses a formidable challenge for inference from empirical data, as certain moments, such as the mean and variance, may be lost during preprocessing steps. The hyperparameter and noise estimation can also be challenging for SBI. Moreover, there is no established rule for determining the number of simulations for training, aside from relying on z-score values during in-silico testing, as it depends on available computational resources.

Various sequential methods, such as SNPE [125], SNLE [144], and SNRE [145], have been proposed to reduce computational costs of SBI by iteratively refining the fit to specific targets. These approaches aim for more precise parameter estimation by progressively adjusting the model based on each new data set or subset, potentially enhancing the accuracy of the fit at the reduced computational effort. The choice of method depends on the specific characteristics and requirements of the problem being addressed [146]. Our previous study indicates that for inferring whole-brain dynamics of epilepsy spread, the SNPE method outperforms alternative approaches [68]. Nevertheless, sequential methods can become unstable, with simulators potentially diverging and causing probability mass to leak into regions that lack prior support [68]. In this study, we used single-round training to benefit from an amortization strategy. This approach brings the costs of simulation and network training upfront, enabling inference on new data to be performed rapidly (within seconds). This strategy facilitates personalized inference at the subject level, as the generative model is tailored by the SC matrix, thereby allowing for rapid hypothesis evaluation specific to each subject (e.g., in delineating the epileptogenic and propagation zones). Note that model comparison across different configurations or model structures, as well-established in dynamic causal modeling [147, 148], has yet to be explored in this context.

Deep learning algorithms are increasingly gaining traction in the context of whole-brain modeling. The VBI toolkit leverages a class of deep generative models, called Normalizing Flows (NFs; [127, 128]), to model probability distributions given samples drawn from those distributions. Using NFs, a base probability distribution (e.g., a standard normal) is transformed into any complex distribution (potentially multi-modal) through a sequence of invertible transformations. Variational autoencoders (VAEs; [149, 150]) is a class of deep generative models to encode data into a latent space and then decode it back to reconstruct the original data. Recently, Sip et al. [151] introduced a method using VAEs for nonlinear dynamical system identification, enabling the inference of neural mass models and region- and subject-specific parameters from functional data. VAEs have also been employed for dimensionality reduction of whole-brain functional connectivity [42], and to investigate various pathologies and their severity by analyzing the evolution of trajectories within a low-dimensional latent space [52]. Additionally, Generative Adversarial Networks (GANs; [152, 153]) have demonstrated remarkable success in mapping latent space to data space by learning a manifold induced from a base density [154]. This method merits further exploration within the context of whole-brain dynamics. To fully harness the potential of deep generative models in large-scale brain network modeling, integrating VAEs and GANs into the VBI framework would be beneficial. This will elucidate their strengths and limitations within this context and guide future advancements in the field.

In summary, VBI offers fast simulations, taxonomy of feature extraction, and deep generative models, making it a versatile tool for model-based inference from different neuroimaging modalities, helping researchers to explore brain (dys)functioning in greater depth. This advancement not only enhances our theoretical understanding but also holds promise for practical applications in diagnosing and treating neurological conditions.

Information Sharing Statement

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

Glossary of technical terms

Bayesian Rule: A fundamental belief updating principle that calculates the probability of a hypothesis given new evidence.

Deep neural density estimators: A class of artificial neural network-based approaches that are used to learn and approximate the underlying probability distribution from a given dataset.

Control (generative) parameters: The bifurcation parameters, setting, or configuration within a generative model that controls the synthesis of data and potentially represents causal relationships.

Generative model: A statistical, machine learning, or mechanistic model that represents the underlying data distribution to generate new data resembling the original dataset.

Likelihood: The conditional probability of observing the evidence given a particular hypothesis.

Markov chain Monte Carlo: A family of stochastic algorithms used for uncertainty quantification by drawing random samples from probability distributions, in which the sampling process does not require knowledge of the entire distribution.

Prior: The initial probability assigned to a hypothesis before considering new evidence.

Probability distribution: The statistical description of potential outcomes of random events, where a numerical measure is assigned to the possibility of each specific outcome.

Probability density estimation: The process of inferring the underlying probability distribution of a random event based on observed data.

Posterior: The updated probability of a hypothesis after taking into account both prior beliefs and observed evidence.

Simulation-based inference: A machine learning approach that involves generating synthetic data through forward simulations to make inferences about complex systems, often when analytic or computational solutions are unavailable.

Whole-brain or virtual brain models: Computational models informed by personalized anatomical data, i.e., a set of equations describing regional brain dynamics placed at each node, which are then connected through structural connectivity matrix.

6. Supplementary figures

Figure S1.

Figure S2.

Figure S3.

Figure S4.

Figure S5.

Figure S6.

Figure S7.

Figure S8.

Figure S9.

Figure S10.

Figure S11.

Figure S12.

Acknowledgements

This research has received funding from EU’s Horizon 2020 Framework Programme for Research and Innovation under the Specific Grant Agreements No. 101147319 (EBRAINS 2.0 Project), No. 101137289 (Virtual Brain Twin Project), No. 101057429 (project environMENTAL), and government grant managed by the Agence Nationale de la Recherche reference ANR-22-PESN-0012 (France 2030 program). We acknowledge the use of Fenix Infrastructure resources, which are partially funded from the European Union’s Horizon 2020 research and innovation programme through the ICEI project under the grant agreement No. 800858. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Additional information

Author contributions

Conceptualization: V.J., and M.H. Methodology: A.Z., M.W., L.D., and M.H. Software: A.Z, M.W., L.D., and M.H. Investigation: A.Z. Visualization: A.Z., Supervision: V.J., S.P., and M.H. Funding acquisition: V.J., M.H. Writing - original draft: A.Z., and M.H. Writing - review & editing: A.Z, M.W., L.D., S.P., V.J, and M.H.

Abbreviations

• VBI: virtual brain inference

• BOLD: blood-oxygen-level-dependent

• fMRI: functional magnetic resonance imaging

• EEG: Electroencephalography

• MEG: Magnetoencephalography

• sEEG: Stereoelectroencephalography

• SC: structural connectivity

• FC: functional connectivity

• FCD: functional connectivity dynamic

• PSD: power spectral density

• SBI: simulation-based inference

• MAF: masked autoregressive flow

• NSF: neural spline flow

• MCMC: Markov chain Monte Carlo