1 Introduction

Diffusion-weighted Magnetic Resonance Imaging (dMRI) is a promising technique for characterizing brain microstructure in-vivo using a paradigm called microstructure imaging [Novikov et al., 2018a, Alexander et al., 2019, Jelescu et al., 2020]. Traditionally, microstructure imaging quantifies histologically meaningful features of brain microstructure by fitting a forward (biophysical) model voxel-wise to the set of signals obtained from images acquired with different sensitivities, yielding model parameter maps [Alexander et al., 2019].

Most commonly used techniques rely on a non-linear curve fitting of the signal and return the optimal solution, i.e. the best parameters guess of the fitting procedure. However, this may hide model degeneracy, that is all the other possible estimates that could explain the observed signal equally well [Jelescu et al., 2016]. Another crucial consideration in model fitting is accounting for the uncertainty in parameter estimates. This uncertainty serves various purposes, including assessing result confidence [Jones, 2003], quantifying noise effects [Behrens et al., 2003] or assisting in experimental design [Alexander, 2008].

Instead of attempting to remove the degeneracies, which has been the focus of a large number of studies [Palombo et al., 2023, De Almeida Martins et al., 2021, Slator et al., 2021, Jelescu et al., 2022, Warner et al., 2023, Uhl et al., 2024, Mougel et al., 2024, Olesen et al., 2022, Palombo et al., 2020, Howard et al., 2022, Jones et al., 2018, Vincent et al., 2020, Henriques et al., 2021, Afzali et al., 2021, Lampinen et al., 2023, Zhang et al., 2012-07, Guerreri et al., 2023, Gyori et al., 2022, Novikov et al., 2018b], we propose to highlight them and present all the possible parameter values that could explain an observed signal, providing users with more information to make more confident and explainable use of the inference results.

Posterior distributions are powerful tools to characterize all the possible parameter estimations that could explain an observed measurement, their uncertainty, and existing model degeneracy [Box and Tiao, 2011]. Bayesian inference allows for the estimation of these posterior distributions, traditionally approximating them using numerical methods, such as Markov-Chain-Monte-Carlo (MCMC) [Metropolis et al., 1953]. In quantitative MRI, these methods have been used for example to estimate brain connectivity [Behrens et al., 2003], optimize imaging protocols [Alexander, 2008], or infer crossing fibers by combining multiple spatial resolutions [Sotiropoulos et al., 2013]. However, these classical Bayesian inference methods are computationally expensive and time consuming. They also often require adjustments and tuning specific to each biophysical model [Harms and Roebroeck, 2018].

Harnessing a new deep learning architecture for automatic signal feature selection and efficient sampling of the posterior distributions using simulation-based inference [Cranmer et al., 2020, Lueckmann et al., 2017, Papamakarios et al., 2019], here we propose μGUIDE: a general Bayesian framework to estimate posterior distributions of tissue microstructure parameters from any given biophysical model/signal representation. μGUIDE extends and generalises previous work [Jallais et al., 2022] to any forward model and without acquisition constraints, providing fast estimations of posterior distributions voxel-wise. We demonstrate μGUIDE using numerical simulations on three biophysical models of increasing complexity and degeneracy and compare the obtained estimates with existing methods, including the classical MCMC approach. We then apply the proposed framework to dMRI data acquired from healthy human volunteers and participants with epilepsy. μGUIDE framework is agnostic to the origin of the data and the details of the forward model, so we envision its usage and utility to perform Bayesian inference of model parameters also using data from other MRI modalities (e.g. relaxation MRI) and beyond.

2 Results

2.1 Framework overview

The full architecture of the proposed Bayesian framework, dubbed μGUIDE, is presented in Fig. 1. μGUIDE allows to efficiently estimate full posterior distributions of tissue parameters. It is comprised of two modules that are optimized together to minimize the Kullback-Leibler divergence between the true posterior distribution and the estimated one for every parameters of a given forward model. The ‘Neural Posterior Estimator’ (NPE) module [Papamakarios et al., 2017] uses normalizing flows [Papamakarios et al., 2021] to approximate the posterior distribution, while the ‘Multi-Layer Perceptron’ (MLP) module is used to reduce the data dimensionality and ensure fast and robust convergence of the NPE module.

Figure 1:

The full posterior distribution contains a lot of useful information. To summarize and easily visualize this information, we propose three measures that quantify the best estimates and the associated confidence levels, and a way to highlight degeneracy. The three measures are the Maximum A Posteriori (MAP), which corresponds to the most likely parameter estimate; an uncertainty measure, which quantifies the dispersion of the 50% most probable samples using the interquartile range, relative to the prior range; and an ambiguity measure, which measures the Full Width at Half Maximum (FWHM), in percentage with respect to the prior range. Fig. 2 presents those measures on exemplar posterior distributions.

Figure 2:

We show exemplar applications of μGUIDE to three biophysical models of increasing complexity and degeneracy from the dMRI literature: Ball&Stick [Behrens et al., 2003] (Model 1); Standard Model (SM) [Novikov et al., 2018a] (Model 2); and extended-SANDI [Palombo et al., 2020] (Model 3).

2.2 Evaluation of μGUIDE on simulations

2.2.1 Comparison with MCMC

