Autocalibrated cardiac tissue phase mapping with multiband imaging and k‐t acceleration

To develop an autocalibrated multiband (MB) CAIPIRINHA acquisition scheme with in‐plane k‐t acceleration enabling multislice three‐directional tissue phase mapping in one breath‐hold.


| INTRODUCTION
Tissue phase mapping (TPM) is a phase-contrast (PC) MRI technique for the velocity measurement of myocardial wall motion, 1,2 and constitutes a valuable tool for the diagnosis of abnormal ventricular function in a range of cardiovascular medical conditions. [3][4][5][6][7] Tissue phase-mapping protocols commonly combine 2D time-resolved gradient-echo PC-MRI sequences with one, two or three-directional velocity encoding. Consequently, scan times are substantially increased with respect to 2D gradient-echo cine sequences. In particular, three-directional velocity encoding extends scan times beyond feasible breathhold duration. A common strategy consists of using navigator scans to track the diaphragm motion during free breathing. However, even when navigators are being used, the image quality is worse due to respiratory motion. 8 Furthermore, a decrease in peak velocities was observed in TPM data acquired during free-breathing versus breath-hold. 9 Therefore, it is desirable to accelerate TPM acquisitions so that they can be conducted within a single breath-hold.
In addition to using standard parallel imaging techniques such as GRAPPA 10 and SENSE 11 for in-plane acceleration of TPM, 12 undersampling in both the phase encoding and the temporal domain (k-t acceleration 13 ) has been used. 14 A concern when applying k-t acceleration for velocity quantification consists of preserving the velocity time courses when high acceleration factors are being used. 15 Although the time required to acquire a single slice can be reduced using k-t acceleration methods, the acquisition of multiple slices, a strategy that is typically performed to cover basal, midventricular and apical motion of the myocardium, still requires separate breath-holds. Such TPM protocols increase the overall scan time and are susceptible to the patient's potential incapability to reproduce breath-holds, which in turn may lead to spatial mismatch of the slices. This could be avoided by acquiring several slices at the same time within a single breath-hold.
Simultaneous multislice imaging, 16 also called multiband (MB), in combination with CAIPIRINHA encoding, 17,18 is a promising technique to address this problem, as it allows to excite and acquire multiple slices at the same time. The benefits of MB have been demonstrated in cardiac MRI, [19][20][21][22][23][24][25][26][27] and the technique was recently extended to TPM with one-directional velocity encoding. 28 Two recent studies have also demonstrated the benefit of autocalibrated MB, 28,29 a variant of MB-CAIPIRINHA encoding that does not require external reference scans for the reconstruction, thus limiting the number of required breath-holds to only one.
In this study, we seek to address the unmet need of a TPM acquisition scheme capable of performing three-directional velocity encoding of the myocardial motion velocities at multiple slices in a single breath-hold. To achieve this, autocalibrated MB and in-plane k-t acceleration were combined. A k-t SENSE reconstruction algorithm that simultaneously tackles k-t and multiband acceleration was developed and optimized such that the temporal fidelity of the reconstructed velocity times courses is maximized. By using the proposed methods, an up to 11-fold net acceleration was achieved, allowing three-directional velocity encoding of the myocardium at two/three slices with a temporal resolution of 38.4 ms and inplane spatial resolution of 2.3 mm, to be acquired in a single breath-hold of 18 heartbeats.

