Model-based multi-parameter mapping

Abstract

Quantitative MR imaging is increasingly favoured for its richer information content and standardised measures. However, computing quantitative parameter maps, such as those encoding longitudinal relaxation rate (R₁), apparent transverse relaxation rate (R₂^*) or magnetisation-transfer saturation (MT_sat), involves inverting a highly non-linear function. Many methods for deriving parameter maps assume perfect measurements and do not consider how noise is propagated through the estimation procedure, resulting in needlessly noisy maps. Instead, we propose a probabilistic generative (forward) model of the entire dataset, which is formulated and inverted to jointly recover (log) parameter maps with a well-defined probabilistic interpretation (e.g., maximum likelihood or maximum a posteriori). The second order optimisation we propose for model fitting achieves rapid and stable convergence thanks to a novel approximate Hessian. We demonstrate the utility of our flexible framework in the context of recovering more accurate maps from data acquired using the popular multi-parameter mapping protocol. We also show how to incorporate a joint total variation prior to further decrease the noise in the maps, noting that the probabilistic formulation allows the uncertainty on the recovered parameter maps to be estimated. Our implementation uses a PyTorch backend and benefits from GPU acceleration. It is available at https://github.com/balbasty/nitorch.

Graphical abstract

Fig. 1

Qualitative analysis of two toy problems related to MR relaxometry. The first toy problem (first row) has an objective function of the form ${(x-\exp (-y\left)\right)}^2$, while the other toy problem (second row) as one of the form ${(x-\exp (-\exp \left(y\right)\left)\right)}^2$. In each case, the first column shows the objective function in black along with quadratic approximations corresponding to Gauss-Newton (red) and the proposed second order method (blue). The second column shows the magnitude of the gradient-hessian ratio. The black line corresponds to a purely quadratic objective function, whose optimum is found in one step ($\left|g\left(x\right)\right|/h\left(x\right)=|x-x^{*}|$). When the ratio is under this black line, monotonic convergence is ensured. The third column shows convergence of the Gauss-Newton (red) and proposed (blue) iterative schemes when attacking from either a convex or non-convex lobe, with the optimum shown in dashed black. In each case, convergence is shown both in terms of objective value per iteration and parameter value per iteration. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2

Qualitative analysis of two-dimensional toy problems related to MR relaxometry. The first toy problem (top half) has an objective function of the form $\ell (y,z)={(x_0-e^{z-t_{0}y})}^2+{(x_1-e^{z-t_{1}y})}^2$, while the other toy problem (bottom half) has one of the form $\ell (y,z)={(x_0-e^{z-t_0{e}^y})}^2+{(x_1-e^{z-t_1{e}^y})}^2$. In each case, the first column shows the objective function along each of the two dimensions (in black) along with quadratic approximations corresponding to Gauss-Newton (red) and the proposed second order methods (blue and green). The second column shows the normalised step size for each method, color-coded such that values below one are blue and values above one are red. The optimum is marked by a red star. Convergence paths from two random initial coordinates (solid and dashed curves) are overlaid. In the second toy problem, two different loading matrices are used: Proposed uses the diagonal of the absolute Hessian $\mathrm{diag}\left(\left|\mathit{H}\right|\right)$ whereas Proposed+ uses its sum across columns $\left|\mathit{H}\right|\mathit{1}$. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3

Negative log-likelihood and step size function of the nonlinear SPGR model. The first column shows the log-likelihood of a single observation with respect to each log-parameter, while keeping the others fixed. The second column shows the step size function (red) obtained with the proposed approximate Hessian, compared with the theoretical perfect step size function $f\left(x\right)=|x-x^{*}|$. As long as the actual step size is below this theoretical function, monotonic convergence is ensured. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4

Optimisation of 1000 simulated voxels for 10,000 iterations using the proposed second order method (left), Gauss-Newton (middle) and Levenberg-Marquardt (right). The top graph shows the negative log-likelihood per iteration for 50 random samples, while the middle graph shows the normalised gain per iteration for these same 50 samples. Gain is only shown for the proposed method. Note that if an iteration was non-monotonic, its gain would be negative, which never happens in our simulations. The bottom graph shows the distribution of the final log-likelihood of all 1000 samples. The dashed line represents the expected log-likelihood of the true parameters. Log-likelihood and iterations are shown in log scale.

Fig. 5

Optimisation of 1000 simulated voxels for 10,000 iterations using the proposed second order method. The optimisation is either performed in the basis of the log-parameters (left), the rate representation (middle) or the time representation (right). The top graph shows the negative log-likelihood per iteration for 50 random samples, while the middle graph shows the normalised gain per iteration for these same 50 samples. The bottom graph shows the distribution of the final log-likelihood of all 1000 samples. The dashed line represents the expected log-likelihood of the true parameters. Log-likelihood and iterations are shown in log scale.

Fig. 6

Top: $R_1$, $R_2^{*}$ and MT_sat maps obtained with all methods. Bottom: z-normalised mean squared error between (unseen) predicted and acquired gradient-echo images. MSE was computed across all voxels in a single volume. Lower is better.

Fig. 7

Decimation of the number of echoes. Maps were computed using the first 2, 4 or 6 echoes. Graphs on the left display the mean squared error (lower is better) relative to the reference maps (SPGR-ML with four repeats of 8 echoes) as a function of the number of echoes used. The corresponding parameter maps are displayed on the right.

Fig. 8

Uncertainty mapping and propagation. Log parameter maps ($\mathbb{E}\left[\log r_1\right]={\tilde{r}}_1$) and their uncertainty ($\mathrm{sd}\left[\log r_1\right]$) were estimated using SPGR+JTV, which includes an edge-preserving prior and SPGR-ML, which has no (or infinitely uninformative) prior. The expected value and standard deviation of $r_1$ and $t_1$ under a Laplace posterior ($\mathbb{E}\left[r_1\right]$) or a peaked posterior ($\exp \left({\tilde{r}}_1\right)$) are shown.