### Data acquisition and model-based reconstruction

The single-shot IR scheme used here has been reported before [11]. For myocardial T1 mapping, data acquisition starts with a non-selective inversion pulse which is triggered to the early diastolic phase with use of a finger pulse signal. After inversion, the signal is continuously acquired for a period of 4 s using a radial FLASH readout with a golden-angle trajectory. To eliminate motion effects during systolic contraction and expansion, only data from the diastolic phase is retrospectively selected for T1 mapping.

The signal from multiple coils is given by

$$ {\mathrm{y}}_j(t)=\int {M}_{t_k}\left(\overrightarrow{r}\right){c}_j\left(\overrightarrow{r}\right){e}^{-i\overrightarrow{r}\overrightarrow{k}(t)}d\overrightarrow{r} $$

(1)

with *c*_{j} the jth coil sensitivity map, \( \overrightarrow{k}(t) \) the chosen k-space trajectory, y_{j}(*t*) the acquired data and \( {M}_{t_k}\left(\overrightarrow{r}\right) \) the magnetization at time *t*_{k} after inversion

$$ {M}_{t_k}={M}_{ss}-\left({M}_{ss}+{M}_0\right)\cdotp {e}^{-{t}_k\cdotp {R}_1^{\ast }} $$

(2)

where *t*_{k} is defined as center of the acquisition window in this study. \( {M}_{ss},{M}_0\ \mathrm{and}\kern0.5em {R}_1^{\ast } \) represent the steady-state signal, equilibrium signal and effective relaxation rate, respectively. After estimation of \( \left({M}_{ss},{M}_0,{R}_1^{\ast}\right) \), T1 can be calculated by

$$ \mathrm{T}1=\frac{M_0}{M_{ss}\cdot {R}_1^{\ast }}\kern0.5em $$

(3)

In Eqs. (1) and (2), both the model parameters \( {\left({M}_{ss},{M}_0,{R}_1^{\ast}\right)}^T\ \mathrm{and}\ \mathrm{all}\ \mathrm{coil}\ \mathrm{sensitivity}\ \mathrm{maps}\ {\left({c}_1,\cdots, {c}_N\right)}^T \) are unknowns, which are directly estimated from k-space using a sparsity constrained model-based reconstruction, i.e.,

$$ \hat{x}=\mathrm{argmin}{\left\Vert F(x)-y\right\Vert}_2^2+\alpha R\left({x}_{\boldsymbol{p}}\right)+\beta Q\left({x}_{\boldsymbol{c}}\right) $$

(4)

Here *F* is the nonlinear forward model mapping all unknowns to the measured data *y*:

$$ F:x\mapsto \left(\begin{array}{c}{P}_1\mathcal{F}\left\{{c}_1\cdot {M}_{t_1}\left({M}_{ss},{M}_0,{R}_1^{\ast}\right)\right\}\\ {}\vdots \\ {}{P}_1\mathcal{F}\left\{{c}_N\cdot {M}_{t_1}\left({M}_{ss},{M}_0,{R}_1^{\ast}\right)\right\}\\ {}{P}_2\mathcal{F}\left\{{c}_1\cdot {M}_{t_2}\left({M}_{ss},{M}_0,{R}_1^{\ast}\right)\right\}\\ {}\vdots \\ {}{P}_n\mathcal{F}\left\{{c}_N\cdot {M}_{t_n}\left({M}_{ss},{M}_0,{R}_1^{\ast}\right)\right\}\end{array}\right) $$

(5)

with *P* the orthogonal projection onto the trajectory and \( \mathcal{F} \) the 2D Fourier transform. The unknowns \( {x}_{\boldsymbol{p}}={\left({M}_{ss},{M}_0,{R}_1^{\ast}\right)}^T \) and *x*_{c} = (*c*_{1}, ⋯, *c*_{N})^{T}. *R*(*x*_{p}) is a L1-Wavelet regularization which exploits joint sparsity in the parameter dimension following the ideas of compressed sensing, while *Q*(*x*_{c}) is a Sobolev norm which is applied to the coil sensitivities to enforce their intrinsic smoothness. α and β are the corresponding regularization parameters. The nonlinear inverse problem in Eq. (4) is solved by the iteratively regularized Gauss-Newton method (IRGNM) [24] where the nonlinear problem is linearized in each Gauss-Newton step and solved by the fast iterative shrinkage-thresholding algorithm (FISTA) [25]. More details of the IRGNM-FISTA algorithm can be found in [20].

### CMR

All CMR studies were conducted on a 3 T system (Magnetom Skyra, Siemens Healthineers, Erlangen, Germany) with approval of the local ethics committee. Phantom measurements employed a 20-channel head/neck coil, while human heart studies used a combined thorax and spine coil with 26 channels. Eight subjects (three female, five male, age 27 ± 3, range 23–32 years; heart rates 62 ± 11 bpm, range 50–80 bpm) with no known illness were recruited. Written informed consent was obtained from all subjects prior to CMR. In vivo T1 measurements were performed within a single breathhold.

The proposed method was experimentally validated at simulated heart rates with a commercial reference phantom (Diagnostic Sonar LTD, Livingston, Scotland, UK) consisting of six compartments with defined T1 values surrounded by water. The gold standard T1 map for the phantom was estimated using an IR spin-echo method [26] with 9 IR scans (TI = 30, 530, 1030, 1530, 2030, 2530, 3030, 3530, 4030 ms), TR/TE = 4050/12 ms, FOV 192 × 192 mm^{2}, matrix size 192 × 192, and a total acquisition time of 2.4 h.

