Quantification of balanced SSFP myocardial perfusion imaging at 1.5 T: Impact of the reference image

To investigate the use of a high flip‐angle (HFA) balanced SSFP (bSSFP) reference image (in comparison to conventional proton density [PD]–weighted reference images) for conversion of bSSFP myocardial perfusion images into dynamic T1 maps for improved myocardial blood flow (MBF) quantification at 1.5 T.


| INTRODUCTION
Cardiac MR (CMR) perfusion imaging is routinely used for ischemia assessment in patients with suspected coronary artery disease and is currently recommended by international guidelines. 1,2 A conventional CMR perfusion acquisition consists of a dynamic single-shot saturation-recovery (SR) sequence, acquired in 3-4 short-axis slices of the left ventricular myocardium. Visual assessment of first-pass perfusion images is typically performed to determine ischemic regions of the myocardium, which exhibit delayed and/or reduced signal enhancement relative to healthy regions of the myocardium. However, the diagnostic accuracy of visual assessment is dependent on the experience of the reader, 3 and can lead to underestimation of ischemic burden in the case of microvascular disease or multivessel coronary artery disease. 4 Therefore, quantitative assessment of myocardial blood flow (MBF) may be valuable, as it allows for an objective and absolute assessment of myocardial ischemia. 5 Cardiac MR MBF quantification has undergone significant development over the last few decades, and recent publications have shown that voxel-wise MBF quantification can lead to increased or comparable diagnostic accuracy over visual assessment. 3,5,6 The quantification of MBF using CMR has some benefits over PET, which is the current gold standard for perfusion quantification. Cardiac MR perfusion does not expose the patient to ionizing radiation, and high in-plane spatial resolution (< 2 × 2 mm 2 ) is attainable, which improves the detection of subendocardial defects and enables derivation of transmural perfusion gradients. 7,8 However, high-resolution SR images can suffer from low SNR. The relatively high myocardial SNR of balanced SSFP (bSSFP) perfusion imaging at 1.5 T is therefore advantageous to maximize the precision of resultant MBF maps. 9,10 Voxel-wise MBF quantification requires deconvolution of the arterial input function (AIF) (often selected as left-ventricular blood pool signal) and (voxelwise) myocardial enhancement curves. A number of parametrized models have been proposed to constrain the deconvolution. [11][12][13][14] The deconvolution process assumes that AIF and myocardial enhancement curves are linearly related to gadolinium concentration (Gd). However, the dose required for sufficient myocardial CNR results in a nonlinear relationship between signal enhancement and Gd. 15 For high Gd, such as the peak Gd typically reached in the blood pool, the significant nonlinearity results in signal saturation or "clipping." To avoid signal saturation of the AIF, dual-bolus 16,17 and dual-sequence 18 protocols are typically used for quantitative CMR perfusion imaging. Any remaining nonlinearity of the AIF and myocardial enhancement curves can then be corrected by using signal-calibration techniques to determine dynamic T 1 maps, which can then be converted to Gd using the relaxivity of the contrast agent. Cernicanu et al developed such an approach by deriving a signal equation for a gradient-echo (GRE) readout scheme and correcting for coil sensitivity and equilibrium magnetization with a proton density (PD)-weighted image acquired at the beginning of the sequence. 19 This approach has been extended to hybrid EPI 20 and bSSFP readout schemes. 21 For bSSFP, a MBF intersegment variability (0.14 ± 0.09 mL/g/min vs PD-GRE: 0.21 ± 0.09 mL/g/min and PD-bSSFP: 0.20 ± 0.10 mL/g/min; p ≤ .046).

Conclusion:
We have demonstrated the feasibility of using a HFA-bSSFP reference image for MBF quantification of bSSFP perfusion imaging at 1.5 T. Results from simulations demonstrate that the HFA-bSSFP reference image results in improved precision and reduced sensitivity to effective FA compared with conventional techniques using a PD reference image.
Preliminary in vivo data acquired at rest also demonstrate improved precision and intersegment variability using the HFA-bSSFP technique compared with PD techniques; however, a clinical study in patients with coronary artery disease under stress conditions is required to determine the clinical significance of this finding.