| Sequence description
A prospectively electrocardiographic-triggered cardiac PC 2D cine MR sequence with three-directional velocity encoding 30 was extended to include k-t acceleration 13 and autocalibrated MB CAIPIRINHA. 28,29 This combination is not straightforward, because for each technique certain conditions with respect to the sampling patterns and the RF phases have to be fulfilled. Before presenting the rationale for the combination of the two techniques, a brief introduction to autocalibrated MB CAIPIRINHA and k-t acceleration is provided.
In autocalibrated MB CAIPIRINHA, the simultaneously excited slices appear shifted relative to one another through phase modulation of the RF-MB excitation pulses. 17 The arrangement of these pulses along the phase-encoding and temporal (cardiac phase) dimension is depicted in Figure 1. For the case of MB2 ( Figure 1A), pulse A is calculated by the | 2431 FERRAZZI, BASSENGE, Et Al. slice 1, and of π for slice 2. The pair of RF phases applied to the two slices is referred to as (0, 0) for pulse A and (0, π) for pulse B. Similarly, in MB3 ( Figure 1B), pulses A, B, and C are calculated by the summation of three SB pulses with RF phases (0, 0, 0), 0, 2π 3 , 4π 3 and 0, 4π 3 , 2π 3 , respectively. Because the transmit phases of the RF pulses are linearly increased as k-space lines are acquired, CAIPIRINHA encoding imparts slice-dependent shifts along the phase-encoding direction, which in turn favors parallel imaging reconstruction. 17 Using the encoding structure of Figure 1, the shift applied to slice 2 in the case of MB2 is equal to half of the FOV FOV 2 , whereas the first slice (slice 1) is not spatially shifted. In the case of MB3, slices 2 and 3 are shifted by FOV 3 and − FOV 3 relative to slice 1, which is located in the center of the FOV.
To extract the coil sensitivities from the actual MB acquisition, the temporal domain is exploited in autocalibrated MB. A key condition of this technique is that the RF phases of the MB pulses are cycled in time according to the schemes of Figure 1A,B along the horizontal axis. In doing so, consecutive cardiac phases can be grouped to generate slice unaliased coil-calibration data, 28,29 which can subsequently be used to unfold the slices.
In this work, a k-t scheme with 5-fold nominal acceleration (effective in-plane acceleration R k-t = 3.67) was combined with autocalibrated multiband factors R MB of 2 and 3. In Figure 1, colored versus lightly colored locations in k-t space correspond to acquired versus not acquired data points as the result of k-t undersampling. Note that the undersampling schemes for MB2 and MB3 are different. These were optimized independently, similarly to what has been proposed in Tsao et al. 31 The central k-space lines instead correspond to a fully sampled region.
To ensure that the steady-state condition across consecutive heartbeats is fulfilled, 28 the acquisition order of the k-space lines is chosen such that only identical RF pulses are transmitted within each heartbeat interval. As a result, the signal evolution across consecutive heartbeats is the same, as the magnetization vector at every slice location is subject to the same excitation history (apart from a global phase offset to achieve CAIPIRINHA encoding, which does not provoke F I G U R E 1 Radiofrequency pulses plotted as a function of cardiac phases and k-space lines for the proposed k-t acquisition schemes. In this representation, two pulses (A and B) are generated for multiband 2 (MB2) (A, left), and three (A, B, and C) for MB3 (B, right). A, Case MB2. Pulse A is calculated by the summation of two single-band (SB) RF pulses with RF phases set to 0 for the excitation of both slices (1 and 2), whereas pulse B is calculated by the summation of two SB RF pulses with RF phase of 0 for slice 1, and of π for slice 2. Colored versus lightcolored locations represent acquired versus not acquired data points due to k-t acceleration. In this example, a fully sampled autocalibrated region consisting of six lines plus one side of the outer, subsampled part of k-space is shown. The order in which the k-space lines are acquired as a function of consecutive heartbeats is depicted in Supporting Information Video S1. B, Same as (A) for the case of MB3. In this case, A, B, and C are formed by the summation of SB pulses with RF phases (0, 0, 0), 0, 2π a difference in signal behavior across slices or heartbeats). In this study, RF-phase sequences of type ABAB… for MB2 and ACBACB… for MB3 were selected (see Supporting Information Video S1). Note, however, that other designs are possible.

