Joint T1 and T2 Mapping with Tiny Dictionaries and Subspace-Constrained Reconstruction

TL;DR

This paper introduces an adaptive method for joint T1-T2 mapping that generates tiny dictionaries by approximating the Bloch-response manifold, enabling efficient subspace reconstruction with minimal loss of accuracy.

Contribution

It proposes a novel approach that decouples dictionary size from accuracy by adaptively refining sampling grids based on local approximation error.

Findings

Excellent agreement with traditional template matching

Reduces dictionary sizes by one to two orders of magnitude

Effective in phantom and in vivo studies

Abstract

Purpose: To develop a method that adaptively generates tiny dictionaries for joint T1-T2 mapping. Theory: This work breaks the bond between dictionary size and representation accuracy (i) by approximating the Bloch-response manifold by piece-wise linear functions and (ii) by adaptively refining the sampling grid depending on the locally-linear approximation error. Methods: Data acquisition was accomplished with use of an 2D radially sampled Inversion-Recovery Hybrid-State Free Precession sequence. Adaptive dictionaries are generated with different error tolerances and compared to a heuristically designed dictionary. Based on simulation results, tiny dictionaries were used for T1-T2 mapping in phantom and in vivo studies. Reconstruction and parameter mapping were performed entirely in subspace. Results: All experiments demonstrated excellent agreement between the proposed mapping technique and template matching using heuristic dictionaries. Conclusion: Adaptive dictionaries in combination with manifold projection allow to reduce the necessary dictionary sizes by one to two orders of magnitude.

Equations

$$E_{tn}=\sum _j{J}_{tj}(X_{jn}-x_j^{\text{ref}})-(Y_{tn}-y_t^{\text{ref}})$$

$$n^{\ast }=\arg \underset{n}{\max }\sqrt{\sum _t{E}_{tn}^2}$$

$$j^{\ast }=\arg \underset{j}{\max }\left|(X_{j{n}^{\ast }}-x_j^{\text{ref}})/x_j^{\text{ref}}\right|.$$

$$X_{\bullet 1}$$

$$\text{ and }{X}_{\bullet 2}$$

$$y^{\text{ref}}$$

$$\{\hat{x},\hat{\rho }\}$$

$$=\arg \underset{x,\rho }{\min }{\Vert A(x,{\rho }^{-1}{)}^{\top }-J{x}^{\text{ref}}+y^{\text{ref}}\Vert }_2^2,\text{where }A:=[J|-y]$$

$${\alpha }^{\ast }=\arg \underset{\alpha }{\min }{\Vert y-{\mathcal{P}}_{\vec{k}}\mathcal{F}S{\Phi }_K\alpha \Vert }_2^2+\lambda R(\alpha )$$

Full text

Joint T1 and T2 Mapping

with Tiny Dictionaries

and Subspace-Constrained Reconstruction

Volkert Roeloffs

Biomedizinische NMR, MPI für biophysikalische Chemie, 37070 Göttingen, Germany

Martin Uecker

Institute for Diagnostic and Interventional Radiology, University Medical Center, 37075 Göttingen, Germany

German Centre for Cardiovascular Research (DZHK), partner site Göttingen, Germany

Jens Frahm

Biomedizinische NMR, MPI für biophysikalische Chemie, 37070 Göttingen, Germany

German Centre for Cardiovascular Research (DZHK), partner site Göttingen, Germany

Running head: T1 and T2 Mapping with Tiny Dictionaries

Correspondence to:

Dr. V. Roeloffs

Biomedizinische NMR, MPI für biophysikalische Chemie

37070 Göttingen, Germany

volkert.roeloffs@mpibpc.mpg.de

Approximate word count: 148 (abstract) 2900 (body)

Number of pages: 14

Number of figures: 5

Number of tables: 0

Submitted to Magnetic Resonance in Medicine as a Note

Date of submission: 2018-12-22

Date of revision: XXX

Abstract

Purpose: To develop a method that adaptively generates tiny dictionaries for joint $T_{1}$-$T_{2}$ mapping.