K E Y W O R D S
myocardial blood flow, myocardial perfusion, perfusion quantification, T 1 mapping look-up table is generated using Bloch simulations to model the signal evolution for a range of Gd. The normalized SR signal intensity values are then matched to the signal dictionary to estimate Gd for each dynamic pixel value. 21 The PD image used for normalization of the bSSFP perfusion images may be either a low flipangle (FA) bSSFP or low-FA GRE image, 21 with the latter demonstrating improved homogeneity of MBF maps at 3 T. 22 Most of the previous work evaluating techniques used for nonlinearity correction has focused on the accuracy of AIF and tissue Gd curves. 19,21,22 Because deconvolution is an inherently noise-sensitive technique, it is also imperative to optimize the SNR of the myocardial curves. Although spatial and temporal filtering can improve the SNR of the myocardial enhancement curves, this can result in the loss of information on the presence and extent of ischemia. 23,24 In this work, a high-FA (HFA) bSSFP reference image is proposed for conversion of bSSFP perfusion images into dynamic myocardial T 1 maps, for MBF quantification at 1.5 T. This is in contrast to previous methods, which use a PD reference image for normalization to ensure that the signal intensity is independent of T 1 . 21,22 Therefore, the T 1 independence of the proposed HFA-bSSFP method will be evaluated to ensure there is no significant bias introduced, compared with conventional PD techniques. It is hypothesized that using a HFA reference image (with the FA matched to the SR images) will provide improved precision and reduced intersegment variability compared with conventional PD reference-image techniques, due to the higher SNR 25 and the more comparable off-resonance response profiles of the HFA reference and SR images, respectively. Part of this work has been presented as conference proceedings. 26

| Pulse sequence
The prototype dual-sequence perfusion acquisition used to compare the proposed (HFA-bSSFP) and conventional PD reference-image techniques in vivo is illustrated in Figure  1. The sequence consists of a modified reference image acquisition followed by a conventional SR dual-sequence perfusion acquisition. The reference image acquisition includes a single low-resolution AIF reference image (PD-AIF) and three reference images for each high-resolution myocardial slice: a low-FA (5°) PD-GRE image, a low-FA (8°) PD-bSSFP image, and a HFA (50°) bSSFP image. Each slice is acquired separately for the HFA-bSSFP acquisition to minimize cross-talk. The PD-AIF, PD-GRE, PD-bSSFP, and HFA-bSSFP readouts are preceded by a rest period of 3, 3, 3, and 5 heartbeats, respectively. These rest periods were chosen based on Bloch simulations (using native myocardial T 1 time = 1200 ms and rest heart-rate range of 60-100 bpm) to achieve near-full recovery (> 97%) of longitudinal magnetization before each referenceimage acquisition and minimize bias between methods. During each RR interval of the conventional perfusion acquisition, a low-resolution SR image is acquired with a GRE readout for estimation of the AIF, and three highresolution SR images are acquired with bSSFP readout for estimation of voxel-wise myocardial enhancement curves for three slices. , and high flip-angle [HFA] bSSFP acquisitions following 3-rest, 3-rest, 3-rest, and 5-rest heartbeats, respectively, to allow for magnetization recovery), followed by a conventional saturation-recovery dual-sequence perfusion acquisition. Images with a red frame signify arterial input function (AIF) acquisitions, and "…" symbolizes multiple heartbeats

| Conversion of signal intensity to T 1 and Gd
The T 1 fitting software was developed in MATLAB (MathWorks, Natick, MA) using a signal dictionary approach. The signal dictionary is created for a range of T 1 values, and an exhaustive search of the dictionary is performed using least-squares fitting to find the best T 1 match, as described subsequently.

