This method, which relies on maximum likelihood estimation (MLE), extracts biomarkers by performing a voxelwise model fit every time new data is acquired. An appropriate signal model $M$ is required, parameterised by $n_{y}$ parameters of interest; for each combination of $y$, the probability of observing the acquired data $x$ is known as the likelihood $L$ of those parameters:

Developments in qMRI acquisition and analysis have led to increased (i) image spatial resolution (i.e. greater $n_{v}$) and (ii) model complexity (i.e. greater $n_{y}$), such that conventional MLE fitting has become increasingly computationally expensive.

Deep learning approaches address this by reframing $n_{v}$ independent problems into a single global model fit: learning the function $\mathcal{F}$ that maps any $x$ to its corresponding $y_{gt}$:

Deep neural networks aim to approximate this function by composing a large but finite number of building-block functions, parametrised by $n_{p}$ network parameters $p$ (“weights”):

In this context, model fitting (“training”), is performed over network parameters $p$ and involves maximising $\hat{\mathcal{F}}$’s mean performance over a large set of training examples; the trained network is defined by the best-fit parameters $\hat{p}$. This fitting problem, whilst more computationally expensive to solve than any individual voxel ($n_{v}=1$) MLE, is only tackled once; once $\hat{\mathcal{F}}$ is learnt, it can be applied at negligible cost to new, unseen, data. This promise of rapid, zero-cost parameter estimation has led to the development of two broad classes of DL-based parameter estimation methods.

Supervised_GT methods approximate $\mathcal{F}$ by minimising the difference between a large number of noise-free training labels (groundtruth parameter values) and corresponding network outputs (noise-free parameter estimates); training loss is calculated in the parameter space $Y$:

These methods produce higher bias, lower variance parameter estimation than conventional MLE fitting (Grussu et al., 2021; Gyori et al., 2022) and, by adjusting $W$, can be tailored to selectively boost estimation performance on a subset of the parameter space $Y$.

In contrast, Self-supervised methods compute training loss within the signal space $X$, by minimising the difference between network inputs (noisy signals) and a filtered representation of network outputs (noise-free signal estimates):

Under Gaussian noise conditions, single-voxel Self-supervised loss (i.e. minimising the sum of squared differences between a noisy signal and its noise-free signal estimate) is indistinguishable from the corresponding objective function in conventional fitting.

In contrast, under the Rician noise conditions encountered in MRI acquisition(Gudbjartsson and Patz, 1995), Self-supervised training loss no longer matches conventional fitting. Indeed, the sum of squared errors between noisy signals and noise-free estimates is not an accurate difference metric in the presence of Rician noise.

To summarise: existing supervised DL techniques are associated by high estimation bias, low variance, and end-user flexibility; in contrast, self-supervised methods have lower bias, higher variance, but are limited by the fact their loss is calculated in the signal space $X$.

In light of this, we propose Supervised_MLE, a novel parameter estimation method which combines the advantages of Supervised_GT and Self-supervised methods. This method is contrasted to existing techniques in Fig 1.

This method mimics Self-supervised’s low-bias performance by learning a regularised form of conventional MLE, but does so in the parameter space $Y$, within a supervised learning framework. This addresses the limitations of Self-supervised: Rician noise modelling is incorporated, and parameter loss weighting is not limited by sampling scheme $z$.

Our method learns $\hat{\mathcal{F}}$ by training on noisy signals paired with conventional MLE labels. These labels act as proxies for the groundtruth parameters we wish to estimate:

where $y_{MLE,i}$ is the maximum likelihood estimate associated with the $i^{th}$ training sample.

Our method offers one final advantage over Self-supervised approaches. In addition to the parameter estimation improvements relating to noise model correction and parameter loss weighting, it naturally interfaces with Supervised_GT. In so doing, it presents the opportunity to combine low-bias and low-variance methods into a single, tunable hybrid approach, by a simple weighted sum of each method’s loss function:

The intravoxel incoherent motion (IVIM) model (Le Bihan et al., 1986) was investigated as an exemplar 4-parameter non-linear qMRI model which poses a non-trivial model fitting problem and is well-represented in the DL qMRI literature (Bertleff et al., 2017; Barbieri et al., 2020; Mpt et al., 2021; Mastropietro et al., 2022; Rozowski et al., 2022):

Network architecture was harmonised across all network variants, and represents a common choice in the existing qMRI literature (Barbieri et al., 2020): 3 fully connected hidden layers, each with a number of nodes matching the number of signal samples $z$ (i.e. b-values), and an output layer with a number of nodes matching the number of model parameters. Wider (150 nodes per layer) and deeper (10 hidden layers) networks were investigated and found to have equivalent performance, during both training and testing, at the cost of increased training time. All networks were implemented in Pytorch 1.9.0 with exponential linear unit activation functions (Clevert et al., 2015); ELU performance is similar to ReLU, but is more robust to poor network weight initialisation.

Training datasets were generated at SNR $=[15,30]$ to investigate parameter estimation performance at both high and low noise levels. At each SNR, 100,000 noise-free signals were generated from uniform IVIM parameter distributions ($S_{0}\in[0.8,1.2]$, $f\in[0.1,0.5]$, $D_{slow}\in[0.4,3.0]10^{-3}mm^{2}/s$, $D_{fast}\in[10,150]10^{-3}mm^{2}/s$, representing realistic tissue values), sampling them with a real-world acquisition protocol (Zhao et al., 2015) ($b=[0,10,20,30,50,80,100,200,400,800]$ $s/mm^{2}$), and adding Rician noise. Training data generative parameters were drawn from uniform, rather than in-vivo, parameter distributions to minimise bias in network parameter estimation(Gyori et al., 2022). Data were split 80/20 between training and validation. MLE labels were calculated using a bound-constrained non-linear fitting algorithm, implemented with scipy.optimize.minimize, using either Rician log-likelihood (for Supervised_{MLE, Rician}) or sum of squared errors (for Supervised_{MLE, Gaussian}) as fitting objective function. This algorithm was initialised with groundtruth values (i.e. generative $y$) to improve fitting robustness and avoid local minima. Training/validation samples associated with ‘poor’ MLE labels (defined as lying on the boundary of the bound-constrained estimation space) were held out during training and ignored during validation.

Network training was performed using an Adam optimizer (learning rate = 0.001, betas = (0.9, 0.999), weight decay=0) as follows: Supervised_GT (at SNR 30) was trained 16 times on the same data, each time initialising with different network weights, to improve robustness to local minima during training. From this set of trained networks, a single Supervised_GT network was selected on the basis of validation loss. The trained weights of this selected network were subsequently used to initialise all other networks; in this way, any differences in network performance could be solely attributed to differences in training label selection and training loss formulation. In the case of supervised loss formulations, the inter-parameter weight vector $W$ was chosen as the inverse of each parameter’s mean value over the training set, to obtain equal loss weighting across all four IVIM parameters.

Networks were tested on both synthetic and real qMRI data. The synthetic approach offers (i) known parameter groundtruths to assess estimation against, (ii) arbitrarily large datasets, and (iii) tunable data distributions, but is based on possibly simplified qMRI signals. This approach was used to assess parameter estimation performance in a controlled, rigorous manner; real data was subsequently used to validate the trends observed in silico.

Synthetic data was generated with sampling, parameter distributions, and noise levels matching those used in network training. The IVIM parameter space in which the networks were trained was uniformly sub-divided 10 times in each dimension, to analyse estimation performance as a function of parameter value. At each point in the parameter space, 500 corresponding noisy signals were generated and used to test network performance, accounting for variation under noise repetition.