| Image acquisition
The study was approved by the local ethics board, and written informed consent was obtained from all subjects. Acquisitions were performed at 3 T (Magnetom Verio; Siemens, Erlangen, Germany). A vendor-provided electrocardiographic device was used to synchronize the acquisitions to the cardiac cycle. A 16-channel anterior cardiac coil combined with a 12-channel posterior spine coil was used to acquire the data.
Data were acquired in short-axis orientation with readout direction right-left and phase-encoding direction anteriorposterior in 4 subjects (S1 to S4) using FOV = 400 × 300 mm 2 , slice thickness = 6 mm, interslice distance = 20 mm, in-plane resolution = 2.3 mm, TR = 4.8 ms, TE = 3.1 ms, flip angle = 10°, and a readout bandwidth of 810 Hz/ pix with asymmetrical echo (20%). All data were acquired using three-directional velocity encoding (VENC z = 50 cm/s through plane, VENC xy = 30 cm/s in plane). The temporal resolution was 38.4 ms, and the number of acquired k-space lines per heartbeat was 2. Twelve fully sampled lines were acquired at the center of k-space, and the total acquisition used 18 heartbeats. The effective acceleration R 0 = R k-t * R MB was 7.34 and 11.01 for MB2 and MB3 acquisitions, respectively. In each scanning session, the number of cardiac phases was adjusted individually to match the subject's heartbeat duration. For each MB data set, SB k-t accelerated reference data with identical parameters to the MB k-t scans were acquired during separate breath-holds at equivalent slice locations. The number of SB k-t acquisitions matched the MB factor (i.e. two acquisitions for the MB2 k-t data, and three for MB3 k-t).
In a fifth scan (subject S5), fully sampled (FS) SB slices with one-directional velocity encoding (v z only) were acquired as ground truth to investigate temporal fidelity issues (see section 2.5). The short-axis SB slices were acquired at six different locations as part of six separate breath-holds: from apex (data set SL1) toward base (data set SL6) in steps of 1.4 cm. Acquisition parameters were closely matched to those of the actual MB k-t accelerated data described previously, including the same spatial and temporal resolutions. To match the temporal resolution, the receiver bandwidth was adjusted from 810 to 660 Hz/pix, and the number of k-space lines acquired in each heartbeat interval was set to four (from two).

| Reconstruction algorithm
Data were reconstructed using a hybrid space 32 temporally hierarchical k-t SENSE 13 formulation that is solved by conjugate gradient as follows: Here, j indexes the hierarchy of reconstructions; x is the spatio/temporal signal to be reconstructed; S represents the static coil sensitivities calculated as described in section 2.4; and A j is the encoding operator, thus comprising the MB RF phases, the sampling encoding k-t lattice shown in Figure  1A,B, and the Fourier transform, obtained at level j by grouping together the samples of N j cardiac phases. In y j , the measured data are arranged according to the problem at temporal resolution of level j; λ is a regularization weight; R j is the net acceleration factor at level j; and w j is an estimate of the covariance of the signal in the temporal Fourier domain, which is arranged diagonally into W j to weight the k-t regularization at the next iteration. Finally, F t is the temporal Fourier transform and H σ is a Hann-shaped spatial filter of scale σ. Thus, the method iteratively extracts finely (time-resolved) calibrated information to aid in conditioning the increasingly accelerated reconstructions (starting from the identity matrix W 0 = I d ). The hierarchical plan uses N j = {5, 4, 3, 2, 1, 1}.
In this work, λ and σ were optimized based on forward simulation of the imaging process using ground-truth velocities (see section 2.5). Finally, after MB k-t reconstruction was performed, asymmetrical echo was handled using a projection onto convex sets (POCS) reconstruction algorithm. 33,34

| Coil calibration
The static sensitivities of the receive coils (term S of Equation 1) were estimated directly from the SB k-t and MB k-t accelerated data sets. In the first case, a synthetic fully sampled non-time-resolved k-space data set at full spatial resolution was obtained from the SB k-t accelerated acquisitions by averaging the acquired phase-encode lines from all cardiac phases plus subsequent temporal normalization. For the MB k-t acquisitions, and taking as an example the case MB2, this procedure was repeated twice. The first time, only the odd cardiac phases of Figure 1A were considered, and the second time the even cardiac phases. This generated two fully sampled synthetic MB2 data sets (C 1 and C 2 ) at the full resolution scale, and because of the temporal alternation of the MB pulses in Figure 1A, these images have predefined phase offsets in image space. 28 Finally, this predicted phase behavior can be exploited to generate sensitivity information at both slice locations (S 1 and S 2 ) by applying Hadamard decoding 35 (ie, S 1 = C 1 + C 2 and S 2 = C 1 − C 2 ). The same rationale was extended to generate static sensitivities for the case of MB3.