| Signal dictionary creation
A signal dictionary is created using Bloch simulations to model the signal evolution for a pair of reference and SR images. The signal of the reference image is the same for all dictionary entries, assuming typical native myocardial T 1 / T 2 times at 1.5 T of 1200 ms/50 ms, whereas the signal of the SR image is generated for a range of myocardial T 1 times (0-1200 ms) in steps of 0.1 ms, reflecting a specific Gd in the myocardium. To improve accuracy of the signal dictionary, a Gd-specific T 2 time is used for each dictionary entry. The T 2 is calculated according to Eqs. 1 and 2 using the native myocardial T 1 / T 2 times defined previously and R 1 /R 2 relaxivities of 4.7/6.8 L/mmol-s for the contrast agent (corresponding to gadobutrol [Gadovist; Bayer, Leverkusen, Germany] 27 ): The recorded signal value for each image in the signal dictionary is the absolute value of the transverse magnetization at the center of k-space. The signal dictionary entries represent values normalized to the equilibrium magnetization (ie, an initial longitudinal magnetization of 1 was used).

| T 1 fitting and conversion to Gd
The fitting process is summarized in Figure 2. Because the signal dictionary is normalized, each pair of measured reference and SR signal values are first individually scaled to each signal dictionary entry (similar to Refs 28 and 29), using the following scaling factor (Sf ): where SI Meas,Ref and SI Meas,SR represent the measured signal intensity of the reference image and the SR image, respectively; and SI Dict,Ref and SI Dict,SR represent the modeled signal intensity of the reference image and the SR image, respectively. The scaled signal-intensity values are then calculated as follows: The L2-norm between SI Scaled Meas and SI Dict is calculated for each dictionary entry, to determine the best match. The corresponding T 1 estimate is then converted to Gd using Eq. 1.

| Myocardial blood flow quantification
Myocardial blood flow is estimated using Fermiconstrained deconvolution, restricted to the first pass in the blood pool, as described previously 11,30,31 using inhouse software developed in MATLAB. The Gd estimates from the AIF were estimated using a similar approach as described previously for the myocardium but considering a typical native-blood T 1 at 1.5 T of 1600 ms. No spatial smoothing nor temporal smoothing was applied to the data before deconvolution.

| Numerical simulations
Numerical simulations were carried out in MATLAB to characterize accuracy and precision of dynamic T 1 and MBF estimates using the proposed HFA-bSSFP and conventional PD-GRE and PD-bSSFP techniques. All Bloch simulations used the nominal FAs, TR, TE, saturation time, and echo train length listed in section 2.4.1.

| Characterization of dynamic T 1 estimates
For T 1 simulations, a signal dictionary was first generated for fitting with the following fixed parameters: effective FA = nominal FA, native myocardial T 1 = 1200 ms, native myocardial T 2 = 50 ms, saturation efficiency (SE) = 1, with SE defined as where Mz res is the residual longitudinal magnetization immediately after the saturation pulse is applied, and Mz 0 is the normalized equilibrium magnetization (= 1). Signal values of reference/SR images were then simulated for a range of dynamic myocardial T 1 times (100-1200 ms in steps of 20 ms) and for the parameter ranges described subsequently, to characterize accuracy and precision of dynamic T 1 estimates.

Accuracy study
Simulations were performed for a range of effective FA of the excitation pulse (0.7-1.0 of nominal FA in steps of 0.005), native myocardial T 1 (600-1600 ms, in steps of 20 ms), native myocardial T 2 (30-70 ms, in steps of 1 ms), and SE (0.9-1.1, in steps of 0.001). Each pair of reference/ SR myocardial signal values generated was fitted to the signal dictionary to generate a T 1 estimate, as described in section 2.2.2. T 1 accuracy was calculated as the difference between estimated and reference dynamic T 1 values.

Precision study
Monte-Carlo simulations (N = 10 000) were performed to study the precision of dynamic T 1 estimates over a range of SNR levels corresponding to the expected SNR range in vivo of the HFA-bSSFP image (40-80). Rician noise was added to the simulated signal values. Each pair of noisy reference/SR signal values was fitted to the signal dictionary to generate a T 1 estimate (section 2.2.2). T 1 precision was calculated as the SD (over N) of the estimated T 1 values.