We performed a comparison between the posterior distributions obtained using μGUIDE and MCMC, a classical Bayesian method. Fig. 3A shows posterior distributions on three exemplar simulations with SNR = 50 using the model 2 (SM), obtained with 15000 samples. Sharper and less biased posterior estimations are obtained using μGUIDE. Fig. 3B presents histograms for each model parameter of the bias between the ground truth value used to simulate a signal, and the MAP of the posterior distributions obtained with either μGUIDE or MCMC, on 200 simulations. Results indicate that the bias is similar or smaller using μGUIDE. Overall, μGUIDE posterior distributions are more accurate than the ones obtained with MCMC.

Figure 3:

Moreover, it took on average 29.3 s to obtain the posterior distribution using MCMC on a GPU (NVIDIA GeForce GT 710) for one voxel, while it only took 0.02 s for μGUIDE. μGUIDE is about 1500 times faster than MCMC, which makes it more suitable for applying it on large datasets.

2.2.2 The importance of feature selection

Fig. 4 shows the MAP extracted from the posterior distributions versus the ground truth parameters used to generate the diffusion signal with μGUIDE and manually-defined summary statistics for the three models. Less biased MAPs with lower ambiguities and uncertainties are obtained with μGUIDE, indicating that the MLP allows for the extraction of additional information not contained in the summary statistics, helping to solve the inverse problem with higher accuracy and precision. μGUIDE generalizes the method developed in [Jallais et al., 2022] to make it applicable to any forward model and any acquisition protocol, while making the estimates more precise and accurate thanks to the automatic feature extraction.

Figure 4:

2.2.3 μGUIDE highlights degeneracies

Fig. 5 presents the posterior distributions of microstructure parameters for the three models obtained with μGUIDE on exemplar noise-free simulations. Blue curves correspond to non-degenerate posterior distributions, while the red ones present at least one degeneracy for one of the parameters. As the complexity of the model increases, degeneracy in the model definitions appear. This figure showcases μGUIDE ability to highlight degeneracy in the model parameter estimation.

Figure 5:

Tables 1 and 2 present the number of degenerate cases for each parameter in the three models, on 10000 simulations. Table 1 considers noise-free simulations and the training and estimations were performed on CPU. Table 2 reports results on noisy simulations (Rician noise with SNR = 50), with training and testing performed on a GPU (NVIDIA GeForce RTX 4090). The time needed for the inference and to estimate the posterior distributions on 10000 simulations, define if they are degenerate or not, and extract the MAP, uncertainty, ambiguity, are also reported. The more complex the model, the more degeneracies.

Table 1:

Table 2:

2.3 Application of μGUIDE to real data

After demonstrating that the proposed framework provides good estimates in the controlled case of simulations, we applied μGUIDE to both a healthy volunteer and a participant with epilepsy. The estimation of the posterior distributions is done independently for each voxel. To easily assess the values and the quality of the fitting, we are plotting the maximum-a-posterior, ambiguity and uncertainty maps, but the full posterior distributions are stored and available for all the voxels. Voxels presenting a degeneracy are highlighted with a red dot.

2.3.1 Healthy volunteer

We applied μGUIDE to a healthy volunteer, using the Ball&Stick, SM and extended-SANDI models. Fig. 6 presents the parametric maps of an exemplar set of model parameters for each model, alongside their degeneracy, uncertainty and ambiguity. The Ball&Stick model presents no degeneracy, the SM presents some degeneracy, mostly in voxels with high likelihood of partial voluming with cerebrospinal fluid (CSF) and at the white matter - gray matter boundaries. The extended-SANDI model is the model showing the highest number of degenerate cases, mostly localized within the white matter areas characterized by complex microstructure, e.g. crossing fibers. This result is expected, as the complexity of the models increases, leading to more combinations of tissue parameters that can explain an observed signal. Measures of ambiguity and uncertainty allow to quantify the confidence in the estimates and help interpreting the results.

Figure 6:

2.3.2 Participant with epilepsy

Fig. 7 demonstrates μGUIDE application to a participant with epilepsy, using the SM. Noteworthy, the axonal signal fraction estimates within the epileptic lesion show low uncertainty and ambiguity measures hence high confidence, while orientation dispersion index estimates show high uncertainty and ambiguity suggesting low confidence, cautioning the interpretation.

Figure 7:

3 Discussion

3.1 Applicability of μGUIDE to multiple models

The μGUIDE framework offers the advantage of being easily applicable to various biophysical models and representations, thanks to its data-driven approach for data reduction. The need to manually define specific summary statistics that capture the relevant information for microstructure estimation from the multi-shell diffusion signal is removed. This also eliminates the acquisition constraints that were previously imposed by the summary statistics definition [Jallais et al., 2022]. The extracted features contain additional information compared to the summary statistics (see Appendix A), resulting in a notable reduction in bias (on average 5.2 fold lower), uncertainty (on average 2.6 fold lower) and ambiguity (on average 2.7 fold lower) in the estimated posterior distributions. Consequently, μGUIDE improves parameters estimation over current state-of-the-art methods (e.g. Jallais et al. [2022]), showing for example reduced bias (on average 5.2 fold lower) and dispersion (on average 6.4 fold lower) on the MAP estimates for each of the three example models investigated (see Fig. 4).