For IR radial FLASH, continuous data acquisition was performed with a tiny golden angle (18.71°) [27] after non-selective inversion. Because there is no intermediate image reconstruction, model-based reconstructions offer a flexible choice of temporal resolution, i.e., they allow a combination of an arbitrary (small) number of radial spokes for each k-space frame. However, as long as the T1 accuracy is not compromised, a certain degree of temporal discretization (data binning) is recommended to reduce the computational demand [19, 20]. In this study, 17 spokes formed one k-space and resulted in a temporal resolution of 45 ms. According to the subjects’ heart rates, the resulting number of k-space frames were 48 ± 9, range 33–57 for reconstructions in this study. Single-shot myocardial T1 maps of the mid-ventricular slices were acquired at a nominal in-plane resolution of 1.0 × 1.0 mm^{2} and 8 mm slice thickness using a FOV 256 × 256 mm^{2} in combination with a resolution of 512 complex data points per radial spoke (two-fold oversampling). Other parameters were TR/TE = 2.67/1.67 ms, nominal flip angle 6°, bandwidth 850 Hz/pixel and total acquisition time 4 s.

To access reproducibility of the proposed method, the single-shot sequence was performed 3 times on each subject: The first two measurements were repeated one after the other, while the third one was done with a 5-min break, during which time the subject was taken out of the scanner. For comparisons, single-shot T1 maps were also estimated using the frame-based nonlinear inversion (NLINV) reconstruction with subsequent pixel-wise fitting as described in [11] without and with spatial filtering by a modified nonlocal means filter [28] from the same datasets. Further, a 5(3)3 MOLLI sequence provided by the vendor was applied for reference using a FOV of 360 × 306.6 mm^{2}, in-plane resolution 1.41 × 1.41 × 8 mm^{3}, TR/TE = 2.24/1.12 ms, nominal flip angle 35°, bandwidth 1085 Hz/pixel and total acquisition time 11 heart beats.

### Implementation

All data was processed off-line. Multicoil raw data were first corrected for gradient delays [29] and then compressed to 10 virtual channels using a principal component analysis (PCA). A convolution-based gridding [30] without density compensation was used to interpolate the radial samples onto a Cartesian grid on which all successive iterations were performed. All the computations were done in Berkeley advanced reconstruction toolbox (BART) [31] on a 40-core 2.3 GHz Intel Xeon E5–2650 PC with a RAM size of 500 GB.

The parameter maps \( {\left({M}_{ss},{M}_0,{R}_1^{\ast}\right)}^T\ \mathrm{were}\ \mathrm{initialized}\ \mathrm{with}\ {\left(1.0,1.0,1.5\right)}^T \) and all coil sensitivities zeros for all reconstructions. 10 Gauss-Newton steps were employed to ensure convergence. Similar to [20], regularization parameters α and β were initially set to 1 and subsequently reduced by a factor of 3 in each Gauss–Newton step. A minimum value of α was used to control the noise at higher Gauss–Newton steps. The chosen value of *α*_{min} was defined by optimizing signal to noise ratio (SNR) without compromising quantitative accuracy or delineation of structural details. With the above settings, the whole computation took around 6 h using the CPUs. However, with a reduced number (e.g., 6) of virtual coils, computations could be run on a GPU, which took 10 to 20 min per dataset

### Data analysis

Results in this work are reported as mean ± standard deviation (SD). For the assessment of myocardial T1 values, the regions of interest (ROIs) in the inter-ventricular septum were carefully selected to exclude the blood pool using arrShow [32] tool in MATLAB (MathWorks, Natick, Massachusetts, USA) and performed by two independent observers. Similar to [8, 33], the precision of T1 estimation was evaluated using coefficient of variation (CV = SD_{ROI}/Mean_{ROI} × 100%). The reproducibility error was calculated by \( \sqrt{\left({\sum}_{i=1}^{n_s}\mathrm{T}{1}_{\mathrm{diff}}^2(i)\right)/{n}_s}, \) where T1_{diff}(*i*) is the T1 difference between different measurements, *n*_{s} is the number of subjects. Further, a repeated measures analysis of variance (ANOVA) with Bonferroni post hoc test was used for comparisons and a *P* value < 0.05 was considered significant.

In addition, edge sharpness was quantitatively measured for both the proposed model-based reconstruction and MOLLI. It was done by fitting each septal T1 line profile (starting from the blood pool to the middle of the myocardial septum) to a parameterized sigmoid function [34]: \( s\left(\mathrm{x}\right)=\frac{\mathrm{a}}{1+{\mathrm{e}}^{-\mathrm{k}\cdot \left(\mathrm{b}-\mathrm{x}\right)}}+c \), where x is the length (unit: millimeter) along the line profile and (a, b, c, k)^{T} are the fitting parameters: a determines the vertical range, b determines the center location, c defines the vertical offset and k quantifies the growth rate or sharpness of the edges (The higher |k|, the sharper the edges). The above nonlinear least square fitting was then performed in MATLAB (MathWorks) using the Levenberg-Marquardt algorithm with a stopping criteria similar to [11].