| Characterization of MBF estimates
To characterize MBF estimates, realistic AIF and myocardial enhancement curves (in units of Gd) were first generated as detailed in the Appendix for a range of MBF values (1-5 mL/g/min in steps of 0.05 mL/g/min), covering both stress and rest conditions. Each data point on each myocardial Gd curve was then converted to T 1 using Eq. 1 to create reference myocardial T 1 curves for each MBF value.

Accuracy study
Using the reference myocardial T 1 curves, Bloch simulations were used to generate a pair of reference/SR signal values for each T 1 data point. These simulations were performed for a range of effective FA of the excitation pulse, native myocardial T 1 , native myocardial T 2 , and SE (using the same ranges as for the T 1 accuracy study). A simulated myocardial T 1 curve was then estimated from the simulated signal values (ie, one T 1 estimated for each pair of simulated reference/SR signal values). These myocardial T 1 curves were then converted into dynamic Gd curves and used to calculate an estimated MBF value as described in section 0. The MBF accuracy was calculated as the difference between the reference and estimated MBF.

Precision study
Monte-Carlo simulations (N = 10 000) were performed to study the precision of MBF estimates, similarly as described for dynamic T 1 estimates. Briefly, each data point on each reference myocardial T 1 curve was converted to F I G U R E 2 For each pixel of each saturation-recovery image, an exhaustive search over the signal dictionary is performed to estimate T 1 on a pixel-by-pixel basis using least-squares fitting. The result is a "dynamic" T 1 map for each saturation-recovery image. T 1 values are then converted to gadolinium concentration (Gd) to compute a Gd map. The dynamic T 1 and Gd maps shown above represent the peak myocardial enhancement frame a pair of reference/SR signal values using Bloch simulations. Rician noise was added to each pair of signal values, followed by T 1 fitting, Gd conversion and MBF fitting, as described previously. The MBF precision was measured as the SD (over N) of the MBF estimates.

| In vivo evaluation
All imaging studies were performed at 1.5 T (MAGNETOM Aera; Siemens Healthcare, Erlangen, Germany; VE11C software) with a 32-channel spine array and an 18-element body coil.

| Data acquisition
The pulse sequence described in section 0 was acquired at rest in 8 patients (6 male, 2 female, mean age 60 ± 10 years) referred for a clinical contrast-enhanced CMR scan, following approval by the National Research Ethics Service (15/NS/0030) and written informed consent. The clinical indications for the scan included assessment of left-ventricular function, 2 etiology of cardiomyopathy, 2 exclusion of secondary causes of hypertension, 2 morphology assessment before atrial ablation, 1 and myocardial inflammation. 1 The three short-axis slices were planned on the systolic phase of 2-chamber, 3-chamber, and 4-chamber cine images to ensure coverage of the base, mid, and apical regions. The sequence was acquired in free breathing, including all reference and SR images. A 0.075-mmol/kg dose of gadobutrol (Gadovist) was administered at a rate of 4 mL/s and followed by a 25-mL flush of normal saline, after all reference images had been acquired. A sixpulse train saturation pulse optimized for the range of B 0 and B 1 inhomogeneity expected in the left ventricle at 1.5 T 32 was used to ensure optimal saturation efficiency. Offline motion correction was performed using precompiled (prototype) C++ code, provided by the MR scanner manufacturer. 33 An SR single-shot acquisition (SASHA) sequence was acquired in a breath-hold after the perfusion sequence (following a ~15-second pause to ensure full magnetization recovery), to provide reference postcontrast myocardial T 1 times as close as possible to the end of the perfusion scan. The following sequence parameters were used for SASHA: TR/ TE/ α = 2.7 ms/1.1 ms/70°, FOV = 360 mm × 270 mm 2 , slice thickness = 10 mm, acquired in-plane resolution = 1.4 × 1.8 mm 2 , in-plane acceleration = GRAPPA 2, partial Fourier = 7/8. Slice geometry was matched to the midslice of the perfusion acquisition.