In this study, we presented applications of μGUIDE to brain microstructure estimation using three well-established biophysical models, with increased complexity: the Ball&Stick model, the Standard Model, and an extended-SANDI model. However, our approach is not limited to brain tissue nor to diffusion-weighted MRI and can be extended to different organs by employing their respective acquisition encoding and forward models, such as NEXI for exchange estimates [Jallais et al., 2024]; mcDESPOT for myelin water fraction mapping using quantitative MRI relaxation [Deoni et al., 2008]; VERDICT in prostate imaging [Panagiotaki et al., 2014]; or even adapted to different imaging modalities (e.g., EEG and MEG), where there is a way to link (via modelling or simulation) the observed signal to a set of parameters of interest. This versatility underscores the broad applicability of our proposed approach across various biological systems and imaging techniques.

It is important to note that μGUIDE is still a model-dependent method, meaning that the training process is based on the specific model being used. Additionally, the number of features extracted by the MLP needs to be predetermined. One way to determine the number of features is by matching it with the number of parameters being estimated. Alternatively, a dimensionality-reduction study using techniques like t-distributed stochastic neighbour embedding (tSNE) [Van der Maaten and Hinton, 2008] can be conducted to determine the optimal number of features.

3.2 μGUIDE: an efficient framework for Bayesian inference

One notable advantage of μGUIDE is its amortized nature. With this approach, the training process is performed only once, and thereafter, the posterior estimations can be independently obtained for all voxels. This amortization enables efficient estimations of the posterior distributions. μGUIDE outperforms in terms of speed conventional Bayesian inference methods such as MCMC, showing a “--1500 fold acceleration. The time savings achieved with μGUIDE make it a highly efficient and practical tool for estimating posterior distributions in a timely manner.

This unlocks the possibility to process with Bayesian inference very large datasets in manageable time (e.g. approximately 6 months to process 10k dMRI datasets) and to include Bayesian inference in iterative processes that require the repeated computation of the posterior distributions (e.g., dMRI acquisition optimization Alexander [2008]).

In the dMRI community, the use of SBI methods to characterize full posterior distributions as well as quantify the uncertainty in parameter estimations was first introduced in Jallais et al. [2022] for a grey matter model. An application to crossing fibers has recently been proposed by Karimi et al. [2024]. Those approaches use different density estimators. This work and Jallais et al. [2022] rely on Masked Autoregressive FLows (MAF) Papamakarios et al. [2017], while the work by Karimi et al. [2024] is based on Mixture Density Networks (MDN) [Bishop, 1994]. MAFs have been found to show superior performance compared to MDNs [Goncalves et al., 2020, Patron et al., 2022].

3.3 μGUIDE quantifies confidence to guide interpretation

Quantifying confidence in an estimate is of crucial importance. As demonstrated by our pathological example, changes in the tissue microstructure parameters can help clinicians decide which parameters are the most reliable and better interpret microstructure changes within diseased tissue. On large population studies, the quantified uncertainty can be taken into account when performing group statistics and to detect outliers.

Multiple approaches have been used to try and quantify this uncertainty. Gradient descent often provides a measure of confidence for each parameter estimate. Alternative approaches use the shape of the fitted tensor itself as a measure of uncertainty for the fiber direction [Koch et al., 2002, Parker and Alexander, 2003]. Other methods also rely on bootstrapping techniques to estimate uncertainty. Repetition bootstrapping for example depends on repeated measurements of signal for each gradient direction, but imply a long acquisition time and cost, and are prone to motion artifacts [Lazar and Alexander, 2005, Jones, 2003]. In contrast, residual bootstrapping methods resample the residuals of a regression model. Yet, this approach is heavily dependent on the model and can lead to overfitting [Whitcher et al., 2008, Chung et al., 2006]. In general, resampling methods can be problematic for sparse samples, as the bootstrapped samples tend to underestimate the true randomness of the distribution [Kauermann et al., 2009]. We propose to quantify the confidence by estimating full posterior distributions, which also has the benefit of highlighting degeneracy.

Variance observed in the posterior distributions can be attributed to several factors. The presence of noise in the signal contributes to irreducible variance, decreasing the confidence in the estimates as the noise level increases. Another source of variance can arise from the choice of acquisition parameters. Different acquisitions may provide varying levels of confidence in the parameter estimates. Under-sampled acquisitions or inadequate b-shells may fail to capture essential information about a tissue microstructure, such as soma or neurite radii, resulting in inaccurate estimates.

μGUIDE can guide users in determining whether an acquisition is suitable for estimating parameters of a given model and viceversa, the variance and bias of the posterior distributions estimated with μGUIDE can be used to guide the optimization of the data acquisition to maximize accuracy and precision of the model parameters estimates.

The presence of degeneracy in the solution of the inverse problem is influenced by the complexity of the model being used and the lack of sufficient information in the data. In recent years, researchers have introduced increasingly sophisticated models to better represent the brain tissue, such as SANDI [Palombo et al., 2020], NEXI [Jelescu et al., 2022] and eSANDIX [Olesen et al., 2022], that take into account an increasing number of tissue features. By applying μGUIDE, it becomes possible to gain insights into the degree of degeneracy within a model and to assess the balance between model realism and the ability to accurately invert the problem. We have recently provided an example of such application for NEXI and SANDIX [Jallais et al., 2024].

3.4 Summary

We propose a general Bayesian framework, dubbed μGUIDE, to efficiently estimate posterior distributions of tissue microstructure parameters. For any given acquisition and signal model/representation, μGUIDE improves parameters estimation and computational time over existing state-of-the-art methods. It allows to highlight degeneracy, and quantify confidence in the estimates, guiding results interpretation towards more confident and explainable diagnosis using modern deep learning. μGUIDE is not inherently limited to dMRI and microstructure imaging. We envision its usage and utility to perform efficient Bayesian inference also using data from any modality where there is a way to link (via modelling or simulation) the observed measurements to a set of parameters of interest.