Real data was acquired from the pelvis of a healthy volunteer, who gave informed consent, on a wide-bore 3.0T clinical system (Ingenia, Philips, Amsterdam, Netherlands), 5 slices, 224 x 224 matrix, voxel size = 1.56 x 1.56 x 7mm, TE = 76ms, TR = 516ms, scan time = 44s per 10 b-values listed in subsection 3.3. For the purposes of assessing parameter estimation methods, we obtained gold standard voxelwise IVIM parameter estimates from a supersampled dataset (16-fold repetition of the above acquisition, within a single scanning session, generating 160 b-values, total scan time = 11m44s). Conventional MLE was performed on this supersampled data to produce best-guess “groundtruth” parameters. During testing, the supersampled dataset was split into 16 distinct 10 b-value acquisitions, each corresponding to a single realistic clinical acquisition. All images were visually confirmed to be free from motion artefacts. The mismatch in parameter distributions between this in-vivo data (highly non-uniform) and the previously-described synthetic data (uniform by construction) limited the scope for validating our in-silico results. To address this, a final synthetic testing dataset was generated from the in-vivo MLE-derived “groundtruth” parameters, and was used for direct comparison between real and simulated data.

Parameter estimation performance was evaluated using 3 key metrics: (i) mean bias with respect to groundtruth, (ii) mean standard deviation under noise repetition, and (iii) root mean squared error (RMSE) with respect to groundtruth. RMSE is the most commonly used metric to evaluate estimation performance (Barbieri et al., 2020; Bertleff et al., 2017), but is limited in its ability to disentangle accuracy and precision; to this end, mean bias and standard deviation were used as more specific measures of network performance.

It is important to note that all methods were assessed with respect to groundtruth qMRI parameters, even those trained on MLE labels. For these methods, the training and validation loss (MLE-based) differed from the reported testing loss (groundtruth-based).

The relative performance of all previously-discussed parameter estimation methods is summarised in Figures 2 and 3. These figures show the bias, variance (represented by its square root: standard deviation), and RMSE of parameter estimates with respect to groundtruth values, reported for each model parameter as a function of its value over the synthetic test dataset; each plotted point represents an average over 500 noise instantiations and a marginalisation over all non-visualised parameters. Marginalisation was required for visualisation of a 4-dimensional parameter space, but was confirmed to be representative of the entire, non-marginalised space, as discussed in §4.5.

In keeping with previously reported results, we show a bias/variance trade-off between different parameter estimation methods. Conventional MLE fitting is provided as a reference (plotted in black). Approaches which, on a theoretical level, approximate conventional MLE (Self-supervised and Supervised_MLE, plotted in red), are generally associated with low bias, high variance, and high RMSE, whereas groundtruth-labelled supervised methods (plotted in blue) exhibit lower variance and RMSE at the cost of increased bias.

Increases in bias, if consistent across parameter space $Y$, do not necessarily reduce sensitivity to differences in underlying tissue properties. However, we show that $Supervised_{GT}$ is associated with bias that varies significantly as a function of groundtruth parameter values. This results in a reduction in information content, visualised as the gradient of the bias plots (top row) in Fig 2. The more negative the gradient, the more parameter estimates are concentrated in the centre of the parameter estimation space $\hat{X}$, and the lower the ability of the method to distinguish differences in tissue properties. This information loss can be seen in Fig 4, which compares $Supervised_{GT}$ to conventional MLE fitting, and shows the compression in $\hat{X}$ over the groundtruth parameter-space $X$.

The above trends, found in simulation, were also observed in real-world data. Fig 5 shows the bias, variance, and RMSE of parameter estimates with respect to “groundtruth” values (obtained from the supersampled dataset described in §3.5). The $x$ axes of these plots correspond to these reference values. To aid visualisation, 10 uniform bins were constructed along each parameter dimension, into which clinical voxels were assigned based on their “groundtruth” parameter values. Fig 5 plots the mean bias, standard deviation, and RMSE associated with each bin as a function of the bin’s central value, together with the distribution of voxels across the 10 bins.

By calculating the variance of the 16 $b=0$ images, the SNR of this clinical dataset was found to be $\sim$15; Fig 2 is therefore the relevant point of comparison. It can be readily seen that the trends observed in simulated data, described in §4.1, are replicated for $f<0.40$, $D_{slow}<1.5$, and the entire range of $D_{fast}$, namely the regions of parameter-space which are well-represented in the real-world data. Fig 6 confirms that divergence outside of these ranges is due to under-representation in the in vivo test data; the apparent divergences can be replicated in-silico by matching real-world parameter distributions.