Theory: This work breaks the bond between dictionary size and representation accuracy (i) by approximating the Bloch-response manifold by piece-wise linear functions and (ii) by adaptively refining the sampling grid depending on the locally-linear approximation error.

Methods: Data acquisition was accomplished with use of an 2D radially sampled Inversion-Recovery Hybrid-State Free Precession sequence. Adaptive dictionaries are generated with different error tolerances and compared to a heuristically designed dictionary. Based on simulation results, tiny dictionaries were used for $T_{1}$-$T_{2}$ mapping in phantom and in vivo studies. Reconstruction and parameter mapping were performed entirely in subspace.

Results: All experiments demonstrated excellent agreement between the proposed mapping technique and template matching using heuristic dictionaries.

Conclusion: Adaptive dictionaries in combination with manifold projection allow to reduce the necessary dictionary sizes by one to two orders of magnitude.

Key words: subspace reconstruction, multi-parametric mapping, T1 mapping, T2 mapping, model-based, dictionary, quantitative MRI

Introduction

Multi-parametric mapping of MRI-detectable physical or physiological quantities has the potential to detect subtle abnormalities earlier and in a more objective manner than conventional contrast-weighted imaging. However, traditional methods that acquire a set of fully sampled images with varying contrasts and then perform a pixelwise fitting are typically very time-consuming. Model-based methods accelerate the measurement by estimating the quantitative maps directly from undersampled k-space data and remove the need to acquire fully sampled images [1, 2, 3, 4, 5, 6].

Recently, new approaches have been presented that break with simple signal models and employ more sophisticated excitation patterns [7, 8, 9, 10, 11, 12]. One way to deal with the resulting complex signal responses is to generate a bank of signal prototypes or “dictionaries” [2, 7]. However, these dictionaries (i) are typically very large in size, (ii) scale exponentially with the number of parameters, (iii) take long to compute, and (iv) result in a huge number of comparisons at the stage of matching. A variety of ideas have been presented to overcome the associated difficulties, as a simple reduction of the sampling density would result in a reduction of representation accuracy. For example, dictionaries compressed by singular value decomposition (SVD) exploit redundancies to perform the matching process in a reduced-dimensional space [13, 14]. With the use of clustering properties [15] matching can further be sped up as unnecessary comparisons are avoided. Both approaches rely on dictionaries, in which, first of all, sampling positions have been chosen heuristically.

Here, we present a new approach to automatically generate sampling positions in an adaptive way. These positions are then considered a set of support points that approximate the Bloch-response manifold [16] by piece-wise linear functions. Manifold projection in combination with these adaptively designed dictionaries allows reduction of necessary dictionary sizes by one to two orders of magnitude. The new method is applied to accomplish joint $T_{1}$ and $T_{2}$ mapping.

Theory

The MRI signal response to a complex excitation pattern is given by the Bloch equations [17]. If the excitation sequence is sufficiently rich and the signal response sensitive to the parameter of interest, all signal responses lie on a non-linear smooth manifold that is embedded within the higher-dimensional (time-domain) space. The low-dimensional manifold is called the Bloch-response manifold [16] and here used in two ways to break the bond between dictionary size and representation accuracy: First, we approximate the Bloch-response manifold by piece-wise linear functions and consider the dictionary a set of support points. As a consequence, mapping to the parametric domain becomes continuous rather than discretized by the chosen sampling grid. Second, we allow the sampling grid to be refined adaptively during the generation of the dictionary depending on the precision needed. To this end, an initial grid is recursively refined in regions where the locally-linear approximation is not accurate enough.

Piece-wise linear approximation and adaptive sampling

The basic idea of an adaptively refined dictionary generation is to allow coarse sampling in regions with locally-linear signal dependency and fine sampling where non-linear dependencies are present. More specifically, in the vicinity of a reference position $x^{\text{ref}}=(T_{1},T_{2})^{\top}$ in $T_{1}$-$T_{2}$ parameter space, the locally-linear approximation $y-y^{\text{ref}}\approx J(x-x^{\text{ref}})$ holds, where $J$ is the Jacobian matrix (defining the best linear approximation of the nonlinear map at position $x^{\text{ref}}$), $y^{\text{ref}}$ the signal response at the reference position, and $x$ and $y$ neighboring vectors in parameter domain and temporal domain, respectively. Considering neighborhoods $X$ and $Y$, i.e. matrices with columns being neighboring vectors in parameter domain and temporal domain, respectively, approximation errors