4 Methods

4.1 Solving the inverse problem using Bayesian inference

4.1.1 The inference problem

We make the hypothesis that an observed dMRI signal x₀ can be explained (and generated) using a handful of relevant tissue microstructure parameters θ₀, following the definition of a forward model:

The objective is, given this observation x₀, to estimate the parameters ∞₀ that generated it.

Forward models are designed to mimic at best a given biophysical phenomenon, for some given time and scale [Alexander, 2009, Yablonskiy and Sukstanskii, 2010, Jelescu and Budde, 2017, Novikov et al., 2018a, Alexander et al., 2019, Jelescu et al., 2020]. As a consequence, forward models are injection functions (every biologically plausible ∞_i generates exactly one signal x_i), but do not always happen to be bijections, meaning that multiple ∞_i can generate the same signal x_i. It can be impossible, based on biological considerations, to infer which solution ∞_i best reflects the probed structure. We refer to these models as ‘degenerate models’.

Point estimates algorithms, such as minimum least square or maximum likelihood estimation algorithms, allow to estimate one set of microstructure parameters that could explain an observed signal. In the case of degenerate models, the solution space can be multi-modal and those algorithms will hide possible solutions. When considering real-life acquisitions, i.e. noisy and/or under-sampled acquisitions, one also needs to consider the bias introduced with respect to the forward model, and the resulting variance in the estimates [Jones, 2003, Behrens et al., 2003].

We propose a new framework that allows for the estimation of full posterior distributions p(∞|x₀), that is all the probable parameters that could represent the underlying tissue, along with an uncertainty measure and the interdependency of parameters. These posteriors can help interpreting the obtained results and make more informed decisions.

4.1.2 The Bayesian formalism

The posterior distribution can be defined using Bayes’ theorem as follows:

where p(x₀|θ) is the likelihood of the observed data point, p(θ) is the prior distribution defining our initial knowledge of the parameter values, and p(x₀) is a normalizing constant, commonly referred to as the evidence of the data.

The evidence term is usually very hard to estimate, as it corresponds to all the possible realisations of x₀, i.e. . For simplification, methods usually estimate an unnormalized probability density function, i.e.

To approximate these posterior distributions, traditional methods rely on the estimation of the likelihood p(x₀|θ) of the observed data point x₀ via an analytic expression. This likelihood function corresponds to an integral over all possible trajectories through the latent space, that is p(x₀|θ) = ∫ p(x₀, z|∞)dz, where p(x₀, z|∞) is the joint probability density of observed data x₀ and latent variables z. For forward models with large latent spaces, computing this integral explicitly becomes impractical. The likelihood function is then intractable, rendering these methods unusable [Cranmer et al., 2020]. Models that do not admit a tractable likelihood are called implicit models [Diggle and Gratton, 1984].

To circumvent this issue, some techniques have been proposed to sample numerically from the likelihood function, such as Markov-Chain-Monte-Carlo (MCMC) [Metropolis et al., 1953]. Another set of approaches proposes to train a conditional density estimator to learn a surrogate of the likelihood distribution [Papamakarios et al., 2019, Lueckmann et al., 2019], the likelihoodratio [Cranmer et al., 2016, Gutmann et al., 2018] or the posterior distribution [Papamakarios and Murray, 2016, Lueckmann et al., 2017, Papamakarios et al., 2019], allowing to greatly reduce computation times. These methods are dubbed Likelihood-Free Inference (LFI), or SimulationBased Inference (SBI) methods [Cranmer et al., 2020, Tejero-Cantero et al., 2020]. In particular, there has been a growing interest towards deep generative modeling approaches in the machine learning community [Lueckmann et al., 2021]. They rely on specially tailored neural network architectures to approximate probability density functions from a set of examples. Normalizing flows [Papamakarios et al., 2021] are a particular class of such neural networks that have demonstrated promising results for SBI in different research fields [Goncalves et al., 2020, Greenberg et al., 2019].

While this work focuses on the estimate of the posterior distribution using a conditional density estimator, we show a comparison with MCMC, which are commonly used methods in the community. We will therefore introduce this method in the following paragraph.

4.1.3 Estimating the likelihood function

Well-established approaches for estimating the likelihood function are MCMC methods. These methods rely on a noise model to define the likelihood distribution, such as the Rician [Panagiotaki et al., 2012] or Offset Gaussian models [Alexander, 2009]. In this work, we will be using the Microstructure Diffusion Toolbox to perform the MCMC computations [Harms and Roebroeck, 2018], which relies on the Offset Gaussian model. The log-likelihood function is then the following:

where M(∞) is the signal obtained using the biophysical model, σ is the standard deviation of the Gaussian distributed noise, estimated from the reconstructed magnitude images [Dietrich et al., 2007], and m is the number of observations in the dataset.

MCMC methods allow to obtain posterior distributions using Bayes’ formula (Eq. (2)) with the previously defined likelihood function (Eq. (3)) and some prior distributions, which are usually uniform distributions defined on biologically plausible ranges. They generate a multi-dimensional chain of samples which is guaranteed to converge towards a stationary distribution, which approximates the posterior distribution [Metropolis et al., 1953].

