Scientists interested in technical advances in the measurement of longitudinal relaxation time, T₁, and in applications of Bayesian analysis.Several studies have shown the importance of T₁ determination for the characterization of neurological diseases [1-2]. Using the spoiled gradient recalled echo (SPGR) sequence with very short repetition time, TR_SPGR, it has been shown that whole brain T₁-mapping can be obtained within clinically realistic times [3-4]. However, SPGR images suffer from limited signal-to-noise ratio (SNR), which may limit the applicability of this technique to high resolution T₁-mapping. Here, we introduce two Bayesian-based analyses that use spatial information as a prior to improve the quality of voxel-by-voxel T₁-mapping from SPGR data, and compare the results with those derived using a conventional nonlinear least-squares (NLLS)-based algorithm.

SPGR signal model

The mono-component SPGR signal is given by

$$M_{SPGR}^{k}=M_{SPGR}^{0}\sin\alpha\left(\begin{array}{c}\frac{1-E_{1}}{1-E_{1}\cos\alpha}\end{array}\right), [1]$$

where $$$M_{SPGR}^{0}$$$ represents the signal amplitude at zero echo-time, α_k is the k^th excitation FA out of a total of K FAs, and $$$E_{1}=\exp\left(\begin{array}{c}-\frac{TR_{SPGR}}{T_{1}}\end{array}\right)$$$. Here, FAs were assumed well calibrated.

Bayesian Analysis with Spatial Information Collaboration (BASIC) for T₁ mapping

Marginalized likelihood function (MLF): The LF for the SPGR dataset is given by

$$L\left(\begin{array}{c}\mathbf{S}\ \mid\ \mathbf{\Lambda}\end{array}\right)=\left(\begin{array}{c}\sigma_{SPGR}\sqrt{2\pi}\end{array}\right)^{-K}\exp\left(\begin{array}{c}-\frac{\left(\begin{array}{c}\mathbf{S_\mathit{SPGR}}-\mathbf{M_\mathit{SPGR}}\end{array}\right)\left(\begin{array}{c}\mathbf{S_\mathit{SPGR}}-\mathbf{M_\mathit{SPGR}}\end{array}\right)^T}{2\sigma_{SPGR}^2}\end{array}\right), [2]$$

where $$$ \mathbf{\Lambda}=\left(\begin{array}{c}T_{1}\ M_{SPGR}^0\ \sigma_{SPGR}\end{array}\right)$$$, M_SPGR and S_SPGR are the theoretical and measured signal values at different FAs, respectively, and σ_SPGR is the standard deviation (SD) of noise. Assuming noninformative normal-inverse-Gamma distributions on $$$M_{SPGR}^0$$$ and $$$\sigma_{SPGR}$$$, the MLF is given by [5]

$$L\left(\begin{array}{c}\mathbf{S}\ \mid\ {T_{1}}\end{array}\right)=C\left(\begin{array}{c}\left(\begin{array}{c}\mathbf{S_\mathit{SPGR}}\ \mathbf{S_\mathit{SPGR}}^T\end{array}\right)-\frac{\left(\begin{array}{c}\mathbf{g_\mathit{SPGR}}\ \mathbf{S_\mathit{SPGR}}^T\end{array}\right)^{2}}{\left(\begin{array}{c}\mathbf{g_\mathit{SPGR}}\ \mathbf{g_\mathit{SPGR}}^T\end{array}\right)}\end{array}\right)^{-\frac{K}{2}}, [3]$$

which is independent of $$$M_{SPGR}^0$$$ and $$$\sigma_{SPGR}$$$, where C is a proportionality constant independent of T₁ , and $$$\mathbf{g}=\sin\boldsymbol\alpha\left(\begin{array}{c}\frac{1-E_{1}}{1-E_{1}\cos\boldsymbol\alpha}\end{array}\right)$$$.

Spatial prior distribution: The prior is a Gaussian distribution given by

$$P\left(\begin{array}{c}T_{{1}_{i}}\ \mid\ \mu_{{T}_{1}},\sigma_{{T}_{1}}\end{array}\right)=\left(\begin{array}{c}\sigma_{{T}_{1}}\sqrt{2\pi}\end{array}\right)^{-1}\exp\left(\begin{array}{c}-\frac{\left(\begin{array}{c}T_{{1}_{i}}-\mu_{{T}_{1}}\end{array}\right)^{2}}{2\sigma_{T_{1}}^2}\end{array}\right)P\left(\begin{array}{c}\mu_{{T}_{1}},\sigma_{{T}_{1}}\end{array}\right), [4]$$

where $$$T_{{1}_{i}}$$$ is the $$$T_{1}$$$ value for a voxel i, and $$$\mu_{{T}_{1}}$$$ and $$$\sigma_{{T}_{1}}$$$ are, respectively, the mean and SD of $$$T_{1}$$$ calculated from M pixels of a region-of-interest (ROI) around i. $$$P\left(\begin{array}{c}\mu_{{T}_{1}},\sigma_{{T}_{1}}\end{array}\right)=\sigma_{{T}_{1}}^{-1}$$$ is a hyper-prior on $$$\mu_{{T}_{1}}$$$ and $$$\sigma_{{T}_{1}}$$$ taken as a noninformative Jeffrey’s prior.

Posterior probability distribution:

Known $$$\mu_{{T}_{1}}$$$ and $$$\sigma_{{T}_{1}}$$$

