3D multi-echo FLASH scans, acquired at a selection of flip angles, have been widely presented for estimating $$$T_1$$$ at every voxel^1-7. These voxel-by-voxel point estimates of $$$T_1$$$ can then be used to generate $$$T_1$$$ maps. However, these point estimates can be misleading for two reasons. First, given the standard signal model, even the "best" estimate of $$$T_1$$$ for the voxel may still be a very bad fit to the data if the model fails to describe the underlying physical processes in the voxel. Second, it is difficult to determine how well other values of $$$T_1$$$, beyond just the single point estimate, might fit the observed data. To address these issues, we present a practical algorithm that, given raw multi-channel, multi-echo, multi-flip FLASH data, produces both the maximum likelihood estimate (MLE) of $$$T_1$$$ and a confidence interval (CI) at each voxel.

3D FLASH scans were acquired on one organic pineapple using a 1.5 T Avanto (Siemens Healthcare, Erlangen, Germany) and the product 12-channel headcoil. All scans had 1 mm isotropic resolution, 256x256x176 matrix, 4 echoes with TEs of 1.85, 3.77, 5.69 and 7.61 ms, bandwidth of 650 Hz/px, TR of 11 ms, and a duration of 8:16. Four scans were acquired, with flip angles of 5, 10, 20, and 30 degrees. An additional 11 s of scanning with no RF pulses was used to acquire pure noise measurements.

The channel covariance matrix was estimated from the noise measurements. Complex-valued images were reconstructed from each channel and each echo of each FLASH scan and pre-whitened using the covariance. The ML estimate $$$\widehat{T_1}$$$ is the minimizer of

$$\epsilon(T_1) = \min_{\vec{c}, \vec{p}} \sum_{i,j,k} \left|w_{i,j,k} - c_i p_j q_k(T_1) \right|^2$$

where $$$\vec{c}$$$ represents the complex coil sensitivites, $$$\vec{p}$$$ are the complex echo-specific effects (e.g., $$$T_2^*$$$ decay and off-resonance) at this voxel, and $$$q(T_1)$$$ is the FLASH steady state equation. Evaluating $$$\epsilon(T_1)$$$ can be shown to be equivalent to solving for the maximum eigenvalue of an ExE matrix with E the number of echoes.

A very close, but more efficient, approximation solves this in two steps. First, create a CxEF matrix from $$$w_{i,j,k}$$$ where C and F are the number of channels and flips. The first right singular vector of this matrix, multiplied by the first singular value, is an estimate of the "optimal" channel combination. Splitting this vector into an ExF matrix $$$W$$$ gives the new error equation

$$\epsilon'(T_1) = \left| W \right|_F - \frac{\left|W q(T_1) \right|_2^2}{\left| q(T_1) \right|_2^2}$$

The MLE of $$$T_1$$$ is closely approximated by the minimizer of $$$\epsilon'$$$. Empirically, we find that $$$\epsilon'$$$ is monotonic over a wide range of $$$T_1$$$, and so minimize via Brent's Method⁸. Knowing the noise distribution, we can estimate $$$p(\epsilon'(T_1) | T_1)$$$ (the probability of the observed error being pure noise) via Monte Carlo simulation. CIs are then defined via $$$p(\epsilon'(T_1) | T_1) < \alpha$$$. Each $$$\alpha$$$ defines a threshold $$$t$$$ for $$$\epsilon'$$$, and upper/lower ends of the CI are defined as the $$$T_1$$$ where $$$\epsilon'(T_1) - t = 0$$$ found searching upwards/downwards from the MLE⁸. Note if $$$\epsilon'(\widehat{T_1}) > t$$$, the CI is empty. Additionally, the norm of the measurements is compared against the null hypothesis that the voxel is pure noise, and MLE/CI fitting is only performed if this is rejected at $$$p=0.05$$$.

Figure 1 shows the MLE presented as a map image. Figure 2 shows line projections of MLE, upper/lower ends of the 95% CI, voxels flagged as noise, and voxels for which the CI is empty despite having significant signal. We note that regions where the CI is empty are likely due to the model describing the data poorly – there is a "best" $$$T_1$$$ (the MLE), but for any $$$T_1$$$ the residual is larger than can be explained by noise. Empty CIs provide an important diagnostic for the reliability of the MLE. Similarly, wide CIs provide an indication of uncertainty in the MLE. Unfortunately narrow CIs do not necessarily indicate precision; the CI is only expected to cover the true $$$T_1$$$ with $$$p=\alpha$$$.

Fitting these values took 34 minutes using 8 cores of a server with Intel Xeon E5-2695 CPUs, or roughly 176 μs per voxel. The work could be more highly parallelized for further acceleration.

We have presented a novel algorithm for efficiently computing the MLE and a CI for $$$T_1$$$ given multi-flip, multi-echo data. Our results show that the CI can be used to help interpret the validity of the MLE, particularly detecting problems with model specification, even when the MLE $$$T_1$$$ map appears reasonable.