The need to compute the signal following the forward model at each iteration makes these sampling methods computationally expensive and time consuming. In addition, they require some adjustments specific to each model, such as the choice of burn-in length, thinning and the number of samples to store. Harms and Roebroeck [2018] recommend to use the Adaptive Metropolis-Within-Gibbs (AMWG) algorithm for sampling dMRI models, initialized with a maximum likelihood estimator obtained from non-linear optimization, with 100 to 200 samples for burn-in and no thinning. The recommended number of samples is model-dependent. Authors recommendations can be found in their paper.

4.1.4 Bypassing the likelihood function

An alternative method was proposed to overcome the challenges associated with approximating the likelihood function and the limitations of MCMC sampling algorithms. This approach involves directly approximating the posterior distribution by using a conditional density estimator, i.e. a family of conditional probability density function approximators denoted as q_ϕ(∞|x). These approximators are parameterized by ϕ and accept both the parameters θ and the observation x as input arguments. Our posterior approximation is then obtained by minimizing its average Kullback-Leibler divergence with respect to the conditional density estimator for different choices of x, as per [Papamakarios and Murray, 2016]:

which can be rewritten as

where C is a constant that does not depend on ∞. Note that in practice we consider a N-sample Monte-Carlo approximation of the loss function:

where the N data points (θ_i, x_i) are sampled from the joint distribution with θ_i ~ p(θ) and x_i ~ p(x|θ_i). We can then use stochastic gradient descent to obtain a set of parameters ϕ which minimizes L^N.

If the class of conditional density estimators is sufficiently expressive, it can be demonstrated that the minimizer of Eq. (6) converges to p(θ|y) when N → ∞ [Greenberg et al., 2019]. It is worth noting that the parametrization ϕ, obtained at the end of the optimization procedure, serves as an amortized posterior for various choices of x. Hence, for a particular observation x₀, we can simply use q_ϕ(θ|x₀) as an approximation of p(θ|x₀).

4.2 μGUIDE framework

The full architecture of the proposed Bayesian framework, dubbed μGUIDE, is presented in Fig. 1. The analysis codes underpinning the results presented here can be found upon publication in the Cardiff University data catalogue and on Github: https://github.com/mjallais/uGUIDE (both CPU and GPU are supported).

μGUIDE is comprised of two modules that are optimized together to minimize the Kull-back-Leibler divergence between the true posterior distribution and the estimated one for every parameters of a given forward model. The ‘Neural Posterior Estimator’ (NPE) module uses normalizing flows to approximate the posterior distribution, while the ‘Multi-Layer Perceptron’ (MLP) module is used to reduce the data dimensionality and ensure fast and robust convergence of the NPE module. The following Sections provide more details about our implementation of each module.

4.2.1 Neural Posterior Estimator

In this study, the Sequential Neural Posterior Estimation (SNPE-C) algorithm [Papamakarios and Murray, 2016, Greenberg et al., 2019] with a single round is employed to train a neural network that directly approximates the posterior distribution. Thus, sampling from the posterior can be done by sampling from the trained neural network. Neural density estimators have the advantage of providing exact density evaluations, in contrast to Variational Autoencoders (VAE) Kingma and Welling [2019] or generative adversarial networks (GAN) Goodfellow et al. [2014], which are better suited for generating synthetic data.

The conditional probability density function approximators used in this pro ject belong to a class of neural networks called normalizing flows [Papamakarios et al., 2021]. These flows are invertible functions capable of transforming vectors generated from a simple base distribution (e.g. the standard multivariate Gaussian distribution) into an approximation of the true posterior distribution. An autoregressive architecture for normalizing flows is employed, implemented via the Masked Autoregressive Flow (MAF) Papamakarios et al. [2017], which is constructed by stacking five Masked Autoencoder for Distribution Estimation models (MADE) Germain et al. [2015]. An explanation of how MAF and MADE work is provided in Appendix B.

To test that the predicted posteriors for a given model are not incorrect we use posterior predictive checks, which is described in more details in Appendix C.

4.2.2 Handling the large dimensionality of the data with Multi-Layer Perceptron

As the dimensionality of the input data x grows, the complexity of the corresponding inverse problem also increases. Accurately characterizing the posterior distributions or estimating the tissue microstructure parameters becomes more challenging. As a consequence, it is often necessary to rely on a set of low-dimensional features (or summary statistics) instead of the raw data for the inference task process [Blum et al., 2013, Fearnhead and Prangle, 2012, Papamakarios et al., 2019]. These summary statistics are features that capture the essential information within the raw data, allowing to reduce the size of the input vector. Learning a set of sufficient statistics before estimating the posterior distribution makes the inference easier and offers many benefits (see e.g. the Rao-Blackwell theorem).

A follow-up challenge lies in the choice of suitable summary statistics. For well-understood problems and data, it is possible to manually design these features using deterministic functions that condense the information contained in the raw signal into a set of handful summary statistics. Previous works, such as Novikov et al. [2018b] and Jallais et al. [2022], have proposed specific summary statistics for two different biophysical models. However, defining these summary statistics is difficult and often requires prior knowledge of the problem at hand. In the context of dMRI, they also rely on acquisition constraints and are model-specific.