[TABLE]

can be defined for each time point $t$ for each vector $n$ in these neighborhoods. The closer the neighboring vectors are to the reference position, the smaller the approximation errors become (smoothness of the Bloch-response manifold).

This motivates our proposed strategy for adaptive dictionary generation: Starting with an initial neighborhood, the entire dictionary can be built by recursively splitting into downsized neighborhoods until the total approximation error $E^{\text{tot}}=\sqrt{\sum_{t,n}E_{tn}^{2}}$ fulfills the stopping criterion $E^{\text{tot}}~{}<~{}\varepsilon~{}\sqrt{\sum_{t,n}Y_{tn}^{2}}$ which is controlled by the predefined error tolerance $\varepsilon$.

While downsized neighborhoods could be generated in different ways (quadtree- or binary space partitioning, isotropic downscaling of neighborhoods, etc.), we propose to split the current neighborhood only into a single direction at a time.

More specifically, we identify the index $n^{*}$ of the neighbor exhibiting the largest root-sum-squares error

[TABLE]

and select a coordinate axes as split direction by computing the largest relative parameter deviation with respect to the reference position

[TABLE]

The current neighborhood is then split into the direction of the $j^{*}$-th coordinate axis yielding two half-sized neighborhoods both subject to recursive splitting.

The final result of the recursive building process is a lookup table linking $N$ model signals in the dictionary $D_{i}\in\{1,2,...,N\}$ to their corresponding position in parameter space. Here, for the proposed manifold projection, the Jacobian matrix is stored additionally for each position.

Generating neighborhoods

For joint $T_{1}$-$T_{2}$ mapping, the embedded Bloch-response manifold is two-dimensional. Consequently, a minimum of two neighbors have to be generated for a new reference position in parameter domain. Here, these two neighbors are generated according to

[TABLE]

The integers $p$ and $q$ reflect the recursion depths in $T_{1}$ and $T_{2}$ splitting direction, respectively, and are increased as long as the approximation error exceeds the prescribed threshold.

Manifold projection

As the final approximation error for all neighborhoods is smaller than the error tolerance, it is guarantied that the non-linear Bloch-response in the vicinity of each entry $D_{i}$ can linearly be approximated by the respective Jacobian matrix. Consequently, the manifold projection is realized in two steps. First, an appropriate entry $y^{\text{ref}}$ in the dictionary has to be identified which is realized by pattern matching similar to Refs. [7, 18, 14]:

[TABLE]

Second, the linear function in the region around this reference position has to be inverted to project the reconstructed response signal $y$ onto the piece-wise linear manifold in the parametric domain. The final projection can be cast into a least-squares optimization problem of the form

[TABLE]

and can be solved by the Moore-Penrose pseudo inverse $A^{+}$ to yield the final quantitative result $(\hat{T_{1}},\hat{T_{2}},\hat{\rho}^{-1})^{\top}=A^{+}(Jx^{\text{ref}}-y^{\text{ref}})$ which includes the proton density $\rho$. This scaling constant is added as an additional unknown as all entries in the dictionary have been generated with unit proton density.

Methods

Heuristic design and template matching

To evaluate accuracy and size of the adaptively generated dictionaries, the heuristically designed dictionary in Ref. [7] is chosen as a benchmark. Ma and coworkers partitioned the $T_{1}$-$T_{2}$ space into 4 regions with different sampling densities (see Figure 1E). The parametric representation of a signal response is found by identifying the best matching entry in the dictionary $D$ (template matching) and assigning the corresponding parameter values from the lookup table to this pixel. This procedure is identical to the first step of our proposed manifold projection.