| Simulation
A simulation was performed to assess the temporal fidelity of the reconstructed velocity time courses and slice-separation performances, and to find suitable parameters for the reconstruction framework. The forward simulation of the imaging process used the FS one-directional velocity-encoded SB slices acquired in subject S5 (see section 2.2). The FS SB slices SL1 to SL6 were retrospectively undersampled in the k-t domain to form synthetic SB k-t data sets and combined in groups of two or three to create artificial MB2/MB3 k-t data sets according to the forward model of Equation (1). The synthetic SB k-t and MB2/MB3 k-t data sets were then reconstructed using the proposed reconstruction framework. Note that to investigate temporally smoothing effects solely deriving from the proposed MB k-t reconstruction scheme, the POCS algorithm was not used in the case of the simulated data. For the case of MB2 k-t, a total of nine slice combinations with slice gaps of 1.4 cm (SL1SL2 to SL5SL6) and 2.8 cm (SL1SL3 to SL4SL6) were simulated. For MB3 k-t, six combinations were assessed (SL1SL2SL3 to SL4SL5SL6 with slice gap of 1.4 cm, and SL1SL3SL5 to SL2SL4SL6 with slice gap of 2.8 cm).
To evaluate temporal fidelity, the RMS error (RMSE) and the average velocity peak reduction (peak red.) at the diastolic peak (cardiac phase 13) of the slice-averaged velocity time courses was computed as introduced in section 2.6. This was performed for artificial SB k-t and MB2/MB3 k-t data sets, each in comparison to FS ground-truth velocities. The error metrics were then averaged across all combinations of slices, and the procedure was repeated for a wide range of the regularization parameters and σ.
Once suitable values for and σ were identified, relative error maps (in %) (which comprise the MB slice leakage and k-t acceleration residual folding artifacts plus noise amplification) were evaluated by computing the absolute value of the complex error between FS data sets against SB k-t and MB2/ MB3 k-t reconstructions normalized by the complex FS data.

| Velocity analysis
Velocities v x , v y , and v z were calculated as the phase difference of the velocity-encoded data. Then, for both MB and SB PC data, the myocardial tissue motion of the left ventricle was analyzed using a custom MATLAB (MathWorks, Natick, MA) tool, 36 which allowed eddy current correction 37 and the manual segmentation of the myocardium. To provide a more representative view of the myocardial tissue motion, the inplane velocities v x and v y were transformed into a polar coordinate system as previously described. 30 In particular, the in-plane motion was characterized by radial (v r ) and tangential (v ) components with origin in the center of mass of the myocardial left ventricle. Conventions are (1) positive longitudinal velocities (v z ) for movements toward the apex, (2) positive radial velocities for myocardial contraction, and (3) positive tangential velocities to characterize clockwise rotations of the myocardium as seen from the apex. Pixelwise velocities were averaged within myocardial masks to obtain slice-averaged velocity time courses. Using these data from all slices and subjects, Bland-Altman plots were constructed for SB versus MB data.
For the statistical analysis, Lilliefors tests were performed to assess whether the velocity data were normally distributed. When this criterion was fulfilled (P ≥ .05), paired t-tests were performed between SB and MB data. Otherwise (P < .05), nonparametric Wilcoxon signed-rank tests were used.
Finally, the American Heart Association (AHA) 12 (for MB2, consisting of two slices) and 16 (for MB3, consisting of three slices) segment model was used to measure the velocities within independent subregions of the heart as described in Ruth et al. 36 Identical statistical tests to those described for the myocardial mask results were performed.