In this work, the proposed framework aims to be applicable to any forward model and be as general as possible. We therefore propose to learn the summary statistics from the high-dimensional input signals x using a neural network. This neural network is referred to as an embedding neural network. The observed signals are fed into the embedding neural network, whose outputs are then passed to the neural density estimator. The parameters of the embedding network are learned together with the parameters of the neural density estimator, leading to the extraction of optimal features that minimize the uncertainty of p(ϕ|x). Here, we propose to use a MLP with three layers as a summary statistics extractor. The number of features N_f extracted by the MLP can be either defined a priori or determined empirically during training.

4.2.3 Training μGUIDE

To train μGUIDE we need couples of input vectors x and corresponding ground-truth values for the model parameters that we want to estimate, θ. The input x can be real or simulated data (e.g. dMRI signals); or a mixture of these two. We train μGUIDE by stochastically minimizing the loss function defined in Eq. (6) using the Adam optimizer [Kingma and Ba, 2015] with a learning rate of 10⁻³ and a minibatch size of 128. We use 1 million simulations for each model, 5% of which are randomly selected to be used as a validation set. Training is stopped when the validation loss does not decrease for 30 consecutive epochs.

4.2.4 Quantifying the confidence in the estimates

The full posterior distribution contains a lot of useful information about a given model parameter best estimates; uncertainty; ambiguity and degeneracy. To summarize and easily visualize this information, we propose three measures that quantify the best estimates and the associated confidence levels, and a way to highlight degeneracy.

We start by checking whether a posterior distribution is degenerate, that is if the distribution presents multiple distinct parameter solutions, appearing as multiple local maxima (Fig. 2). To that aim, we fit two Gaussian distributions to the obtained posterior distributions. A voxel is considered as degenerate if the derivative of the fitted Gaussian distributions changes signs more than once (i.e. multiple local maxima), and if the two Gaussian distributions are not overlapping (the distance between the two Gaussian means is inferior to the sum of their standard deviations).

For non-degenerate posterior distributions, we extract three quantities:

1. The Maximum A Posteriori (MAP), which corresponds to the most likely parameter estimate.

2. An uncertainty measure, which quantifies the dispersion of the 50% most probable samples using the interquartile range, relative to the prior range.

3. An ambiguity measure, which measures the Full Width at Half Maximum (FWHM), in percentage with respect to the prior range.

Fig. 2 presents those measures on exemplar posterior distributions.

4.3 Application of μGUIDE to biophysical modelling of dMRI data

We show exemplar applications of μGUIDE to three biophysical models of increasing complexity and degeneracy from the dMRI literature. For each model, we compare the fitting quality of the posterior distributions obtained using the MLP and manually defined summary statistics.

4.3.1 Biophysical models of dMRI signal.

Model 1: Ball&Stick [Behrens et al., 2003]. This is a two-compartment model (intra-neurite and extra-neurite space) where the dMRI signal from the brain tissue is modeled as a weighted sum, with weight f_in, of signals from water diffusing inside the neurites, approximated as sticks (i.e. cylinders of zero radius) with diffusivity D_in, and water diffusing within the extra-neurite space, approximated as Gaussian diffusion in an isotropic medium with diffusivity D_e. The direction of the stick is randomly sampled on a sphere. This model has the main advantage of being nondegenerate. We define the summary statistics as the direction-averaged signal (six b-shells, see Section 4.3.3).

Model 2: Standard Model (SM) [Novikov et al., 2018a]. Expanding on model 1, this model represents the dMRI signal from the brain tissue as a weighted sum of the signal from water diffusing within the neurite space, approximated as sticks with symmetric orientation dispersion following a Watson distribution and water diffusing within the extra-neurite space, modelled as anisotropic Gaussian diffusion. The microstructure parameters of this two-compartment model are the neurite signal fraction f , the intra-neurite diffusivity D_a, the orientation dispersion index ODI, and the parallel and perpendicular diffusivities within the extra-neurite space and We use the LEMONADE [Novikov et al., 2018b] system of equations, which is based on a cumulant decomposition of the signal, to define six summary statistics.

with r_s and D_s the soma radius and diffusivity respectively, and α_m the mth root of , with J_n(x) the Bessel functions of the first kind. We use the six summary statistics defined in [Jallais et al., 2022], which are based on a high and low b-value signal expansion.

Prior distributions p(θ) are defined as uniform distributions over biophysically plausible ranges. Signal fractions are defined within the interval [0,1], diffusivities between 0.1 and 3 μm² ms^-1, ODI between 0.03 and 0.95, and C_s between 0.15 and 1105 μm² (which correspond to r_s ∈ [1; 15] μm and fixed D_s = 3 μm² ms^-1).

The SM imposes the constraint < To generate samples uniformly distributed on the space defined by this condition, we are using two random variables u₀ and u₁, both sampled uniformly between 0 and 1, and then relate them to and using the following equations:

The extended-SANDI model requires for the signal fractions to sum to 1, that is f_n+f_s+f_k = 1. To uniformly cover the simplex f_n + f_s + f_e = 1, we define two new parameters k₁ and k₂ , uniformly sampled between 0 and 1, and use the following equations to get the corresponding signal fractions:

To ensure comparability of results, we extract the same number of features N_f using the MLP as the number of summary statistics for each model. We therefore use N_f = 6 for the Ball&Stick, the SM and the extended-SANDI models.

4.3.2 Validation in Numerical Simulations

