Myocardial tissue characterization by T1-T2 mapping is promising for diagnosing myocardial infarction, ischemia, and more. This is typically performed using ECG triggering and breath-holding, which is uncomfortable and unreliable in patients. This work describes a novel method for non-ECG free-breathing joint T1-T2 mapping using cardiovascular MR multitasking, modeling the underlying 6D multidimensional image—which has 2 spatial dimensions + 4 time dimensions (cardiac, respiratory, T1, and T2)—as a low-rank tensor. T1 and T2 measurements in acute myocardial infarction patients agreed with reference methods and predicted late gadolinium enhancement with 100% sensitivity and 92% specificity.
T1-T2 weighting was generated using T2prep/IR preparation. T1 recovery between prep pulses was sampled using FLASH readouts (flip angle=5°, TR/TE=3.6/1.6 ms) throughout the 2.5-sec recovery period. T2prep duration was cycled through 12, 20, 30, 40, and 50 ms to provide multiple T2 weightings (Fig. 1). Sampling was performed by interleaving golden-angle radial data (imaging data $$$\mathcal{D}_\mathrm{img}$$$) and the 0° spoke (LRT subspace training data $$$\mathcal{D}_\mathrm{tr}$$$) every other readout.
$$$\quad$$$We represent the underlying multidimensional cardiovascular image as a 5-way tensor $$$\mathcal{A}$$$, with dimensions indexing voxel location $$$\mathbf{x}$$$, cardiac phase $$$t_\mathrm{c}$$$, respiratory position $$$t_\mathrm{r}$$$, time since preparation $$$t_\mathrm{T1}$$$ (representing T1 recovery), and T2prep duration $$$t_\mathrm{T2}$$$ (representing T2 decay). Spatiotemporal image correlation makes this tensor low-rank3. In Tucker form, this is expressed as:$$\mathcal{A=G}\times_1\mathbf{U_x}\times_2\mathbf{U}_\mathrm{c}\times_3\mathbf{U}_\mathrm{r}\times_4\mathbf{U}_\mathrm{T1}\times_5\mathbf{U}_\mathrm{T2},$$where $$$\mathcal{G}$$$ is the core tensor, and each factor matrix $$$\mathbf{U}$$$ contains basis functions describing each dimension. We recover $$$\mathcal{A}$$$ by learning $$$\mathcal{G}$$$ and the temporal factor matrices from $$$\mathcal{D}_\mathrm{tr}$$$, followed by recovery of $$$\mathbf{U_x}$$$ from $$$\mathcal{D}_\mathrm{img}$$$ using the other known factors.
$$$\quad$$$Although the training data in $$$\mathcal{D}_\mathrm{tr}$$$ are frequently collected, they do not fully cover $$$(t_\mathrm{c},t_\mathrm{r},t_\mathrm{T1},t_\mathrm{T2})$$$-space. Therefore, before extracting $$$\mathcal{G}$$$ and the temporal factor matrices, we complete $$$\mathcal{D}_\mathrm{tr}$$$ according to:$$\hat{\mathcal{D}}_\mathrm{tr}=\arg\min_{\mathcal{D}_\mathrm{tr}\in\Psi}\|\mathbf{d}_\mathrm{tr}-\Omega_\mathrm{tr}(\mathcal{D}_\mathrm{tr})\|_2^2+\lambda\|\mathbf{D}_{\mathrm{tr},(n)}\|_*+R(\mathcal{D}_\mathrm{tr}),$$where $$$\Psi$$$ is a subspace estimated from the SVD of a dictionary of T1 recovery curves2, where $$$\Omega_\mathrm{tr}$$$ returns the vector of sampled training data $$$mathbf{d}_\mathrm{tr}$$$, and where subscript $$$(n)$$$ denotes mode-$$$n$$$ unfolding of the tensor. $$$R(\cdot)$$$ is a regularization functional, which for this work was chosen as the temporal total variation (TV) along $$$t_\mathrm{c}$$$ and $$$t_\mathrm{r}$$$.4 $$$\mathcal{G}$$$ and the temporal factor matrices are then extracted from the high-order SVD of $$$\hat{\mathcal{D}}_\mathrm{tr}$$$.
$$$\quad$$$We then recover $$$\mathbf{U_x}$$$ according to:$$\mathbf{\hat{U}_x}=\arg\min_\mathbf{U_x}\|\mathbf{d}_\mathrm{img}-E(\mathbf{U_x\Phi})\|_2^2+R(\mathbf{U_x}),$$where $$$E(\cdot)$$$ combines spatial encoding and sampling, $$$R(\cdot)$$$ is a regularization functional (in this work, spatial TV), and where $$$\mathbf{\Phi}=\mathbf{G}_{(1)}(\mathbf{U}_\mathrm{T2}\otimes\mathbf{U}_\mathrm{T1}\otimes\mathbf{U}_\mathrm{r}\otimes\mathbf{U}_\mathrm{c})^T$$$.
$$$\quad$$$Data were collected on a 3T Siemens Verio from n=5 patients with AMI. Each subject underwent: 1) native T1 mapping using MOLLI 5(3)36; 2) native T2 mapping using T2prep-SSFP mapping7; 3) native T1-T2 mapping using multitasking; and 4) LGE imaging. MOLLI, T2prep-SSFP, and LGE images were collected at diastole during an end-inspiration breath-hold. CMR multitasking images were continuously acquired without ECG triggering or breath-holding over the course of 85 sec, resulting in 6D cardiac- and respiratory-resolved images (16 cardiac bins, 5 respiratory bins, 344 inversion times, 5 T2prep durations); T1 and T2 values were fit at end-diastole and end-inspiration to match the reference method motion states. Segmentwise analysis was performed in six myocardial segments per subject.
Representative mapping results are in Fig. 2. Comparison of multitasking T1 and T2 values (T1: 1177±147 ms; T2: 45.6±8.3 ms) to MOLLI (T1: 1206±31 ms) and T2prep-bSSFP mapping (T2: 44.3±6.7 ms) were positively correlated (T1: r=0.90, T2: r=0.68) (Figs. 3a,4a). Multitasking T1 showed a –29 ms bias (p=0.02) vs. MOLLI; T2 had no statistically significant bias (p=0.30) (Figs. 3b,4b).
$$$\quad$$$Neither T1 mapping method on its own was predictive of LGE status, with the area under the curve (AUC) of the receiver operating characteristic (ROC) curve = 0.55 for MOLLI and 0.50 for multitasking (Fig. 3c,d); both T2 mapping methods were predictive of LGE status, with the ROC AUC = 0.71 for T2prep-bSSFP mapping and 0.84 for multitasking (Fig. 4c,d). A 2D feature space of multitasking T1 and T2 values (Fig. 5) provided even greater diagnostic accuracy, with the pictured decision boundary yielding 100% sensitivity and 92% specificity.