| Simulation
The simulation results are summarized in Figure 2. The RMSE (top) and velocity peak red. (bottom) for SB k-t (left), MB2 k-t (center), and MB3 k-t data sets (right) are reported as a function of λ and σ. Scores are calculated in comparison to FS ground-truth velocities averaged over all simulated combinations of slices. Although RMSE and velocity peak red. of the simulated SB k-t data sets remain, on average, fairly constant for a wide range of regularization parameters (mean ± twice the SD 0.06 ± 0.02 cm/s for RMSE and 0.77% ± 0.62% for velocity peak red. across all tested parameter combinations), our results indicate that there is a greater dependence for the cases of MB2/MB3 k-t. However, in such cases, a value of λ = 7 and σ = 10 produced average RMSE and peak red. ± twice the SD across all combinations of slices (and not across the parameter space λ and σ) of 0.08 ± 0.22 cm/s and 1.03% ± 6.47% for MB2, and 0.13 ± 0.58 cm/s and 2.21% ± 13.45% for MB3 (note that, for the peak red., a positive value denotes a lower velocity value with respect to ground truth). Figures 3A and 4A show the reconstruction results from the simulated data using λ = 7 and σ = 10 for one MB2 and one MB3 slice combinations, respectively. In each panel, magnitude (left) and longitudinal velocities (right) data of FS (top), SB k-t (middle), and MB2/MB3 k-t data (bottom) are shown at peak diastole. From the results of the simulation, it is evident how higher acceleration levels lead to noisier reconstructions. However, with such regularization, the longitudinal velocities averaged over the myocardium are reproduced well across all imaging modalities ( Figures 3B and 4B). The M mode visualization (space vs cardiac phases) of these reconstructions (magnitude and velocity maps) along the green lines in Figures 3A and 4A is depicted in Figure 5 for MB2 (left) and MB3 (right), respectively. Reconstructions at higher regularization leading to reduced noise amplification but smoother velocity time courses and reduced peak velocities are possible. For example, Supporting Information Figure S1 shows the same data as in Figure 4A (case MB3 k-t) reconstructed using λ = 23 and σ = 20.
Because the purpose of this study was to reliably reproduce velocity time courses, λ and σ were set to 7 and 10 to reconstruct the actual MB2 and MB3 k-t accelerated data. Figure 6 reports the magnitude maps of the complex relative error of the data from Figures 3 and 4 of the FS slices against SB k-t and MB2/MB3 k-t data sets. Maps have a "noise-like" structure without apparent residual folding artifacts deriving from k-t and MB acceleration. However, it is characterized by higher values over the myocardium than in the blood due to the lower signal intensities. Average relative error values within the myocardium ± twice the SD across all data sets were 16.59% ± 7.68% for SB k-t, 18.13% ± 12.17% for MB2 k-t, and 19.66% ± 15.59% for MB3 k-t reconstructions. Figure 7 reports the reconstructions of the SB k-t (top) and MB2 k-t data (bottom) from one exemplary subject (S2) at peak diastole. The first column reports magnitude data and, from left to right, the velocity maps v x , v y , and v z before eddycurrent correction is applied. Figure 8 reports the equivalent case for MB3 (subject S3). As mentioned in section 2.3, the POCS algorithm was used to perform all reconstructions of the prospectively accelerated scans. However, for a better comparison with the data used in the simulation, the reconstruction before the POCS algorithm is applied for the case MB3 k-t (Figure 8) is reported as Supporting Information Figure S2. Finally, M mode visualization of these reconstructions along the green lines in Figures 7 and 8 is depicted in Supporting Information Figures S3 and S4 for MB2 and MB3, respectively. Figure 9A shows slice-averaged longitudinal, radial, and tangential velocities from 1 exemplary subject (S1) acquired with MB2. Velocity time courses displayed in blue and yellow were obtained with SB k-t and MB2 k-t acquisitions. In both cases, velocities are consistent with the literature 30,38 ; in particular, the longitudinal velocities v z show shortening toward the apex in systole and lengthening in diastole. As expected, the radial velocities v r exhibit a rapid contraction of the myocardium during systole, followed by dilation during diastole. A more complex behavior is seen in the tangential velocities v , where changes in the rotational direction are observable. Figure 9B shows Bland-Altman plots of the SB versus the MB slice-averaged longitudinal, radial, and tangential velocity values of all cardiac phases, slices, and subjects that were scanned for the case of MB2. Means ± twice the SDs of the differences Δv z , Δv r , and Δv were quantified as 0.13 ± 1.98 cm/s, −0.06 ± 0.88 cm/s, and −0.05 ± 0.99 cm/s. For the case of MB3, a similar behavior of the velocity vectors is observable in subject S4 ( Figure 10A). Figure 10B shows Bland-Altman plots for the case of MB3 across all subjects and slices. Means ± twice the SDs of Δv z , Δv r , and Δv were quantified as 0.31 ± 3.21, 0.1 ± 1.06, and −0.06 ± 1.27 cm/s. After confirming that Δv z and Δv for the case of MB3 were not normally distributed (while all other velocity error scores were), statistical analyses showed that differences between the SB and MB averaged velocities of the myocardial tissue motion at every cardiac phase were not significant (P-values in Figures 9B and  10B), despite the fact that MB and SB data were acquired during separate breath-holds. However, when looking at individual subregions of the AHA model, 3 of all possible 36 cross comparisons (AHA subregions times three velocity directions) had statistically significant differences in the F I G U R E 3 Simulation results for MB2. A, MB2 reconstructions.