For a meaningful comparison, all adaptive dictionaries share the boundary conditions of Ref. [7], namely $T_{1}\in[$0.1\text{,}\mathrm{s}$,$5\text{,}\mathrm{s}$]$, $T_{2}\in[$0.02\text{,}\mathrm{s}$,$3\text{,}\mathrm{s}$]$, and the physical constraint $T_{1}\geq T_{2}$. Here, offsets in the $B_{0}$ field are excluded.

MRI

All MRI studies were performed at a field strength of $3\text{\,}\mathrm{T}$ (Magnetom Prisma, Siemens Healthineers, Erlangen, Germany) using a 64-channel head coil. Volunteers without known illness were recruited and written informed consent was obtained before MRI according to the regulations of the local ethics committee.

Data acquisition was accomplished with use of an Inversion-Recovery (IR) Hybrid-State Free Precession (HSFP) experiment [11] to sensitize the MRI response signal to $T_{1}$ and $T_{2}$ relaxation. The flip angle pattern is originally optimized for maximal mapping efficiency at $T_{1}=$781\text{,}\mathrm{ms}$,T_{2}=$65\text{,}\mathrm{ms} with $T_{R}=$4.5\text{\,}\mathrm{ms}. Due to specific absorption rate constraints of the slice-selective excitation pulses, we prolonged the repetition time in this study to $T_{R}=$5\text{,}\mathrm{ms} which scales the expected maximal efficiency to be achieved at $T_{1}=$868\text{\,}\mathrm{ms}$,T_{2}=$72\text{\,}\mathrm{ms}. The flip angle pattern was implemented with 2D radial sampling using a tiny-Golden-Angle scheme [19] with \Psi_{N=10}\approx$$$. Phantom and brain imaging was performed with a spatial resolution of 0.75\text{\times}0.75\text{\times}4\text{,}\mathrm{mm}$in a total acquisition time of$T_{\text{ACQ}}=4.3\text{\,}\mathrm{s}$$. Signal time courses and corresponding gradients were computed using the analytical expression for HSFP (compare eq. 7 in Ref [[11](#bib.bib11)]). Slice profile effects were explicitly taken into account by evaluating this expression with different B_{1}$ strengths. To this end, the time course of the implemented RF excitation pulse (bandwidth-time-product 2) was Fourier transformed and one of the two symmetric lobes discretized into 20 factors that scale each flip angle in the excitation pattern. The slice-profile compensated signal time course was eventually obtained by averaging the 20 individual time courses.

Subspace and reconstruction

To further reduce dictionary size and to minimize noise amplification, we formulate the reconstruction as a subspace-constrained [3, 5] linear inverse problem. A subspace size of $K=4$ was chosen heuristically and the subspace basis was determined by performing a SVD [3, 13, 20] on either the full adaptive dictionary (phantom and brain study) or on a set of signals from a uniform grid in $T_{1}$-$T_{2}$ parameter space (numerical simulation). The latter strategy ensures a dictionary-independent basis. With this choice, the following minimization problem is solved:

[TABLE]

where $y$ denotes the radial raw data, $\mathcal{P}_{\vec{k}}$ the projection onto the sampled k-space trajectory, $\mathcal{F}$ the Fourier transform, $S$ multiplication with the (predetermined) coil sensitivity profiles, $\Phi_{K}$ the temporal basis, and $\alpha$ the unknown subspace coefficients. Coil sensitivity profiles $S$ were predetermined by ESPIRIT [21] using the gridding solution of the first subspace coefficient, and spatial correlations across subspace coefficients were exploited by a locally-low-rank regularizer $R$ [20, 22].

Similar to our previous work [22], gridding, gradient delay correction, and precomputation of the transfer point-spread-function is performed by custom MatLab routines, while image reconstruction was performed by a customized version of BART [23] using the ADMM optimizer [24] ($\rho=0.01$, 100 iterations) and locally-low-rank regularization ($\lambda=0.0003$ and block size $8\times 8$).

The linear subspace transformation $\Phi^{\top}_{K}$ is also applied to each entry in the dictionary $D$ and each Jacobian matrix $J$, such that the manifold projection is simply performed with their subspace representations $\hat{D}=\Phi^{\top}_{K}D$ and $\hat{J}=\Phi^{\top}_{K}J$.

Phantom design

For a quantitative validation, a home-brew phantom was designed consisting of 9 gel tubes with distinct $T_{1}$ and $T_{2}$ values. Closely following Ref. [25], $\mathrm{GdCl_{3}}$ was used as a $T_{1}$ modifier and agarose as a $T_{2}$ modifier to generate $T_{1}$ and $T_{2}$ values in the range of typical relaxation times for white and gray matter. To access the power of separability of the proposed method, $T_{1}$ was kept approximately constant while varying $T_{2}$ and vice versa. Ground truth $T_{1}$ relaxometry was realized by four IR single-echo spin-echo data acquisitions [26] (TI=$30\text{\,}\mathrm{ms}530\text{\,}\mathrm{ms}1030\text{\,}\mathrm{ms}1530\text{\,}\mathrm{ms}$) and pixel-wise fitting of the complex data using a freely available custom software package [27]. $T_{2}$ gold standard values were obtained by 5 single-echo spin-echo data acquisitions (TE = $12\text{\,}\mathrm{ms}30\text{\,}\mathrm{ms}73\text{\,}\mathrm{ms}182\text{\,}\mathrm{ms}450\text{\,}\mathrm{ms}$, TR=$4.5\text{\,}\mathrm{s}$) and subsequent fitting of a mono-exponential model to the magnitude data.

$B_{1}$ profile correction

Local deviations in the $B_{1}$ field are known to be a major source of systematic errors in quantitative MRI. This specifically applies to the utilized HSFP sequence as information about $T_{1}$ and $T_{2}$ is encoded in the signal response by traversing the Bloch sphere on a particular path [11]. Imperfect $B_{1}$ field strength leads to deviations from the intended path and mainly results in inaccurate $T_{2}$ values, similar to $B_{1}$ effects in ”MR Fingerprinting” [8, 28]. To correct for this $B_{1}$ deviations, a separate $B_{1}$ map was acquired for the phantom study using a standard sequence of the vendor (Bloch-Siegert method [29]) which matched the spatial resolution of the HSFP sequence. The pixel-wise information about the relative scaling of the nominal $B_{1}$ strength, $rB_{1}$, was used to correct at the stage of the manifold projection. In the spirit of the local-linear approximation, the computed Jacobian matrices can be extended to incorporate the derivative with respect to $rB_{1}$ and Equation 6 is extended to provide signal models for different $B_{1}$ strength.

Code availability

The source code will be made publicly available via https://github.com/ at the time of publication.

Results

Simulation

To investigate the role of the error tolerance $\varepsilon$, adaptive dictionaries were generated for decreasing values of $\varepsilon$ and compared to the heuristic dictionary. With decreasing $\varepsilon$, the number of dictionary entries increases (Figure 1A-E) and the intended adaptivity effect becomes apparent: In each dictionary, the sampling density increases toward the short-$T_{1}$-short-$T_{2}$ region. Note the close similarity between the automatically and the heuristically generated density distributions (Figure 1F).

In Figure 2, these dictionaries have been used to project a probing signal response ($T_{1}=$1.088\text{,}\mathrm{s}, $T_{2}=$0.069\text{\,}\mathrm{s}, $\rho=$1) that is not contained in any of the dictionaries. With decreasing $\varepsilon$, the relative error in $T_{1}$, $T_{2}$, and $\rho$ generally decreases and finally falls below the heuristic error and levels below $1.1\text{\,}\mathrm{\text{\textperthousand}}$. Comparing the dictionary sizes ([Figure 2](#Sx6.F2)B) reveals that the adaptive sampling strategy, depending on the chosen error tolerance, results in dictionaries reduced in size by one to two orders of magnitude compared to heuristic sampling. Based on the excellent accuracy obtained with only 181 dictionary entries, the error tolerance of $\varepsilon=$0.06 was used for both phantom and in vivo studies.

Phantom experiments

Figure 3A shows the four reconstructed subspace coefficient maps obtained for the $T_{1}$-$T_{2}$ phantom using the HSFP flip angle pattern (Figure 3B). The subspace approach allows to store the signal responses for all sampled $T_{1}$-$T_{2}$ parameter combinations (Figure 3C) in a compressed representation with four coefficients per sampling point in the dictionary (basis functions shown in Figure 3D).

The reconstructed subspace coefficients are then mapped pixel-wise to yield the final $T_{1}$, $T_{2}$ and proton density maps (Figure 4). This is done by the proposed manifold projection using the adaptively generated dictionary (Figure 4A), as well as by template matching with the heuristic dictionary (Figure 4B). Quantitative comparison of ROI-wise mean and standard deviations shows excellent agreement between these two methods except for the longest-$T_{1}$-longest-$T_{2}$ tube (upper right grid position). Here, the template matching approach leads to a ”cartoon” artifact in the $T_{2}$ map. The entire compartment is mapped to a constant value of $0.2\text{\,}\mathrm{s}$ resulting in a vanishing standard deviation. Sampling was obviously too coarse in the heuristic dictionary in the region around $T_{2}=$0.2\text{,}\mathrm{s}$$ (see Figure 1E).

The quantitative values are in general agreement with the gold standard measurements, however, for both mapping methods a $T_{2}$ bias is noticeable. Applying the proposed $B_{1}$ profile correction in the manifold projection (Figure 4C) corrects this bias to a large degree.

In Vivo experiments

Figure 5 shows subspace coefficients and parameter maps of a transversal section of the human brain. The parameter maps reveal excellent agreement and demonstrate the efficient use of tiny dictionaries for multi-parametric mapping in vivo.

DISCUSSION

In this work, the locally-linear model for joint $T_{1}$ and $T_{2}$ mapping was built with the analytical Jacobian. Depending on the employed signal model (extended phase graphs (EPG), full Bloch simulation, etc.), this computation can be cumbersome. In this case, a proper replacement for the exact Jacobian matrix is required in the proposed manifold projection. After identification of a specific neighborhood, an approximate Jacobian matrix can easily be obtained by linear regression analysis using the reference position and its neighborhood. With this approximation, the proposed manifold mapping is also applicable to cases in which analytical information on signal derivatives is not easily available.

The quantitative evaluation of the $T_{1}$-$T_{2}$ phantom demonstrated that the proposed $B_{1}$ correction removes the original $T_{2}$ bias to a large degree. However, in particular the short-$T_{2}$ compartments showed a remaining bias which can probably be attributed to neglected effects in the signal model such as finite RF pulses, $T_{2}$-dependent inversion efficiency, and magnetization transfer.

Here, we implemented the HSFP excitation pattern in a 2D (rather than 3D) sequence, so that the signal response becomes a through-slice average. While a proper compensation in the forward model was possible, the excitation pattern is suboptimal in terms of mapping efficiency. Although a rigorous optimization including the slice profile was beyond the scope of this work, it would increase the mapping efficiency.

In contrast to non-linear model-based reconstruction techniques, linear subspace-constrained techniques are inherently tolerant to partial voluming and allow fast reconstruction. However, the choice of the subspace size always becomes a trade-off between noise amplification and model-error. Therefore, it would be highly desirable to combine the following techniques: A non-linear signal model for optimal use of data redundancy, and its embedding in a linear subspace for computational efficiency would eventually make the noise amplification independent of the subspace size. The proposed manifold projection could constitute a key role in such a fused reconstruction technique, but further investigations are necessary.

In conclusion, a novel method to adaptively generate dictionaries for multi-parametric mapping was introduced. The quantitative results for $T_{1}$-$T_{2}$ mapping showed excellent agreement between the proposed manifold projection using adaptive dictionaries and template matching using heuristic dictionaries. The demonstrated ability to perform reconstruction and parameter mapping entirely in subspace justifies the coined term ”tiny dictionaries”. The proposed technique has the potential to overcome problems associated with large dictionaries in quantitative multi-parametric mapping.