We start by validating the proposed method using posterior predictive checks (PPC) and simulated signals from model 2 (see more details in Appendix C). Since PPC alone does not guarantee the correctness of the estimated posteriors, we further validated the obtained posterior distributions comparing them with the Adaptive Metropolis-Within-Gibbs MCMC [Roberts and Rosenthal, 2009]. We generated simulations following the same acquisition protocol as the real data (see Section 4.3.3), and added Rician distributed noise with a Signal-to-Noise-Ratio (SNR) of 50. We then estimated the posterior distributions using both μGUIDE and the MCMC method implemented in the MDT toolbox [Harms and Roebroeck, 2018]. We initialized the sampling using a maximum likelihood estimator. We sampled 15200 samples from the distribution, the first 200 ones being used as burn-in, and no thinning. Similarly, we sampled 15000 samples from the estimated posterior distributions using μGUIDE.

Then, we show that μGUIDE can be applicable to any model. We use model 1 and 3 as examples of simpler (and non-degenerate) and more complex (and degenerate) models than model 2, respectively.

We compared the proposed framework to a state-of-the-art method for posterior estimation [Jallais et al., 2022]. This method relies on manually-defined summary statistics, while μGUIDE automatically extracts some features using an embedded neural network. μGUIDE was trained directly on noisy simulations. The manually-defined summary statistics were extracted from these simulated noisy signals and then used as training dataset for a MAF, similarly to [Jallais et al., 2022].

Finally, we used μGUIDE to highlight degeneracy in all the models. While the complexity of the models increases, more degeneracy can be found. The degeneracy is inherent to the model definition, and is not induced by the noise. μGUIDE allows to highlight those degeneracies and quantify the confidence in the obtained estimates.

The training was performed on N = 10⁶ numerical simulations for each model, computed using the MISST package [lanus et al., 2017] and random combinations of the model parameters, each uniformly sampled from the previously defined ranges, with the addition of Rician distributed noise with SNR equivalent to the experimental data, i.e. 50.

4.3.3 dMRI Data Acquisition and Processing

We applied μGUIDE to dMRI data collected from two participants: a healthy volunteer from the WAND dataset [McNabb et al., 2024] and an age-matched participant with epilepsy, acquired with the same protocol used for the MICRA dataset [Koller et al., 2021]. Data were acquired on a Connectome 3T scanner using a single-shot spin-echo, echo-planar imaging sequence with b-values = [200, 500, 1200, 2400, 4000, 6000] s mm^-2, [20, 20, 30, 61, 61, 61] uniformly distributed directions respectively and 13 non-diffusion-weighted images at 2 mm isotropic resolution. TR was set to 3000 ms, TE to 59 ms, and the diffusion gradient duration and separation to 7 ms and 24 ms respectively. Short diffusion times and TE were achieved thanks to the Connectom gradients, allowing to enhance the signal-to-noise ratio and sensitivity to small water displacements [Jones et al., 2018, Setsompop et al., 2013]. We considered the noise as Rician with an SNR of 50 for both subjects.

Data were preprocessed using a combination of in-house pipelines and tools from the FSL [Andersson et al., 2003, Andersson and Sotiropoulos, 2016, Smith, 2002, Smith et al., 2004] and MRTrix3 [Tournier et al., 2019] software packages. The preprocessing steps included brain extraction [Smith, 2002], denoising [Cordero-Grande et al., 2019, Veraart et al., 2016], drift correction [Vos et al., 2017, Sairanen et al., 2018], susceptibility-induced distortions [Andersson et al., 2003, Smith et al., 2004], motion and eddy current correction [Andersson and Sotiropoulos, 2016], correction for gradient non-linearity distortions [Glasser et al., 2013], and Gibbs ringing artefacts correction [Kellner et al., 2016].

4.3.4. dMRI Data Analysis

Diffusion signals were first normalized by the mean non-diffusion-weighted signals acquired for each voxel. Each voxel was then estimated in parallel using the μGUIDE framework. For each observed signal x_v (i.e. for each voxel), we drew 50000 samples via rejection-sampling from q_ϕ(θ_i|x_v) for each model parameter θ_i, allowing to retrieve the full posterior distributions. If a posterior distribution was not deemed degenerate, the maximum-a-posteriori, uncertainty and ambiguity measures were extracted from the posterior distributions.

The manually-defined summary statistics of the SM are defined using a cumulant expansion, which is only valid for small b-values. We therefore only used the b ≤ 2500 smm^-2 data for this model. In order to obtain comparable results, we restricted the application of μGUIDE to this range of b-values as well. An extra b-shell (b-value = 5000 s mm^-2; 61 directions) was interpolated using mapl [Fick et al., 2016] for the extended-SANDI model when using the method developed by Jallais et al. [2022] based on summary statistics.

The training of μGUIDE was performed as described in 4.3.2 and an example of training dataset and input signal vector is provided in Fig. 8.

Figure 8:

All the computations were performed both on CPU and GPU (NVIDIA GeForce RTX 4090).

Acknowledgements

This work, Maeliss Jallais and Marco Palombo are supported by UKRI Future Leaders Fellowship (MR/T020296/2). We are thankful to Dr. Dmitri Sastin and Dr. Khalid Hamandi for sharing their dataset from a participant with epilepsy, and to Dr. Carolyn McNabb, Dr. Eirini Messar-itaki and Dr. Pedro Luque Laguna for preprocessing the data of the healthy participant from the WAND data. The WAND data were acquired at the UK National Facility for In Vivo MR Imaging of Human Tissue Microstructure funded by the EPSRC (grant EP/M029778/1) and The Wolfson Foundation, and supported by a Wellcome Trust Investigator Award (096646/Z/11/Z) and a Wellcome Trust Strategic Award (104943/Z/14/Z).