Fig 7 contains exemplar parameter maps from the clinical test data, and shows the real-world implications of the trends summarised in Figures 2 and 5: $Supervised_{GT}$’s low-variance, low-RMSE parameter estimation results in artificially smooth IVIM maps biased towards mean parameter values.

Our proposed method occupies the low-bias side of the bias-variance trade-off discussed in §3.5, and offers four broad advantages over the competing method in this space (Self-supervised): (i) flexibility in choosing inter-parameter loss weighting $W$, (ii) incorporation of non-Gaussian (e.g. Rician) noise models, (iii) compatibility with complex, non-differentiable signal models $M$, and (iv) ability to interface with low-variance methods, to produce a hybrid approach tunable to the needs of the task at hand. These advantages are analysed in turn.

We note that RMSE is a poor summary measure of network performance. RMSE is heavily skewed by outliers, and thus favours methods which give parameter estimates consistently close to mean parameter values. Such estimates, as in the case of $D_{fast}$, may contain very little information (Fig 4) despite being associated with low RMSE. Accordingly, we strongly recommend that RMSE be discontinued as a single summary metric for parameter estimation performance: it must always be accompanied by bias, variance, and ideally an analysis of information content.

RMSE’s limitations as a performance metric during testing may also call into question its suitability as a loss metric during training. This work, much like the rest of the DL qMRI literature, employs a training loss (MSE, described in Sections 2.2 and 2.3) which is monotonically related to RMSE. Whilst outside the scope of this work, implementing a non-RMSE-derived training loss (such as mean absolute error) may be worth of future investigation.

The above analysis has been largely based on Figs 2 and 3, which show parameter estimation performance marginalised over 3 dimensions of $X$. This choice, made to aid visualisation, was validated against higher dimensional representations of the same data.

Fig 11 compares Supervised_{MLE, Rician} and Supervised_Groundtruth performance across the entire qMRI parameter space. It can be seen that trends observed in Fig 2 are replicated here; we draw attention to two such examples. Firstly, Fig 2 suggests Supervised_Groundtruth produces lower $f$ standard deviation than Supervised_{MLE, Rician}; Fig 11 confirms this to be the case across all test data. In contrast, Fig 2 suggests that Supervised_Groundtruth produces higher $D_{slow}$ bias at low $D_{slow}$ and lower bias at high $D_{slow}$; Fig 11 confirms a spread of bias differences across the test data: some favouring one method, and others the other. This effect is explored in Fig 12, which compares $D_{slow}$ estimation performance as a function of $f$ and $D_{fast}$ at two specific (non-marginalised) groundtruth $D_{slow}$ values ($0.69,2.71)$. As expected from the marginalised representation in Fig 2, at low $D_{slow}$ Supervised_Groundtruth produces higher bias across the entire $f$-$D_{fast}$ parameter space, whereas at high $D_{slow}$ the opposite is true.

Despite this, it is important to note the limitations of marginalisation. Fig 12 also shows that the relative performance of Supervised_{MLE, Rician} and Supervised_Groundtruth varies across all parameter-space dimensions. Consider $D_{slow}=0.69$, where 2 shows similar marginalised RMSE for these methods. In fact, by visualising this difference as a function of $f$ and $D_{fast}$, we reveal two distinct regions: high $f$/low $D_{fast}$ (where Supervised_{MLE, Rician} produces lower RMSE), and elsewhere (where it produces higher RMSE). This highlights (i) the potential pitfalls of producing summary results by marginalising across entire parameter spaces and (ii) the need to choose parameter-estimation methods appropriate for the specific parameter combinations relevant to the tissues being investigated (Epstein et al., 2021).

This work has focused on voxelwise DL parameter estimation methods: networks which map one signal curve to its corresponding parameter estimate. There are, however, alternatives: convolutional neural network methods which map spatially related clusters (“patches”) of qMRI signals to corresponding clusters of parameter estimates (Fang et al., 2017; Ulas et al., 2019; Li et al., 2022). Our MLE training label approach could be incorporated into such methods, and we leave it to future work to investigate the effect this would have on parameter estimation performance.