| Dynamic T 1 fitting and MBF quantification
Each pair of reference/SR images was converted into a dynamic T 1 map and Gd map as described in section 2.2.2. The AIF images were converted to Gd maps using a similar approach but considering a native-blood T 1 time of 1600 ms. To generate the AIF Gd curve used for deconvolution, the AIF Gd at each time point was averaged over a region of interest manually drawn in the center of the blood pool, ensuring to avoid papillary muscles and partial-volume effects. An MBF map was estimated by deconvolution of the Gd curve for each myocardial voxel and the AIF Gd curve, as described in section 2.2.3. Any data sets in which the motion correction failed to align the reference images with the first-pass dynamics were excluded from T 1 precision, MBF precision, and MBF intersegment variability analysis, and any data sets in which the motion correction failed to align the reference images with the final 30 images of the series were excluded from the T 1 accuracy analysis.

| T 1 accuracy analysis
The myocardium was initially segmented (using a 16-American Heart Association segments model 34 ) on the peak left-ventricular blood pool frame of each motioncorrected perfusion data set, taking care to exclude pixels at the endocardial/epicardial borders to avoid corruption from partial-volume effects and residual motion across the dynamic series. To account for the delay between the end of the perfusion sequence and the SASHA acquisition, a prediction of the delay-corrected T 1 was calculated as follows: First, linear regression was used to determine the slope of average segmental T 1 as a function of time using the final 30 dynamics of the series. Then a delaycorrected T 1 for each myocardial segment was calculated as the extrapolated T 1 value corresponding to the time of the SASHA acquisition. Accuracy of T 1 was calculated as the difference between the average segmental dynamic T 1 and average segmental SASHA T 1 .

| T 1 and MBF precision analysis
T 1 precision analysis was performed using 11 timeframes corresponding to the first pass of the Gd in the myocardium (identified as the timeframe with peak Gd in a septal region of interest ± five timeframes).The intrasegment SD of T 1 estimates was computed for each myocardial segment and for each of the 11 dynamics. For each myocardial segment, the average intrasegment SD across the 11 dynamics was calculated as a surrogate measure of T 1 precision. The intrasegment SD of MBF estimates was calculated for each myocardial segment as a surrogate measure of MBF precision.

| Statistical analysis
The Friedman test was used to determine whether there was a significant difference in accuracy, precision, and spatial variability among the three techniques used for T 1 and MBF fitting. Where a significant difference was found, each pair of techniques were then compared using the Wilcoxon signed-rank test. For all statistical tests, a pvalue < .05 was considered significant.

| Precision
Results from T 1 and MBF precision simulations are presented in Figures 5 and 6, respectively. Precision of T 1 estimates was improved using the HFA-bSSFP technique: For a typical peak T 1 range of 200-400 ms, the SD of T 1 estimates was 30%-80% higher for the PD-GRE technique and 30%-90% higher for the PD-bSSFP technique compared with the HFA-bSSFP technique ( Figure 5D,E, respectively). Precision of MBF estimates was also improved using the HFA-bSSFP technique: The SD of MBF values was 20%-30% higher for the PD-GRE technique and 25%-35% higher for the PD-bSSFP technique compared with the HFA-bSSFP technique for typical rest MBF (1 mL/g/ min) ( Figure 6D,E, respectively), whereas for typical stress MBF (3-5 mL/g/min), the SD of MBF values was 50%-100% higher for the PD-GRE technique and 60%-115% higher for the PD-bSSFP technique, compared with the HFA-bSSFP technique ( Figure 6D,E, respectively). The relative improvement in T 1 and MBF precision using the HFA-bSSFP technique compared with the PD-GRE and PD-bSSFP techniques was independent of SNR, for the range of SNR values investigated.

| In vivo evaluation
No patients had any myocardial segments with late gadolinium enhancement. Motion correction failed for the wash-out period of one patient data set, which was excluded from the T 1 accuracy analysis. The motion correction failed for the first pass of another patient data set,