Appendix

A Correlation analysis between features extracted by μGUIDE and manually-defined summary statistics

Fig. 9 presents the correlation matrices obtained from the correlation between the MLP-extracted features from μGUIDE and the manually-defined summary statistics defined in Section 4.3, considering noisy simulations (SNR=50). For each model, at least one feature extracted by μGUIDE is not or weakly correlated with the summary statistics. Additional information, not contained in the summary statistics, is extracted by the MLP from the input signal, leading to reduced bias, uncertainty and ambiguity in the parameter estimates (see Fig. 4).

Figure 9:

B Masked Autoregressive Flows

The conditional probability density function approximators used in this project belong to a class of neural networks called Normalizing Flows (NF) Papamakarios et al. [2021]. NFs provide a general way of transforming complex probability distributions over continuous random variables into simple base distributions p(z) (such as normal distributions) through a chain of invertible and differentiable transformations f. By applying the change of variable formula, the target distribution q_ϕ (θ|x) can be written as:

where z = f_ϕ(θ; x) is invertible and differentiable (i.e. a diffeomorphism), and J_fϕ (θ; x) is the Jacobian of f_ϕ(θ; x). The forward direction allows for density evaluation, that is learning the mapping between the target and the base distributions, i.e. learning the parameters ϕ. The inverse direction allows to estimate a density estimator q_ϕ (θ | x₀) by sampling points z from the base distribution and applying the inverse transform .

A main requirement is that the flow needs to be expressive enough to approximate any arbitrarily highly complex distribution. An interesting property of diffeomorphisms is that they are closed under composition, which means that a composition of K diffeomorphisms f = f₁ ◦ f₂ ◦ … ◦ f_K is also a diffeomorphism, and the Jacobian determinant is the product of the determinant of each component. Combining multiple transformations allows to increase the expressivity of the general flow. We obtain:

Flows need to be flexible and expressive enough to model any desired distribution but also need to be computationally efficient, that is, computing the associated Jacobian determinants need to be tractable and efficient. Among a number of proposed architectures such as mixture density networks [Bishop, 1994] or neural spline flows [Durkan et al., 2019], we focused on Masked Autoregressive Flows (MAF) [Papamakarios et al., 2017], which has shown state-of-the-art performance as well as the ability to estimate multi-modal posterior distributions [Goncalves et al., 2020, Patron et al., 2022, Papamakarios et al., 2021].

Autoregressive flows are universal approximators and have the form where h_i = c_i(z<i) [Papamakarios et al., 2021]. τ is termed the transformer and is a stritly monotonic function parametrized by h_i, and c_i the i-th conditioner. Each h_i and therefore each can be computed independently in parallel, helping to keep a low computation time. The conditioner constraints each output to depend only on variables with dimension indices less than i, which makes the Jacobian of the flow lower diagonal. Its determinant can then be obtained easily as the product of its diagonal elements.

To efficiently implement the conditioner, this method relies on the Masked Autoencoder for Distribution Estimation (MADE) [Germain et al., 2015] architecture. To create a neural network that obey the autoregressive structure of the conditioner, a fully connected feedforward neural network is multiplied to binary masks, which removes some connections by assigning them a weigh of 0. The binary masks can easily be obtained by following a few simple steps (see Fig. 10 for an illustration):

1. Label the input and output nodes between 1 and D, D being the dimension of the input vector z.

2. Randomly assign each hidden unit a number between 1 and D — 1, which indicates the number of inputs it will be connected to.

3. For each hidden layer, connect the hidden units to units with inferior or equal labels.

4. Connect output units to units with strictly inferior labels.

Figure 10:

For μGUIDE’s implementation, we use a combination of 5 MADEs.

As a result, the MAF architecture only needs a single forward pass through the flow and, combined with the low-cost computation of the determinant, allows for fast training and evaluation of the posterior distributions.

C Posterior Predictive Checks

Posterior Predictive Checks (PPC) are a common safety check to verify inference is not wrong. The idea is to compare input signals with generated signals from samples drawn from the posterior distributions. If the inference is correct, the generated signals should look similar to the input signal.

We performed the following steps:

• Sample N θ_i from the prior distribution: θ_i ~ p(θ)

• Generate the corresponding signals using the forward model: x_i = M(θ_i)

• Perform the inference and estimate the posterior distributions p(θ_i |x_i)

• Sample N_PP samples θ_i,s from p(θ_i|x_i)

• Reconstruct the signals from the sampled θ_i,s using the forward model: x_i,s = M(θ_i,s)

• Compare the obtained x_i,s with x_i.

Fig. 11 presents results on model 2 (SM), on both noise-free and noisy signals (Rician noise with SNR=50) for N =10 random combinations of model parameters, and N_pp = 100. As dMRI data have a high dimensionality, we report the direction-average signal. Plain lines show the signals x_i, and the shaded areas correspond to the area in which the corresponding x_i,s fall. x_i lie within the support of x_i,s, indicating the inference is not wrong. Note that the support of x_i,s is bigger for noisy simulations, reflecting the wider posterior distributions obtained from the inference.

Figure 11: