Simulation-based inference on virtual brain models of disorders

The VBMs emphasize the structural network characteristics of the brain by representing the brain regions as nodes, which are connected via a SC matrix representing white matter fibre tracts (Sanz-Leon et al 2015). The $\mathrm{SC} \in \mathbb{R}^{\mathit{N}_{\mathit{n}} \times \mathit{N}_{\mathit{n}}}$ with N_n parcelled regions used in this study was taken from the open-source brain network modeling tool, the virtual brain (TVB; Sanz Leon et al 2013, Schirner et al 2015). The T1-weighted MRI images were processed to obtain the brain parcellation, whereas diffusion-weighted (DW-MRI) images were used for tractography. With the generated fiber tracts and with the regions defined by the brain parcellation, the connectome was built by counting the fibers connecting all regions. Using Desikan–Killiany parcellation (Desikan et al 2006) in the reconstruction pipeline, the patient’s brain was divided into $N_n = 88$ cortical and subcortical regions. The connectome was normalized so that the maximum value is equal to one.

To ensure a realistic in-silico experiment and therefore assess the function-structure link during the brain disorders via VBMs, we simulated the BOLD fMRI using a whole-brain network model that couples exact mean-field representation of spiking neurons through the weighted edges in the SC matrix. Assuming a Lorentzian distribution on the membrane potentials across decoupled brain regions, with half-width Δ centered at η, the regional brain dynamics can be described by a neural mass model derived analytically in the limit of infinitely all-to-all coupled quadratic integrate-and-fire neurons (Montbrió et al 2015).

By coupling the brain regions via an additive current in the average membrane potential equations, the dynamics of the whole-brain network can be described as follows (Rabuffo et al 2021, Fousek et al 2022):

$$\begin{align} \tau\dot{r_i}(t) & = 2 r_i(t) v_i(t) + \dfrac{\Delta}{\pi \tau} \nonumber \\ \tau \dot{v_i}(t) & = v_i^2(t) - (\pi \tau r_i(t))^2 + J \tau r_i(t) + \eta + I_{stim}(t) + G \sum_{j = 1}^{N_n} SC_{ij}(r_{j}(t)-r_{i}(t))) + \xi(t), \end{align} \tag{ 1 }$$