F I G U R E 4 Simulations of myocardial blood flow (MBF) estimation. A-C,
The MBF estimation-error maps due to discrepancy between signal dictionary excitation (nominal FA) and effective excitation FA (70%-100% of nominal FA). D-F, The MBF estimation-error maps due to discrepancy between signal dictionary saturation efficiency (1.0) and ground-truth saturation efficiency (0.9-1.1). All maps are shown in units of milliliters per gram per minute which was excluded from MBF spatial variability and T 1 / MBF precision analysis. Across all remaining patients, dynamic T 1 and MBF maps were generated successfully. An example case is shown in Figure 7, demonstrating high quality and improved precision for T 1 , Gd, and MBF maps using the HFA-bSSFP technique compared with the PD-GRE and PD-bSSFP techniques. For this patient, the measured T 1 /MBF intrasegment SD was 34 ms/0.21 mL/g/min for the HFA-bSSFP technique, 40 ms/0.27 mL/g/min for PD-GRE, and 48 ms/0.31 mL/g/min for the PD-bSSFP technique.

| DISCUSSION
In this work, a HFA-bSSFP reference image is proposed for conversion of bSSFP perfusion images into dynamic myocardial T 1 maps. The results from simulations demonstrate that this approach results in improved precision of T 1 and MBF estimates compared with conventional techniques, with lower sensitivity to effective excitation FA.
Results from in vivo studies also demonstrated improved precision of dynamic T 1 and MBF estimates using the proposed technique, as well as reduced spatial variability of MBF estimates. One of the primary motivations for using a HFA-bSSFP reference image is the increased SNR of this image. 25 The results from simulation and in vivo studies demonstrate that this approach results in improved precision of dynamic T 1 and MBF estimates compared with conventional techniques, using a PD-GRE or PD-bSSFP reference image for signal-intensity normalization. The noise enhancement resulting from the low SNR of reference images conventionally used for signal-intensity normalization has been acknowledged elsewhere. 19,35 This is typically overcome by using spatial smoothing of the reference image and SR images to improve the precision of MBF estimates. Indeed, spatial smoothing can be combined with the proposed technique to further improve precision of MBF estimates. However, this requires a further processing step and may result in undesirable spatial smoothing at the edges of the myocardium. For the current study, the FA of the HFA reference image was selected to match the FA of the SR images and the associated off-resonance response profile. In the case in which a different imaging FA is used, this approach of using a matched FA for the reference image may be desirable to minimize sensitivity to off-resonance.
An advantage of the conventional PD reference image is that the signal is largely independent of this baseline T 1 . Therefore, it was important to evaluate whether the proposed HFA technique also demonstrated this independence of baseline T 1 and used imaging parameters. This was confirmed using Bloch simulations, which demonstrated that the HFA-bSSFP technique is insensitive to baseline T 1 for the range of expected baseline T 1 values. A native T 1 map could also be acquired before the perfusion acquisition to provide an estimate of baseline T 1 for the signal dictionary.
Simulations also demonstrated that the HFA-bSSFP technique is less sensitive to effective excitation FA. Spatial variation in the effective FA in vivo is expected due to B 1 variation across the myocardium (ranging between ~0.8 and 1.0 at 1.5 T 36 ); therefore, the sensitivity of conventional PD reference-image techniques to effective FA may compromise the accuracy of dynamic T 1 and MBF estimates in some regions of the myocardium. The results of our in vivo study support this hypothesis, as the intersegmental variation of MBF is F I G U R E 7 Example images from a patient study showing reference images, SR images, dynamic T 1 , and Gd maps at peak myocardial Gd and 10 frames after peak Gd, and the resultant MBF maps at rest. The epicardial contour used for analysis is shown on each map/image F I G U R E 8 Midslice segmentwise error in dynamic T 1 estimates in milliseconds averaged over 7 patients. A, Difference between average dynamic T 1 estimate over 30 final frames of the acquisition and SR single-shot acquisition (SASHA) T 1 . B, Difference between extrapolated dynamic T 1 estimate and SASHA T 1 . Center value represents mean error over all segments. Abbreviation: LFA, low flip angle higher for both PD techniques compared with the HFA-bSSFP technique.
A previous study compared the PD-bSSFP and PD-GRE approach for signal-intensity normalization of bSSFP perfusion images at 3 T, 22 concluding that the PD-bSSFP approach results in increased spatial variability of quantitative and semi-quantitative measures of MBF due to increased sensitivity to off-resonance. Our study did not find any significant difference between these techniques in terms of the intersegment spatial variability of MBF estimates, which may be due to the lower B 0 inhomogeneity expected at 1.5 T, combined with high-order shimming and careful adjustment of the shim volume.
It is noted that physiological variations within and between myocardial segments may contribute to the intrasegment SD and intersegment SD of MBF estimates, respectively. Physiological variations in MBF are expected to be small within each myocardial segment (which can be assumed to be supplied by a single epicardial artery) and between myocardial segments (as acquired at rest in patients with no late gadolinium enhancement). More importantly, any contributions of physiological variations in MBF to intrasegment SD or intersegment SD will be the same for all three techniques. Therefore, the impact of physiological variations on the relative comparison of precision and intersegment variations of MBF among the three techniques should be minimal.
Conventionally CMR perfusion protocols require a perfusion sequence acquired at stress followed by a perfusion sequence acquired at rest. Therefore, a residual level of Gd remains in the myocardium before the rest acquisition, resulting in a lower baseline T 1 . An advantage of the conventional PD reference image is that the signal is largely independent of this baseline T 1 . However, results from our simulations have demonstrated that signal dictionarybased T 1 fitting using a HFA reference image is also insensitive to T 1 for the range of T 1 values expected at baseline. A native T 1 map could also be acquired before the perfusion acquisition to provide an estimate of baseline T 1 for the signal dictionary.
The results from the in vivo study showed an increase in dynamic T 1 precision using the HFA-bSSFP technique compared with the PD-GRE technique (intrasegment SD was 16% higher with PD-GRE than HFA-bSSFP); however, the improvement in precision for HFA-bSSFP versus PD-GRE was lower than predicted by simulations (SD was 30%-80% higher with PD-GRE than the HFA-bSSFP technique). It is important to note that the in vivo precision reported is in fact the intrasegment SD of T 1 and MBF. Although the inherent precision of the technique (due to the presence of noise) will contribute to intrasegment SD, there are also contributions from residual motion, physiological variation within each myocardial segment, and artifacts. Therefore, reduction in the intrasegment SD using a HFA-bSSFP technique versus a PD-GRE technique will be less significant than the reduction in the inherent precision, as modeled in simulations.
A number of rest heartbeats were included in the reference-image acquisition section of the sequence to ensure full magnetization recovery before each subsequent reference-image acquisition. The maximum heart rate in our study was 80 bpm, which would result in a maximum T 1 underestimation bias of 2% across all techniques. It is noted that these rest heartbeats would not be required in a clinical application of the sequence, where only one type of reference image is acquired. For the HFA-bSSFP reference-image acquisition, each slice was acquired in a separate heartbeat. The acquisition was designed in this F I G U R E 9 Intersegment SD of MBF in milliliters per gram per minute. A, The SD averaged over all segments and all slices for each patient. B, Boxplots showing patient-wise distribution of SDs averaged over all segments and all slices. On each plot, the red bar represents the median value; the black whiskers represent the minimum and maximum values; and the blue box represents the interquartile range (IQR). The "+" symbols indicate outlier data points (ie, points that are either greater than [Q3 + 1.5 × IQR], or less than [Q1 − 1.5 × IQR]. Significant differences are indicated by an asterisk (p < .05) way to achieve optimal accuracy, even in the worst-case scenario (i.e., breathing motion resulting in two successively acquired slices overlapping). It is acknowledged that acquiring each reference image in a separate heartbeat is less efficient, and may not be necessary to avoid cross-talk, given that there is typically a 10-20-mm slice gap between successive slices. Therefore, this is a conservative strategy, and further investigation into the probability and magnitude of cross-talk effects could be performed to remove this requirement and enable acquisition of all three slices within the same heartbeat.
It is noted that magnetization-transfer effects are not modeled in the Bloch simulations used to create the signal dictionary. Magnetization-transfer effects result in a reduction in signal intensity, particularly at HFAs 37 ; however, the scaling factor used for dictionary matching is expected to compensate to some extent for the mutual reduction in signal intensity of the HFA-bSSFP and SR images. Further studies are needed to fully characterize the effect of magnetization transfer on perfusion quantification using bSSFP imaging.
Our study was performed at 1.5 T, which is currently the most common field strength used for CMR. At 3 T, relaxation times and relaxivity rates will be different, and B 1 inhomogeneity is expected to be more significant. Therefore, further evaluation of the technique would be necessary to determine whether the proposed method is generalizable to 3 T.
There are some limitations of the current study. It was not possible to determine the in vivo accuracy of MBF estimates using HFA-bSSFP, PD-bSSFP, and PD-GRE techniques, as no gold-standard assessment of MBF exists F I G U R E 1 0 Precision of T 1 (calculated as intrasegment SD averaged across 10 dynamics during the first pass) and MBF (calculated as intrasegment SD) measured in 7 patients. T 1 (A) and MBF (B) SD averaged over all segments and all slices for each patient. T 1 (C) and MBF (D) boxplots showing patientwise distribution of SD averaged over all segments and all slices. For a description of the boxplot characteristics, see Figure  9. T 1 (E-G) and MBF (H-J) bullseye plots showing SD averaged over all patients for each segment. The value at the center of each bullseye plot represents the average SD over all segments for CMR perfusion; however, all techniques resulted in similar MBF estimates, allowing comparison of the relative precision and spatial variability of the techniques. It was not possible to compare T 1 estimates at peak Gd, as there is no gold standard T 1 measurement technique for a dynamically changing contrast concentration; however, results from simulations suggest that the HFA-bSSFP technique is accurate for T 1 estimates across the range expected in the myocardium in vivo, and this is further validated by the postcontrast in vivo comparison with SASHA. Although this study demonstrated a significant improvement in precision and reduced spatial variability of MBF estimates using a HFA-bSSFP technique, further evaluation of the technique in a large cohort of patients with ischemia and under stress conditions would be required to assess the clinical significance of these benefits.

| CONCLUSIONS
We have demonstrated the feasibility of using a HFA-bSSFP reference image for MBF quantification of bSSFP perfusion imaging at 1.5 T. Results from simulations demonstrate that the HFA-bSSFP reference image results in improved precision and reduced sensitivity to effective FA compared with conventional techniques using a PD reference image. Preliminary in vivo data acquired at rest also demonstrate improved precision and intersegment variability using the HFA-bSSFP technique compared with PD techniques; however, a clinical study in patients with coronary artery disease under stress conditions is required to determine the clinical significance of this finding.

SUPPORTING INFORMATION
Additional supporting information may be found in the online version of the article at the publisher's website. FIGURE S1 Simulations of T 1 estimation. A-C, T 1 estimation error maps due to discrepancy between signal dictionary native T 1 (1200 ms) and ground-truth native T 1 (600-1600 ms). Grayed-out region signifies invalid combination (dynamic T 1 > native T 1 ). D-F, T 1 estimation-error maps due to discrepancy between signal dictionary-native T 2 (50 ms) and ground-truth native T 2 (30-70 ms). All maps are shown in units of milliseconds FIGURE S2 Simulations of myocardial blood flow (MBF) estimation. A-C, Myocardial blood flow estimation-error maps due to discrepancy between signal dictionary native T 1 (1200 ms) and ground-truth native T 1 (600-1600 ms). D-F, Myocardial blood flow estimation-error maps due to discrepancy between signal dictionary-native T 2 (50 ms) and ground-truth native T 2 (30-

APPENDIX
To generate realistic tissue-response curves for myocardial blood flow (MBF) fitting, a physiological impulse response function was simulated. First, the tissue-residue function was generated using the Fermi function. 1,2 The fitting parameters were set to ensure that the area under the curve (corresponding to the mean transit time) was equal to 6 seconds, which falls within the typical physiological range reported in the literature, 1,3,4 and its value at t = 0 was equal to 1 (indicating that all tracer molecules remain in tissue at t = 0). This residue function was then scaled by the required MBF to calculate the response function as follows: where R F (t) is the response function, and R (t) is the residue function. Finally, a representative myocardial tissue curve in gadolinium concentration units was generated by convolving the arterial input function gadolinium concentration curve from a representative patient data set with this response function.