| Left ventricular motion analysis
Top to bottom: FS, SB k-t, and MB2 k-t reconstructions from the synthetic experiment shown at peak diastole (cardiac phase 13). In this case, the synthetic data set was formed by fully sampled SB data sets SL1 and SL3 (SL1SL3) from subject S5. Magnitude (left) and longitudinal velocities v z (right) are reported. Green lines highlight the positions along which data are displayed in M visualization mode (space vs time) in Figure 5. All data are shown before POCS reconstruction is applied. B, Velocity traces averaged within the myocardium from Figure 3A as a function of the cardiac cycle. Light blue denotes FS; blue denotes SB k-t; and yellow denotes MB2 k-t velocities. Reconstructions are obtained using λ = 7 and σ = 10

F I G U R E 4 Simulation results
for MB3. A,B, Same as the conventions described in Figure 3A,B. To generate the synthetic data set, slices SL2, SL4, and SL6 were used from subject S5. Green lines highlight the positions along which data are displayed in M visualization mode (space vs time) in Figure 5. All data are shown before POCS reconstruction is applied MB2 data compared with the SB k-t equivalents (and 13 of 48 combinations for the case MB3).

| DISCUSSION
In this study, autocalibrated MB CAIPIRINHA with k-t acceleration were combined for highly accelerated TPM acquisitions with increased spatial coverage. A novel multiscale k-t SENSE reconstruction framework is proposed, based on a conjugate gradient that solves simultaneously for both MB and k-t acceleration. The proposed acquisition/reconstruction framework was validated on simulations, and subsequently extended to quantify the myocardial tissue velocities in vivo.
The presented work is subject to some limitations. First of all, an 18-heartbeat period is still too long for reliable use in patients, and the sequence is not compatible with retrospective cardiac triggering. For simplicity, in this study the slice gap was fixed to 20 mm, even though this might be suboptimal in some cases. Asymmetric echo has been used in this work for the MB acquisitions to reduce TR and improve the temporal resolution of the reconstructed time series. Although asymmetrical echo acquisitions are used commonly in phase-contrast imaging, it has been shown that velocity measurements of turbulent flow can be perturbed. 39 Furthermore, we observed in our work that although the POCS reconstruction resulted in an consistent increase in the spatial resolution, its application appeared to partly denoise the phase data (see Supporting Information Figure S2 in comparison to Figure 8). Another limitation of this work is that the signal of the inflowing blood was not suppressed with a "black blood" sequence, which is typically used in cardiac TPM to reduce flow artifacts and improve blood/myocardium contrast. In fact, combining the saturation of the blood signal with a MB excitation is not straightforward and is still subject to ongoing research. For complete blood saturation not only the blood below/above the bands needs to be addressed, but also F I G U R E 7 In vivo results. Single-band k-t (top) and MB2 k-t reconstructions (bottom) of subject S2 at peak diastole. First column reports magnitude data and, from left to right, the velocity maps v x , v y , and v z . All velocity maps are shown before eddy current correction is applied. Green lines highlight the positions along which data are displayed in M visualization mode (space vs time) in Supporting Information Figure S3. All data are shown after POCS reconstruction is applied F I G U R E 8 In vivo results. Same as in Figure 7 for the case of MB3 (subject S3 at diastole). Green lines highlight the positions along which data are displayed in M visualization mode (space vs time) in Supporting Information Figure S4. All data are shown after POCS reconstruction is applied blood in between the bands. Although MB excitation in combination with conventional black-blood imaging has been realized for T 2 mapping, 40 MB black-blood saturation pulses have not been demonstrated so far to the best of the authors' knowledge. As the aim of this work was to investigate the combination of in-plane k-t and through-plane MB acceleration, and its potential to reliably reproduce velocity times courses, the focus was set on the reconstruction obtained using rather moderate velocity-encoding values of 50 cm/s and 30 cm/s to limit flow artifacts. However, the work is expected to further benefit from additional blood saturation, and the investigation of combining both techniques remains to be demonstrated in future cardiac TPM studies.
In this study, an autocalibration scheme was used to estimate the static coil sensitivities. However, the explicit acquisition of such information is possible without a substantial increase in total scan time before/after the main acquisition block. Nevertheless, an autocalibration scheme such as the one proposed in this study has advantages, as (1) it increases patient comfort as no additional breathholds are required, (2) it has matched resolution between reference data and accelerated acquisitions, and (3) there is spatial alignment between the coil sensitivities and MB data.
For the reconstruction, the proposed algorithm extracts time-resolved calibration data to aid conditioning of successive iterations at finer temporal scales. A k-t SENSE regularization term (parameter λ) together with spatial filtering (parameter σ) of the calibration data are introduced. These parameters were chosen based on forward simulations of the imaging process. The optimizations were carried out such that temporal fidelity was maximized, leading to excellent RMSE and peak velocity reduction scores for MB2, and good reproducibility for MB3. However, this approach led to noisy reconstructions, especially when high accelerations were used (case MB3, Figure 3B). As shown in Supporting Information Figure S1, higher regularization terms can be used that yield images with lower noise, at the expense of smoothed velocity time courses with reduced peak velocities. It is clear that such level of trade-off depends on the specific application under consideration, and that individual optimizations often characterized by competing objectives/requirements might need to be carried out each time.
When comparing the actual accelerated SB k-t and MB k-t velocity time courses averaged within the myocardium, statistically significant differences were not detected. However, when performing the same analysis considering each individual AHA subregion, statistically significant differences were detected in 8% of all possible comparisons for the case of MB2, and 27% for MB3. The reasons for this regional variability can be attributed to various factors: single-band k-t and MB k-t data were in fact acquired F I G U R E 9 Motion traces for MB2 k-t. A, Slice-averaged longitudinal (v z ), radial (v r ), and tangential (v ) myocardial velocities as a function of the heartbeat duration for subject S1. Yellow time courses are measured with MB2 k-t, whereas blue plots represent the SB k-t equivalents acquired at the same slice locations as separate breath-holds. B, Bland-Altman plots of the SB versus MB2 velocity time courses across all subjects and slices. In this representation, different colors signify different subjects. Means ± twice the SDs of the differences Δv z , Δv r , and Δv are shown by solid and dashed lines, respectively. as separate breath-holds, and the AHA model masks were drawn independently on each data set. However, the increased variability between MB3 and MB2 and the noisy reconstructions of Figure 4A suggest that there is an upper limit beyond which acceleration is not possible. Although good reproducibility of the velocities was obtained, it appears that at this stage the MB3 k-t data reconstructed using the proposed framework are not clinically viable. There might be alternative acquisition/reconstruction frameworks to explore. For example, an attractive approach would consist of combining MB3 imaging with compressed sensing 41 acceleration. An initial study in the context of cardiac perfusion has suggested that this is feasible, 20 and the TPM framework so far developed for the myocardial velocities could be used as a benchmark for validation.

