This probabilistic approach is designed to estimate the posterior distribution of control parameters in virtual brain models from whole-brain neuroimaging recordings. (A) The process begins with constructing a personalized connectome using diffusion tensor imaging and a brain parcellation atlas. (B) The personalized virtual brain model is then assembled. Neural mass models describing the averaged activity of neural populations, in the generic form of , are placed to each brain region, and interconnected via the structural connectivity matrix. Initially, the control parameters are randomly drawn from a simple prior distribution. (C) Next, the VBI operates as a simulator that uses these samples to generate time series data associated with neuroimaging recordings. (D) We extract a set summary statistics from the low-dimensional features of the simulations (FC, FCD, PSD) for training. (E) Subsequently, a class of deep neural density estimators is trained on pairs of random parameters and their corresponding data features to learn the joint posterior distribution of the model parameters. (F) Finally, the amortized network allows us to quickly approximate the posterior distribution for new (empirical) data, enabling us to make probabilistic predictions that are consistent with the observed data.
Parameter descriptions for capturing whole-brain dynamics using Wilson-Cowan neural mass model.
Parameter descriptions for capturing whole-brain dynamics using Jansen-Rit neural mass model. EP: excitatory populations, IP: inhibitory populations, PSP: post synaptic potential, PSPA: post synaptic potential amplitude.
Parameter descriptions for capturing whole-brain dynamics using Stuart-Landau oscillator.
Parameter descriptions for capturing whole-brain dynamics using 2D Epileptor neural mass model.
Parameter descriptions for capturing whole-brain dynamics using Montbrió model.
Parameter descriptions for capturing whole-brain dynamics using Wong-Wang model.
Parameter descriptions for the Balloon-Windkessel model to map neural activity to the BOLD signals detected in fMRI.
Flowchart of the VBI Structure.
This toolkit consists of three main modules: (1) The Simulation module, implementing various whole-brain models, such as Wilson-Cowan (WCo), Jansen-Rit (JR), Stuart-Landau (SL), Epileptor (EPi), Montbrió (MPR), and Wong-Wang (WW), across different numerical computing libraries (C++, Cupy, PyTorch, Numba). (2) The Features module, offering an extensive toolbox for extracting low-dimensional data features, such as spectral, temporal, connectivity, statistical, and information theory features. (3) The Inference module, providing neural density estimators (such as MAF and NSF) to approximate the posterior of parameters.
Bayesian inference on the control parameters of the whole-brain network model of Wilson-Cowan, generating the beta oscillation.
The inferred parameters are , where subscripts S, G, and C denote the subthalamic nucleus, globus pallidus, and cortex, respectively. We conducted 10k random simulations for training of MAF density estimator. (A) The observed and predicted neural activities. (B) The observed and predicted spectral power densities up to 20 Hz. (C) The estimated posterior distributions (in red), along with the prior (in blue) and ground-truth values (in green). (D) The posterior shrinkage versus z-scores for measuring the accuracy and reliability of the estimates.
Bayesian inference on heterogeneous control parameters in the whole-brain network of Jansen-Rit model.
The set of inferred parameters is , with the global scaling parameter G and average numbers of synapse between neural populations per Ci region, given i ∈ {1, 2,…, Nn = 88} parcelled regions. Summary statistics of power spectrum density (PSD) were used for training, with a budget of 50k simulations. (A) and (B) illustrate the observed and predicted neural activities, respectively. (C) and (D) show the observed and predicted data features, such as PSD. (E) and (F) display the posterior distribution of Ci per region, and global coupling G, respectively. The ground truth and prior are represented by vertical green lines and a blue distribution, respectively. (G) shows the inference evaluation the shrinkage and z-score of the estimated posterior distributions.
Bayesian inference on global scaling parameter G and the averaged velocity V of signal transmission using the whole-brain network model of Stuart-Landau oscillators.
The set of estimated parameters is , and the summary statistics of PSD signals with a budget of 2k simulations were used for training. (A) illustrates exemplary observed and predicted neural activities (in green and red, respectively). (B) shows the observed and predicted PSD signals (in green and red, respectively). (C) and (D) display the posterior distribution of averaged velocity V and global coupling G, respectively. The true values and prior are shown as vertical green lines and a blue distribution, respectively. (E) shows the joint posterior distribution indicating a high correlation between posterior samples. (F) illustrates the sensitivity analysis based on the eigenvalues of the posterior distribution.
Bayesian inference on the spatial map of epileptogenicity across brain regions in the VEP model.
The set of inferred parameters is , as the global scaling parameter and spatial map of epileptogenicity with i ∈ {1, 2,…, Nn = 88} parcelled regions. (A) The observed seizure envelope generated by the Epileptor model, given two regions as epileptogenic zones (in red) and three regions as propagation zones (in yellow), while the rest are healthy (in green). (B) The predicted seizure envelope, by training MAF model on a dataset containing 10k simulations, using only the total power and seizure onset per region as the data features. (C) and (D) show the observed and predicted data features, respectively. (E) and (F) show the posterior distributions of heterogeneous control parameters ηi, and global coupling parameter G, respectively. (G) The posterior z-scores versus posterior shrinkages for estimated parameters.
Bayesian inference on heterogeneous control parameters in the whole-brain dynamics using Montbrió model.
The set of inferred parameters is , as the global scaling parameter and excitability per region, with i ∈ {1, 2,…, Nn = 88} parcelled regions. VBI provides accurate and reliable posterior estimation using both spatio-temporal and functional data features for training, with a budget of 500k simulations. (A) and (B) illustrate the observed and predicted BOLD signals, respectively. (C) and (D) show the observed (upper triangular) and predicted (lower triangular) data features (FC and FCD), respectively. (E) and (F) display the posterior distribution of excitability parameters ηi per region, and global coupling G, respectively. The true values and prior are shown as vertical green lines and a blue distribution, respectively. (G) shows the inference evaluation by the shrinkage and z-score of the posterior distributions.
Bayesian inference on the parametric mean-field model of Wong-Wang (pDMF model), with linear coefficients (aw, bw, cw, aI, bI, cI, aσ, bσ, cσ) ∈ ℝ9, reparameterizing the recurrent connection strength wi, external input current Ii, and noise amplitude σi for each region.
Summary statistics of spatio-temporal and functional data features were used for training, with a budget of 50k simulations. (A) The diagonal panels display the ground-true values (in green), the uniform prior (in blue), and the estimated posterior distributions (in red). The upper diagonal panels illustrate the joint posterior distributions between parameters, along with their correlation (ρ, in the upper left corners), and ground-truth values (green stars). High-probability areas are color-coded in yellow, while low-probability areas are shown in black. (B) The observed and predicted BOLD time series (in green and red, respectively). (C) The observed and predicted data features, such as FC/FCD matrices. (D) The inference evaluation by calculating the shrinkage and z-score of the estimated posterior distributions.
(A) Benchmarking the simulation cost of the whole-brain network of Montbrió model, with 2.5M time step on Intel(R) Core(TM) i9-10900 CPU 2.80 GHz (in red) and NVIDIA RTX A5000 GPU (in blue). GPUs deliver substantial speedups up to 100x over multi-core CPUs. (B) Comparing MAF and NSF density estimators, MAF is typically 2-4 times faster than NSF during the training process.
Inference on generative parameters in the Wilson-Cowan model (Equation 1) of beta oscillation.
The set of inferred parameters is θ = wSG, wGS, wCS, wSC, wGG, wCC ∈ ℛ6. (A) The diagonal panels show the estimated posterior distributions (in red), true values (in green), and the uniform prior (in blue). The upper diagonal panels depict the joint posterior distributions between parameters, including their correlations (ρ, in the upper left corners), with true values marked by green stars. High-probability areas are color-coded in yellow, while low-probability areas are shown in black. Panels (B), (C) display the observed and predicted activity timeseries, respectively. Panels (D), (E) show the observed and predicted PSD, respectively, used as data features for training. (F) Evaluates the inference by presenting the shrinkage and z-scores of the estimated posterior distributions.
Inference on heterogeneous generative parameters in the whole-brain network of Jansen-Rit model (Equation 2).
The set of inferred parameters is , consisting of the global scaling parameter G and the average number of synapses between neural populations Ci per region, with i ∈ {1, 2,…, Nn = 88} parcelled regions. For panels A to C, only the area under the curve of the power spectral density (PSD) in the [0,30] Hz band was used, while for panels D to F, additional statistical features of the PSD were included in the training, with a total budget of 50k simulations. (A) and (D) illustrate the observed and predicted EEG signals. (B) and (E) show the observed and predicted PSD signals. (C) and (F) display the posterior distribution of Ci for each region. The true values are indicated by vertical green lines, and the prior distribution is shown in blue.
Bayesian inference on the spatial map of epileptogenicity across brain regions in the VEP model (Equation 4), with slow time scale seperation, τ = 10 ms.
(A) The observed seizure envelope generated by the 2D Epileptor model given two regions as epileptogenic zones (in red) and three regions as propagation zones (in yellow), while the rest are healthy (in green). (B) The predicted seizure envelope, by training MAF density estimator on a dataset containing 10k simulations, using only the total power per region as the data features. (C) and (D) show the observed and predicted data features, respectively. (E) and (F) show the posterior distribution of heterogeneous ηi, and global coupling parameter G, respectively. (F) The posterior z-scores versus posterior shrinkages for estimated parameters.
Inference on heterogeneous generative parameters in the whole-brain network of Montbrió model (Equation 5) involves estimating a parameter set denoted as .
Here, G is the global scaling parameter, and ηi represents the excitability of each brain region, where i ∈ {1, 2,…, Nn = 88} corresponds to the parcelled regions. Relying on only functional data features (FC/FCD) for training proves insufficient to accurately estimate the posterior distribution of heterogeneous excitabilities, even with a budget of 500k simulations. In (A), the posterior estimates are shown in red, the true values of ηi in green, and the prior distribution in blue. Panels (B) and (C) show the observed and predicted FC/FCD matrices, respectively. The simulation to generate training data took approximately 24 hours, utilizing 10 GPUs concurrently. The training process lasted around 8 hours, while sampling required just a few seconds.
Evaluating the reliability and accuracy of estimation (parameter G) with respect to the number of simulation used for training the MAF/NSF density estimator by (A) posterior shrinkage, and (B) z-score versus the number of simulations.
Exemplary summary statistics from the FC/FCD matrices in the whole-brain network of Montbrió model (Equation 5).
Inference of the global coupling parameter in the whole-brain network of Montbrió model (Equation 5) using functional data features, using a budget of only 100 simulations.
(A) The uniform prior distribution of G in blue, the estimated posterior in red and the true value is shown in green. (B) and (C) The observed and predicted FC/FCD matrices, respectively.
The distribution of observed and predicted functional connectivity (FC) in panel and the functional connectivity dynamics in panel (B), which are demonstrated in Figure 8.
Inference of the global coupling parameter in the whole-brain model of Montbrió (Equation 5) using TVB simulator, using summation of FC and FCD matrix elements as low dimensional data features, constrained by a budget of only 100 simulations.
(A) The uniform prior distribution of G in blue, the estimated posterior in red and the true value is shown in green. (B) and (C) show the observed and predicted FC/FCD matrices, respectively.
The non-identifiability in estimating G and J in the Wong-Wang model (Equation 6).
A total of 16k simulations and functional data features (see Figure S11) were used for training. The upper row illustrates the FCD variance and the sum of FC matrices plotted against the parameters G and J, respectively, for three selected observations marked by red circles. The lower panels display the parameter estimations corresponding to the selected observations. Due to structural non-identifiability in the product of G and J, the posterior distributions appear diffuse, showing a nonlinear dependency (banana shape) in the joint posterior plots.
Exemplary summary statistics from the FC/FCD matrices in the whole-brain network of Wong-Wang model (Equation 6).