$$$\mu_{{T}_{1}}$$$ and $$$\sigma_{{T}_{1}}$$$ can be estimated from voxel-by-voxel NLLS fit. In this case, the posterior distribution of T₁ is given by

$$P\left(\begin{array}{c}T_{{1}_{i}}\ \mid\ \mathbf{S_{\mathit{i}}},\mu_{{T}_{1}},\sigma_{{T}_{1}}\end{array}\right)\propto L\left(\begin{array}{c}\mathbf{S_{\mathit{i}}}\ \mid\ T_{{1}_{i}}\end{array}\right)P\left(\begin{array}{c}T_{{1}_{i}}\ \mid\ \mu_{{T}_{1}},\sigma_{{T}_{1}}\end{array}\right), [5]$$

from which the estimate $$$\widehat{T}_{{1}_{i}}$$$ can be derived as the posterior mean through

$$\widehat{T}_{{1}_{i}}=\int T_{{1}_{i}}\ P\left(\begin{array}{c}T_{{1}_{i}}\ \mid\ \mathbf{S_{\mathit{i}}},\mu_{{T}_{1}},\sigma_{{T}_{1}}\end{array}\right)dT_{{1}_{i}}. [6]$$

We will refer to this method as BASIC₁-T₁.

Unknown $$$\mu_{{T}_{1}}$$$ and $$$\sigma_{{T}_{1}}$$$

If $$$\mu_{{T}_{1}}$$$ and $$$\sigma_{{T}_{1}}$$$ are unknown, T₁'s of all M voxels of the ROI and $$$\mu_{{T}_{1}}$$$ and $$$\sigma_{{T}_{1}}$$$ become nuisance parameters. In this case, the posterior distribution for the ROI is given by

$$P\left(\begin{array}{c}T_{{1}_{1:M}},\mu_{{T}_{1}},\sigma_{{T}_{1}}\ \mid\ \mathbf{S_{\mathit{1:M}}}\end{array}\right)\propto \prod_{j=1}^M L\left(\begin{array}{c}\mathbf{S_{\mathit{j}}}\ \mid\ T_{{1}_{j}}\end{array}\right)P\left(\begin{array}{c}T_{{1}_{j}}\ \mid\mu_{{T}_{1}},\sigma_{{T}_{1}}\end{array}\right), [7]$$

After marginalization over $$$\sigma_{{T}_{1}}$$$ of the prior distribution given in Equation 4, Equation 7 becomes

$$P\left(\begin{array}{c}T_{{1}_{1:M}},\mu_{{T}_{1}}\ \mid\ \mathbf{S_{\mathit{1:M}}}\end{array}\right)\propto \prod_{j=1}^M L\left(\begin{array}{c}\mathbf{S_{\mathit{j}}}\ \mid\ T_{{1}_{j}}\end{array}\right)P\left(\begin{array}{c}T_{{1}_{j}}\ \mid\mu_{{T}_{1}}\end{array}\right), [8]$$

where $$$P\left(\begin{array}{c}T_{{1}_{j}}\ \mid\mu_{{T}_{1}}\end{array}\right)=\frac{\Gamma(1/2)}{2\sqrt{\pi}\mid T_{{1}_{j}}-\mu_{{T}_{1}}\mid}$$$, and $$$\Gamma(.)$$$ is the gamma function [5].

Finally, $$$\widehat{T}_{{1}_{i}}$$$ can be derived through

$$\widehat{T}_{{1}_{i}}=\int ...\int {T}_{{1}_{i}}\ P\left(\begin{array}{c}T_{{1}_{1:M}},\mu_{{T}_{1}}\ \mid\ \mathbf{S_{\mathit{1:M}}}\end{array}\right)d{T}_{{1}_{1:M}}\ d\mu_{{T}_{1}}. [9]$$

Given the high-dimensional nature of Equation 9, we used a Markov chain Monte Carlo-based algorithm to calculate the integral [6]. We will refer to this method as BASIC₂-T₁.

The performance of BASIC_1,2-T₁ was evaluated on T₁-weighted images generated at FAs of 4^o, 8^o and 18^o with TR_SPGR =8ms from a numerical phantom constructed with different T₁ values (Fig.1). Analysis was performed at $$$SNR=M_{SPGR}^{0}/\sigma_{SPGR}$$$ of 100, 200 and 500. Results were compared to those obtained with stochastic region contraction (SRC), an NLLS-based algorithm that does now require specific initial parameter estimates [7]. For the quantitative comparison, relative error maps were calculated. To verify that the performance of BASIC_1,2-T₁ cannot be obtained by simple image filtering, T₁-maps were calculated from smoothed T₁-weighted images with the mean filter of size 5x5 corresponding to the ROI size used in BASIC_1,2-T₁.Fig.1 shows T₁ parameter and error maps derived using SRC-NLLS and BASIC_1,2-T₁ analyses. At high SNR, all methods performed well. At low-to-moderate SNRs, SRC-NLLS results show substantial random variation in estimated T₁, while BASIC_1,2-T₁ produced much higher quality maps and lower errors. Although T₁-maps derived from smoothed images showed higher homogeneity in homogeneous regions, large errors, with blurring and detail loss, were clearly seen at edges and within small structures. Our approach is based on the fact that BASIC analyses combine voxel-by-voxel fitting with ROI parameter estimation. ROI parameters act as a constraint, while voxel fitting mitigates blurring and detail loss. Estimation of T₁ from SPGR imaging data was markedly improved through use of BASIC analysis. We are currently extending this work to T₂ determination from balanced steady-state imaging data [3] and in-vivo validation.