where v_i and r_i are the average membrane potential and firing rate, respectively, at the ith brain region. The excitability η, the synaptic weight J, the spread of the heterogeneous distribution Δ, and the characteristic time constant τ, can be tuned so that each decoupled node is in a bistable regime, exhibiting a down-state stable fixed point and an up-state stable focus in the 2D phase space (Montbrió et al 2015). 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 (Rabuffo et al 2021). The input current $I_\textrm{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 (Rabuffo et al 2021). Moreover, the parameter G is the network scaling parameter that modulates the overall impact of brain connectivity on the state dynamics (see supplementary, figure S1), $\textrm{SC}_{i,j}$ denotes the symmetrical connection weight between ith and jth nodes with ${i,j} \in \{1, 2, \ldots, N_n \}$, and the dynamical noise $\xi(t) \sim \mathcal N(0, {\sigma}^2)$ follows a Gaussian distribution with mean zero and variance $\sigma^2 = 0.03$.

In this study, the brain network model was implemented in TVB software (Sanz Leon et al 2013, Rabuffo et al 2021) and equipped with BOLD forward solution comprising the Balloon-Windkessel model (Stephan et al 2007) applied to the firing rate. In addition, we have implemented the model in C++, which is around 2 times faster than our Python implementation with Just-in-Time (JIT) compilation, and also CUDA-based GPU, which provides up to 100 times faster computational cost through parallelization (see supplementary, figure S2). The model simulation and parameter estimation were performed on a Linux machine with Intel Core i0-10900 CPU 2.80 GHz and 32 GB RAM.

The approach of virtually modeling the brain connectivity in disorders provides the basis to investigate whether a specific modification in connectome affects the observed brain function in a causal sense. We specifically tested whether the disorder is related to inter-hemispheric and intra-hemispheric connections in SC, by applying two spatial masks on the subject-specific connectome (normalized by the scaling factor G) as follows:

$$\begin{align} \textrm{SC}_\textrm{disordered} = G \textrm{SC}_\textrm{healthy} - \alpha M_\textrm{inter} - \beta M_\textrm{intra}\end{align} \tag{ 2 }$$

where parameters α and β are the unknown normalized intensity of degradation in SC (as the target of parameter estimation in the range [0–1]), according to the spatial masks $M_\textrm{inter}$ and $M_\textrm{intra}$, on inter-hemispheric and intra-hemispheric connections, respectively. In this study, the $M_\textrm{intra}$ mask indicates the afferent and efferent connections originating from the limbic system (see figure 1). The motivation for these masks comes from the fact that they have been shown to capture the leading mechanism of healthy aging (Lavanga et al 2023) and Alzheimer’s Disease (Zimmermann et al 2018, Arbabyazd et al 2023).

Zoom In Zoom Out Reset image size

To quantify the communication between different brain regions and associated dynamics, we computed both static FC and FCD, which is the time-variant representation of FC and reflects the fluctuation of the covariance matrix over time (see supplementary, figure S1). From the simulated BOLD signals generated by TVB, we extracted FC matrix by calculating the Pearson correlation between the BOLD signals of any two brain regions. In order to track the time-dependent changes in the FC, we computed the windowed FCD, with length of sliding window τ = 30 s and step size of $t_{w+1}-t_{w} = 6$ s (Hansen et al 2015, Arbabyazd et al 2020):

$$\begin{align} \textrm{FCD}\left(w_i, w_j\right) = \textrm{corr}\left(\textrm{uppertri}\left(\textrm{FC}\left(w_i\right)\right), \textrm{uppertri}\left(\textrm{FC}\left(w_j\right)\right)\right)\end{align} \tag{ 3 }$$

where, the (i, j) entry of the FCD matrix provided the Pearson correlation between the upper triangular parts of the two FC matrices $\textrm{FC}(w_i)$ and $\textrm{FC}(w_j)$ calculated at each sliding windows $w = 1, 2, \ldots, N_w$. To assess whether the simulated BOLD time-series might have state transitions between up and down states, we quantified the brain’s fluidity by calculating the FCD variance to capture the transitioning behavior (see supplementary, figure S1).

The Bayesian approach offers a principled way for making inference, prediction, and quantifying uncertainty for decision-making process (Gelman et al 1995, Bishop 2006, Box and Tiao 2011). This approach naturally evaluates and incorporates uncertainties in the parameters and observations to drive meaningful conclusions from the data, with a broad range of applications (Bonomi et al 2016, Hashemi et al 2018, Khalvati et al 2019, Guimerà et al 2020, Broderick et al 2023, Sip et al 2023, Wang et al 2023). Parameter estimation within a Bayesian framework is treated as the quantification and propagation of uncertainty through probability distributions placed on the parameters, updated with the information provided by data (Hashemi et al 2020, 2021). Considering the joint probability distribution of the observation x and the unknown control parameters θ, referred to as the generative model. Then by the product rule: $p(\theta, x) = p(\theta \mid x) p(x)$, equivalently, $p(\theta, x) = p(x \mid \theta)p(\theta)$. Hence, given data x and model parameters θ, Bayes rule defines the posterior distribution as

$$\begin{align} p\left(\theta \mid x\right) = \frac{p\left(\theta\right) p\left(x \mid \theta\right)}{p\left(x\right)}\end{align} \tag{ 4 }$$

that combines and actualizes prior information (domain expertise) about unknown parameters with the knowledge acquired from observed data through the so-called likelihood function (data generating process). The prior information $p(\theta)$ is typically determined before seeing the data through beliefs and previous evidence about possible values of the parameters. The likelihood function $p(x \mid \theta)$ represents the probability of some observed outcomes given a certain set of parameters (the information about the parameters provided by the observed data). The denominator $p(x) = \int p(x\mid \theta)p(\theta)d\theta$ represents the probability of the data and it is known as evidence or marginal likelihood. In the context of inference, this term amounts to simply a normalization factor. Note that using Bayesian inference in this study, we aim to estimate the entire posterior distribution of the unknown parameters, i.e. the uncertainty over a range of plausible parameter values that generates the data, rather than a single point estimate for data fitting (e.g. the maximum likelihood estimation; Hashemi et al 2018, or the maximum a posteriori; Friston et al 2014, Razi et al 2015, Vattikonda et al 2021).

The key challenge to perform an efficient Bayesian inference is the evaluation of likelihood function $p(x \mid \theta)$. This is typically intractable for high-dimensional models involving nonlinear latent variables (such as VBMs), as the likelihood evaluation requires an integral over all possible trajectories through the latent space that controls the generative process: $p(x \mid \theta) = \int p(x, z \mid \theta) \textrm{d}z$, where $p(x, z \mid \theta)$ is the joint probability density of data x and unmeasured (hidden) latent variables z, given parameters θ. In particular, when dealing with high-dimensional latent space at whole-brain scales, the computational cost of evaluating the likelihood function can become prohibitive. This makes likelihood-based approaches, such as MCMC sampling, inapplicable (Hashemi et al 2023). SBI or likelihood-free inference performs efficient Bayesian inference for complex models where the calculation of likelihood function is either analytically or computationally intractable (Cranmer et al 2020). Instead of direct sampling from distributions and explicit evaluation of likelihood function, SBI sidesteps this problem by employing deep ANNs to learn an invertible transformation between parameters and data (more precisely, between the features of a simulated dataset and parameters of a parameterized approximation of the posterior distribution).

Taking prior distribution $p(\vec{\theta})$ over the parameters of interest $\vec{\theta}$, a limited number of N simulations are generated for training step $\{(\vec{\theta}_i, \vec{x}_i)\}_{i = 1}^{N}$, where $\vec{\theta_i} \sim p(\vec{\theta})$ and $\vec{x}_i$ is the simulated data given model parameters $\vec{\theta}_i$. In other words, the training data set is an ensemble of N independent and identically distributed samples from the generative model $p(\vec{\theta}, \vec{x}) = p(\vec{\theta}) p(\vec{x} \mid \vec{\theta})$, that can be run in parallel. After the training step, we are able to efficiently estimate the approximated posterior $q_{\phi}(\vec{\theta} \mid \vec{x})$ with learnable parameters φ, so that for the observed data $\vec{x}_\textrm{obs}$: $q_{\phi}(\vec{\theta} \mid \vec{x}_\textrm{obs}) \simeq p(\vec{\theta} \mid \vec{x}_\textrm{obs})$. For more details, see Hashemi et al (2023). Note that due to the amortized strategy, for any new observed data $\vec{x}_\textrm{obs}$, we can readily approximate the true posterior $p(\vec{\theta} \mid \vec{x}_\textrm{obs})$ by a forward pass through the trained ANN (in the order of seconds).

A powerful technique in probabilistic machine learning for sampling and probability density estimation is NSF (Rezende and Mohamed 2015). NFs are a family of generative machine learning models that convert a base distribution into any complex target distribution, where both sampling and density estimation can be efficient and exact (Papamakarios et al 2019a, Kobyzev et al 2020). The aim of density estimation is to offer an accurate description of the underlying probability distribution of an observable data set, where the density itself is unknown. In particular, NFs embedded in the SBI approach, such as (sequential) neural posterior estimation (Lueckmann et al 2017, Gonçalves et al 2020), allow for the direct estimation of joint posterior distributions, and bypass the need for MCMC, while also potentially capturing degeneracy or multi-modalities. Recently advanced models such as MAF (Papamakarios et al 2017) and NSF (Durkan et al 2019) provide efficient and exact density evaluation and sampling from the joint distribution of high-dimensional random variables in a single neural-network pass. These generative models use deep neural networks to learn complex mappings between input data and their corresponding probability densities. They have achieved state-of-the-art performance with diverse applications, which efficiently represent rich structured and multi-modal posterior distributions, capturing complex dependencies and variations in the data distribution (Papamakarios et al 2019a, Kobyzev et al 2020). In this study, we use autoregressive flows (MAF and NSF models), which support invertible nonlinear transformations with low-cost computation (Kobyzev et al 2020). MAF uses affine transformations (Durkan et al 2019) with binary masks, eliminating the need for sequential recursion (Germain et al 2015). NSF leverages splines as a coupling function, enhancing the flexibility of the transformations while retaining exact invertibility (Papamakarios et al 2017). In this study, 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, all as implemented in the public sbi toolbox (Tejero-Cantero et al 2020). By training MAF/NSF on TVB simulations with random parameters, we are able to readily estimate the full posterior of parameters from low-dimensional data features, such as FC/FCD.

We can parameterize the SBI-VBMs in three different ways:

• Pooled structure assumes a common distribution across all individual units (here, the excitability), treating them collectively as a single unit with no variation in the sampling process. For simplicity, pooled (or homogeneous) modeling of excitability assumes a shared value across all brain regions (see supplementary, figure S3(A)). In this case, the set of generative parameters is low-dimensional as $\vec{\theta} = \{G, \alpha, \beta, \eta \} \in \mathbb{R}^{4}$.
• Unpooled structure treats individual units as being sampled independently, allowing for adaptation to diverse characteristics and behaviors within the dataset without assuming a common structure. Unpooled (or heterogeneous) modeling of excitability assumes distinct characteristics for each region without imposing a uniform structure at the whole-brain level (see supplementary, figure S3(B)). This approach leads to more precise and biologically plausible inferences about the brain’s (dys)functioning but significantly increases the dimensionality of the parameter space. In this case, the set of generative parameters becomes high-dimensional as $\vec{\theta} = \{G, \alpha, \beta, \eta_i \} \in \mathbb{R}^{N_n+3}$, with $i \in {1, 2, \ldots, N_n}$, and N_n being the number of brain regions.
• Hierarchical (or partially pooled) structure offers a middle ground between the simplicity of a pooled structure and the complexity of a fully unpooled structure. This balance allows for more flexibility and efficiency in Bayesian modeling while maintaining biological plausibility (see supplementary, figure S3(C)). In this case, the set of generative parameters is optimal as $\vec{\theta} = \{G, \alpha, \beta, \mu_{\eta}, \sigma_{\eta} \} \in \mathbb{R}^{5}$, and $\eta_i \sim \mathcal{N}(\mu_{\eta}, \sigma_{\eta}^{2})$.

Sensitivity analysis is a necessary step to determine which model parameters mostly contribute to variations in the model’s behavior due to changes in the model’s input, i.e. the identifiability analysis (Hashemi et al 2018, 2023). A local sensitivity coefficient measures the influence of small changes in one model parameter on the model output, while the other parameters are held constant. This can be quantified by computing the curvature of objective function through the Hessian matrix (Bates and Watts 1980, Hashemi et al 2018):

$$\begin{align} S\left(\vec{\theta}\right) = \mathcal{D} \left(\dfrac{\partial^2 \mathcal{E}\left({\vec{\theta}}\right)}{\partial \vec{\theta} \partial \vec{\theta}^{\top}} \mid_{{\vec{\theta}}^{\ast}}\right),\end{align} \tag{ 5 }$$

where $\mathcal{D}$ denotes the main diagonal elements of a matrix, $\mathcal{E}$ indicates a fitness function, centered at estimated parameters ${\vec{\theta}}^{\ast}$.

Using the SBI approach, after the training step and posterior estimation for a specific observation, we can also efficiently perform sensitivity analysis by calculating the eigenvalues and corresponding eigenvectors from (Tejero-Cantero et al 2020, Deistler et al 2022):

$$\begin{align} M = \mathbb{E}_{p\left(\vec{\theta}|\vec{x}_\textrm{obs}\right)}\left[\nabla_{\theta}p\left(\vec{\theta}|\vec{x}_\textrm{obs}\right)^T\nabla_{\theta}p\left(\vec{\theta}|\vec{x}_\textrm{obs}\right)\right],\end{align} \tag{ 6 }$$

which then does an eigendecomposition $M = Q \Lambda Q^{-1}$. A strong eigenvalue of the so-called active subspaces (Constantine 2015) indicates that the gradient of the posterior is large, hence, the system output is sensitive to changes along the direction of the corresponding eigenvector.

To measure the reliability of the Bayesian inference using synthetic data, we evaluate the posterior z-scores (denoted by z) against the posterior shrinkage (denoted by s), which are defined as (Betancourt 2014):

$$\begin{align} z & = | \dfrac{\bar \theta-\theta^\ast}{\sigma_\textrm{post}}|,\end{align} \tag{ 7 }$$

$$\begin{align} s & = 1- \dfrac{\sigma^2_\textrm{post}}{\sigma^2_\textrm{prior}},\end{align} \tag{ 8 }$$

where $\bar \theta$ and $\theta^\ast$ are the posterior mean and the true values, respectively, whereas σ_prior, and σ_post indicate the standard deviations of the prior and the posterior, respectively. The posterior z-score quantifies how much the posterior distribution of a parameter encompasses its true value. The posterior shrinkage quantifies how much the posterior distribution contracts from the initial prior distribution. The concentration of estimations towards large shrinkages indicates that all the posteriors are well-identified, while the concentration towards small z-scores indicates that the true values are accurately encompassed in the posteriors.

Figure 1 illustrates the overview of our approach, referred to as SBI-VBMs, to make probabilistic predictions on brain disorders by estimating the joint posterior distribution of the generative parameters. At the first step, the non-invasive brain imaging data such as T1-weighted MRI, and DW-MRI are collected for a specific patient (figure 1(A)). With the generated fiber tracts and the regions defined by the brain parcellation, the connectome (i.e. the total set of links between brain regions) is constructed by counting the fibers connecting all regions. The SC matrix, whose entries represent the connection strength between brain regions, provides the basis for constructing a network model–TVB–which is capable of generating various functional brain imaging data at arbitrary locations (e.g. cortical and subcortical structures). This step constitutes the structural brain network component, which imposes a constraint on brain network dynamics, allowing the hidden dynamics to be inferred from the data.

Then, each brain network node is assigned a computational model of average neuronal activity, a so-called neural mass model (see figure 1(B)). Here, we use an exact mean-field approximation for a population of spiking neurons, which due to its rich dynamics, such as bi-stability, enables us to generate various spatio-temporal dynamics as observed experimentally. Subsequently, the perturbation of SC by spatial masking, the inter- and intra-hemispheric edges can provide an estimation of the level of degradation in the patient connectome caused by disorders (see equation (2)). Notably, the stimulation of certain brain regions can be used to increase the causal evidence in the structure-function relationship. Embedding the disordered SC into the VBM, then complements the generative brain model of disorders (see equation (1)). The set of generative parameters θ comprises the regional excitability parameters η_i , a global coupling parameter G describing the scaling of the brain’s SC, and parameters α and β indicating the levels of interhemispheric deterioration, and intrahemispheric degradation within the limbic system, respectively. With an optimal set of these parameters, VBMs can mimic essential data features observed in fMRI recordings, such as FC/FCD matrices.

Finally, we use SBI (see figure 1(C)) to estimate the posterior distribution of the generative parameters in brain disorders. For this purpose, we first draw the parameters randomly from a prior distribution, i.e. $\vec{\theta}_i \sim p(\vec{\theta})$, for each ith simulation. The VBM simulator takes the parameters $\vec{\theta}$ as input and generates summary statistics $\vec{x}$ of the data as output. By preparing a training dataset through performing VBM simulations with random parameters, a class of deep neural density estimators (such as MAF or NSF model) is then trained on $\{(\vec{\theta}_i, \vec{x}_i)\}_{i = 1}^{N}$ with budget of N simulations, to learn an invertible transformation between data features and parameters of an approximated posterior distribution $p(\vec{\theta} \mid \vec{x})$ (see equation (4)). This is the primary advantage of this approach that the amortized inference at the subject level enables us to evaluate different hypotheses with uncertainty quantification at a negligible computational cost, typically in the order of a few seconds.

Here we show the estimation over the set of unknown generative parameters denoted by $\vec{\theta} = \{G, \alpha, \beta\} \in \mathbb{R}^{3}$. Figure 2 shows the estimated posterior and the sensitivity analysis, illustrating the impact of data features on the inference process. For the training (here using MAF model), the parameters are drawn from a uniform prior in the range: $\theta_i \in \mathcal{U}(0, 1)$. It can be seen that the SBI accurately estimates the parameters G and α using functional data features (FC and FCD) from a budget of 100k simulations (figures 2(A) and (B)). However, we observe diffuse posterior estimation for β mask, by relying only on the functional data features (figure 2(C)). This is quantified by posterior shrinkage, indicating well-identified Bayesian estimation for parameters G and α, but poorly identified β mask when using functional data features (figure 2(D)). This is due to the model’s insensitivity to intra-hemispheric deterioration compared to the whole-network and intra-hemispheric degradation, as calculated from the eigenvalues of the posterior density (figure 2(E)). This is in agreement with the sensitivity analysis conducted using the Hessian matrix, a metric describing the local curvature of a function based on its second partial derivatives (figure 2(F)). The sensitivity of the β mask is nearly zero using Kullback–Leibler divergence or Kolmogorov–Smirnov distance as the fitness function (see supplementary, table S1 and figure S4). Moreover, their values indicate more sensitivity of parameters estimation to FCD than FC (see supplementary, figure S5).

Zoom In Zoom Out Reset image size

To address the practical non-identifiability issue in estimating the β mask, we included the statistical characteristics of the BOLD time-series (such as the mean, variance, skewness, and kurtosis of BOLD time-series) in the training step. As shown in figures 2(A) and (C), the training on spatio-temporal data features results in the better estimation of model parameters, particularly in the estimation of β mask. This is also confirmed by posterior shrinkages, eigenvalues, and Hessian values, as shown in figures 2(D)–(F).

Furthermore, the effect of the simulation budget on the uncertainty quantification when training with the state-of-the-art deep neural density estimators (MAF and NSF) is shown in figure S6. We observed that both MAF and NSF models are compatible in uncertainty quantification (with different simulation budgets), but MAF was 2–4 times faster than NSF during the training process.

When integrating spatio-temporal data in training, some features such as mean and variance may become redundant, especially when using signal processing techniques like normalization or z-scoring. To address this issue, we have perturbed the brain dynamics through stimulation. Figure 3 illustrates the consequence of brain stimulation on the recovery of generative parameters denoted by $\vec{\theta} = \{G, \alpha, \beta\} \in \mathbb{R}^{3}$. Here we used only functional data features for training in SBI.

Zoom In Zoom Out Reset image size

Figures 3(A) and (E) show the simulated BOLD data and the corresponding FC and FCD matrices, in the absence and presence of intervention by stimulation, respectively. It can be seen that the BOLD data display enhanced structural patterns, as evidenced by the higher FC and the increased variance in FCD (the fluidity before stimulation was 0.001, and it increased to 0.004 after stimulation). Here the stimulation is induced by adding a repetitive step current, characterized by a duration of 5 s and an amplitude of $1.2~\mu {\textrm{A cm}^{-2}}$ to the membrane potential variable (through I_stim in equation (1)), within the limbic system.

Figures 3(B) and (F) show that SBI provides accurate recovery of global scaling factor G, before and after stimulation, respectively. This is consistent with our previous results indicating the model’s heightened sensitivity to this parameter. Figures 3(C) and (G) demonstrate the robustness of estimation on the α mask by closely matching the true values, regardless of the stimulation. This indicates the reliability of SBI in estimating the level of degradation between hemispheres, even by perturbation in brain dynamics.

Figures 3(D) and (H) show that SBI provides poor estimation on the β mask when the stimulation is off, but reliable inference when the stimulation is applied, respectively. This demonstrates that the intervention by stimulation is able to address the non-identifiability issue encountered when estimating degradation within brain hemispheres. Perturbing brain dynamics through stimulation increases the likelihood of transitions between brain states, thereby improving the quality of estimation by providing new observable responses for the inference process. Although it may be impractical to administer stimulation to all subjects, VBMs allow us to run simulations in silico and make inferences to evaluate our hypotheses when real experiments are prohibitive.

In the following, we extend the dimensionality of generative parameters to include the excitability parameter. We investigate different forms of model parameterization, such as pooled (homogeneous), unpooled (heterogeneous), and hierarchical (multi-level), to improve the efficiency of SBI-VBMs in higher dimensions, with no need for stimulation (see supplementary, figure S3).

Here we show the estimation over set of unknown generative parameters denoted by $\vec{\theta} = \{G, \alpha, \beta, \eta, J, \Delta \} \in \mathbb{R}^{6}$. In this case, we infer the brain dynamics from parameters located in multiple brain regions with similar or homogeneous characteristics. This form of model parameterization (see supplementary, figure S3(A)) reduces the spatial dimensions, hence the computational complexity of the inference process.

Figure 4 shows the estimated posterior distributions using both spatio-temporal and functional data features extracted from BOLD time-series. The diagonal panels show that the ground-truth values (in green vertical lines) are well under the support of the estimated posterior distributions (in red). Here, we used 100k random simulations from uniform priors (in blue) as: $G, \alpha, \beta \in \mathcal{U}(0, 1)$, $\eta \in \mathcal{U}(-6, -3.5)$, $J \in \mathcal{U}(1, 30)$, and $\Delta \in \mathcal{U}(0.1, 2)$.

Zoom In Zoom Out Reset image size

Due to sufficient budget of simulations and the informativeness of the data features for training, hence the accurate parameter estimation, we can see a close agreement between the data features in the observed and predicted BOLD data (see FC/FCD matrices, shown in the lower diagonal panels). The joint posterior between parameters are shown at the upper diagonal panels, along with their correlation values (ρ shown at the upper left corners). It can be seen that there is a strong correlation between parameters η and J (with ρ = −0.9), also between η and Δ (with ρ = −0.9). This indicates a parameter degeneracy between the excitability and the spread of their distribution, as well the synaptic weight at whole-brain level. This can be due to the homogeneous parameterization used for inference, which ignores all variation among the units being sampled. Using MAF model, the training took around 45 min, whereas generating 10k samples from the posterior took less than 2 s. Supplementary figure S7 demonstrates that excluding the spatio-temporal data features (i.e. statistical characteristics of BOLD time-series) during training leads to poor parameter estimation.

Here we use SBI to estimate posterior distribution in the set of unknown generative parameters denoted by $\vec{\theta} = \{G, \alpha, \beta, \eta_i \} \in \mathbb{R}^{91}$, where η_i with $i \in \{1, 2, \ldots, N_n = 88\}$ is the excitability parameter heterogeneously distributed across brain regions. This form of model parameterization (see supplementary, figure S3(B)) offers a comprehensive explanation for the variation in the excitability across brain regions.

Figure 5 shows the maximum of estimated posterior distribution for three different scenarios in setting the excitabilities across brain regions. It can be seen that SBI provides accurate posterior estimation for various scenarios of excitability maps across brain regions. Here, the training process involved using both spatio-temporal and functional data features. See supplementary, figure S8 for the setups and full estimated posteriors.

Zoom In Zoom Out Reset image size

Importantly, our results indicate that using only the whole-brain functional features (FC and FCD) fails to offer sufficient information for inferring heterogeneous excitability across regions (see supplementary, figure S9).

Regarding the computational cost, using MAF model, the training on 1M simulations took around $10~h$, whereas generating 10k samples from the posterior took less than 60 s. Supplementary figures S10 and S11 demonstrate that both MAF and NSF models were compatible in uncertainty quantification of heterogeneous parameters. Note that here we used 1M random simulations, compared to the 100k simulations used in homogeneous parameterization (figure 4). This indicates the demand for a very large simulation budget during training step to achieve accurate estimation based on a heterogeneous parameterization.

Here we use SBI to estimate posterior distribution in the set of unknown generative parameters denoted by $\vec{\theta} = \{G, \alpha, \beta, \mu_{\eta}, \sigma_{\eta} \} \in \mathbb{R}^{5}$, where $\eta_i \sim \mathcal{N}(\mu_{\eta},\,\sigma_{\eta}^{2})$, indicating that excitabilities are sampled from a Gaussian distribution with mean centered at µ_η , and standard deviation of σ_η . This is a hierarchical parameterization on excitabilities (see supplementary, figure S3(C)), in which they are sampled from a population distribution so that the heterogeneous parameters η_i with $i \in \{1, 2, \ldots, N_n\}$ are now sampled from a generative distribution $\mathcal{N}(\mu_{\eta}, \sigma_{\eta})$. Hence, only two hyperparameters µ_η and σ_η need to be estimated to explain the variation in excitability across the brain, rather than N_n excitability parameters.

Figure 6 shows the estimated posterior by SBI using the both spatio-temporal and functional data features for training. The diagonal panels show that the posterior (in red) accurately encompass the ground-truth values (in green vertical lines) used for generating the observed data. For training, we generated 100k random simulations from uniform priors (in blue) as: $G, \alpha, \beta \in \mathcal{U}(0, 1)$, $\mu_{\eta} \in \mathcal{U}(-6, -4)$, and $\sigma_{\eta} \in \mathcal{U}(0, 1)$. Due to effective form of model parameterization, hence the accurate parameter estimation, we can see a close agreement between the data features in the observed and predicted BOLD data, such as the FC/FCD matrices (shown at the lower diagonal panels). The joint posterior between parameters are shown at the upper diagonal panels, along with their correlation values (ρ at the upper left corners). Using MAF model, the training took around $30~min$, whereas generating 10k samples from the posterior took less than 1 s. In contrast to heterogeneous modeling, the hierarchical structure yields more informative posteriors when using only functional data features, even with a significantly smaller budget of simulations for training (see supplementary, figure S9 versus figure S12). These results indicate that hierarchical approach effectively takes into account the variability across brain regions while naturally incorporating a form of regularization to provide more stable and robust estimates.

Zoom In Zoom Out Reset image size

Lastly, we evaluated the posterior fit in all aforementioned scenarios by plotting z-scores versus shrinkage (see supplementary, figures S13 and S14). We observed that using only functional data features (FC/FCD), we can achieve ideal estimation of the β mask through stimulation and of excitability η through hierarchical parameterization. However, when incorporating all available data features (both spatio-temporal and functional), we achieved an ideal Bayesian inversion for all parameters regardless of parameterization, except for the β mask, which requires stimulation for ideal estimation. These results highlight the importance of brain dynamic perturbation by stimulation and the use of effective Bayesian modeling such as hierarchical when spatio-temporal information is lacking.