| CONCLUSIONS
A novel acquisition/reconstruction technique that combines MB imaging with CAIPIRINHA encoding and k-t acceleration for cardiac TPM was presented. The presented framework allows three-directional velocity encoding of the myocardial wall motion velocities at multiple slices in a single breath-hold.

SUPPORTING INFORMATION
Additional Supporting Information may be found online in the Supporting Information section.

FIGURE S1
Overregularized simulation results for multiband 3 (MB3) k-t. A, MB3 reconstructions. Top to bottom: Fully sampled (FS), single-band (SB) k-t, and MB3 k-t reconstructions from the synthetic experiment in the main manuscript shown at peak diastole (cardiac phase 13). In this case, the synthetic data set was formed by fully sampled SB data sets SL2, SL4, and SL6 from subject S5 (see section 2.2). Magnitude (left) and longitudinal velocitiesv z (right) are reported. B, Velocity traces averaged within the myocardium from Supporting Information Figure S1A as a function of the cardiac cycle. Light blue denotes FS; blue denotes SB k-t; and yellow denotes MB3 k-t velocities. Reconstructions are obtained using λ = 23 and σ = 20 FIGURE S2 In vivo results for MB3 k-t. Same data as in Figure 8 of the main manuscript before projection onto convex sets (POCS) reconstruction is applied FIGURE S3 M visualization mode (space vs time) for the data presented in Figure 7 in the main manuscript. Projections are taken along green lines in Figure 7 FIGURE S4 Same as Supporting Information Figure S3 for the data shown in Figure 8 in the main manuscript VIDEO S1 Data acquisition strategy as a function of consecutive heartbeats in MB2 (A) and MB3 (B). To ensure that the steady-state condition across consecutive heartbeats is fulfilled, the acquisition order of the k-space lines is chosen such that only identical RF pulses are transmitted within each heartbeat interval. For further information, refer to the caption of Figure 1 in the main manuscript