- Open Access
Myocardial perfusion cardiovascular magnetic resonance: optimized dual sequence and reconstruction for quantification
Journal of Cardiovascular Magnetic Resonance volume 19, Article number: 43 (2017)
Quantification of myocardial blood flow requires knowledge of the amount of contrast agent in the myocardial tissue and the arterial input function (AIF) driving the delivery of this contrast agent. Accurate quantification is challenged by the lack of linearity between the measured signal and contrast agent concentration. This work characterizes sources of non-linearity and presents a systematic approach to accurate measurements of contrast agent concentration in both blood and myocardium.
A dual sequence approach with separate pulse sequences for AIF and myocardial tissue allowed separate optimization of parameters for blood and myocardium. A systems approach to the overall design was taken to achieve linearity between signal and contrast agent concentration. Conversion of signal intensity values to contrast agent concentration was achieved through a combination of surface coil sensitivity correction, Bloch simulation based look-up table correction, and in the case of the AIF measurement, correction of T2* losses. Validation of signal correction was performed in phantoms, and values for peak AIF concentration and myocardial flow are provided for 29 normal subjects for rest and adenosine stress.
For phantoms, the measured fits were within 5% for both AIF and myocardium. In healthy volunteers the peak [Gd] was 3.5 ± 1.2 for stress and 4.4 ± 1.2 mmol/L for rest. The T2* in the left ventricle blood pool at peak AIF was approximately 10 ms. The peak-to-valley ratio was 5.6 for the raw signal intensities without correction, and was 8.3 for the look-up-table (LUT) corrected AIF which represents approximately 48% correction. Without T2* correction the myocardial blood flow estimates are overestimated by approximately 10%. The signal-to-noise ratio of the myocardial signal at peak enhancement (1.5 T) was 17.7 ± 6.6 at stress and the peak [Gd] was 0.49 ± 0.15 mmol/L. The estimated perfusion flow was 3.9 ± 0.38 and 1.03 ± 0.19 ml/min/g using the BTEX model and 3.4 ± 0.39 and 0.95 ± 0.16 using a Fermi model, for stress and rest, respectively.
A dual sequence for myocardial perfusion cardiovascular magnetic resonance and AIF measurement has been optimized for quantification of myocardial blood flow. A validation in phantoms was performed to confirm that the signal conversion to gadolinium concentration was linear. The proposed sequence was integrated with a fully automatic in-line solution for pixel-wise mapping of myocardial blood flow and evaluated in adenosine stress and rest studies on N = 29 normal healthy subjects. Reliable perfusion mapping was demonstrated and produced estimates with low variability.
Myocardial perfusion can be evaluated with dynamic cardiovascular magnetic resonance (CMR) during the passage of a bolus of contrast agent. Most commonly, perfusion CMR is evaluated qualitatively, but objective quantitative evaluation would be more desirable. The potential benefits of quantification are: objective assessment, simpler and faster analysis, and the ability to detect disease with a global reduction in flow such as multi-vessel obstructive disease or microvascular disease. Quantification of myocardial blood flow using CMR was first proposed over 20 years ago [1, 2], yet qualitative interpretation of images remains the primary means available to clinicians. The desired output of a quantitative perfusion study is a map of myocardial blood flow in units of ml/g/min.
Quantification requires knowledge of the amount of contrast agent in the myocardial tissue and the arterial input function (AIF) driving the delivery of this contrast agent. Accurate quantification is challenged by the lack of linearity between the measured signal and contrast agent concentration. Ideally these measurements would consist of input and response curves in units of contrast agent concentration. Current, commercially available, myocardial perfusion sequences have not been optimized to yield accurate concentration curves and the observed signal intensity curves are not linearly related to the concentrations of contrast agent, i.e., there is a non-linear relationship between signal intensity and contrast agent concentration, which leads to quantification biases.
The main sources of non-linearity and bias are: spatial signal variations caused by the sensitivity profiles of the surface coils, imperfect saturation of magnetization during contrast bolus passage, T2* decay (and signal loss) caused by high contrast agent concentrations in the blood pool, and the non-linear signal response inherent due to saturation recovery that depends on the parameters of the imaging protocol. It has been proposed that some of the non-linearity of the AIF response curve can be mitigated by imaging the AIF during a separate injection of a bolus with lower concentration (the dual bolus approach), but this approach has some practical drawbacks as it requires multiple injections and acquisitions. Moreover, there are other potential bias sources with his approach, since changes in breathing, etc. between the two measurements may introduce new sources of variation. Consequently, it is desirable to acquire the AIF curve simultaneously with the tissue response curve.
A dual sequence  approach, which separately optimizes the imaging protocols for blood and myocardium has been proposed. This approach may be more easily incorporated into a clinical workflow. The proposed dual sequence was optimized for perfusion quantification and was evaluated using a recently developed fully automatic in-line solution for pixel-wise mapping of myocardial blood flow . This work characterizes the sources of non-linearity and presents a systematic approach to accurate measurements of contrast agent concentration in both blood and tissue of interest.
A saturation recovery (SR) sequence was used for myocardial perfusion imaging during the passage of a bolus of gadolinium based contrast agent as depicted in Fig. 1 which is illustrated for a subject with single vessel disease. Baseline images were acquired prior to bolus administration and continued through the first pass. Typically, images were acquired for 60-90 heartbeats depending on the cardiac output. Proton density (PD) weighted images were acquired at the start.
A multi-slice 2D SR dual imaging sequence is diagrammed in Fig. 2. Low resolution blood pool images used for estimating the AIF were acquired every heartbeat immediately following the R-wave trigger. Higher resolution images were acquired following the AIF and may be sampled every RR or every second RR if greater spatial coverage is desired. The sequence uses a pulse sequel for saturation  for each image. The image readout is single shot using parallel imaging acceleration to reduce the imaging duration. The AIF uses a FLASH readout, whereas the higher resolution myocardial images may be either b-SSFP or FLASH, selected by the user. The measurement begins with the acquisition of PD weighted images used for surface coil intensity correction and normalization of signal values. The PD images are acquired using a low flip angle FLASH readout without SR preparation to minimize artifacts of b-SSFP at low flip angle . An optional chemical shift fat suppression may be used to mitigate artifacts due to the presence of fat around the heart. Fat suppression is used in this study.
Saturation efficiency is very important in quantification since the conversion of signal intensities to gadolinium concentration depends on a known signal recovery and independence of the signal from slice to slice that is achieved by resetting the magnetization to zero for each image. With high saturation efficiency it is also possible to prescribe a mixture of long and short axis views without cross-talk between slices due to the readout. The SR preparation using the 6-pulse design  was chosen since it had excellent saturation efficiency over a wide range of off-resonance and effective transmitter flip angle (FA) which may vary across the heart. A BIR-4 design  achieved excellent saturation performance with a shorter duration but was found empirically to have specific absorption rate (SAR) limitation at higher heart rates particularly at higher field strength such as 3 Tesla. The 6-pulse design consisted of non-selective RF pulses with tailored flip angles separated by gradient spoilers. The voltage of the RF pulses was maximized in order to reduce the pulse duration, and the duration was modulated to achieve the specified FAs. The pulse amplitude corresponded to approx. 27 μT at 1.5 T and was reduced to approx. 11 μT at 3 T to reduce the SAR. To ensure that the sequence would not terminate at run time due to average SAR monitor responding to actual changes in the heart rate, a post-scan acquisition delay of up to 100 s was allowed to increase the averaging interval. The performance of the SR preparation in both blood and myocardium was characterized by simulation for different gadolinium contrast concentrations. At 1.5 T, the SR preparation was 26 ms including 1 ms pre- and 4 ms post-spoiler gradients.
Signal intensities were converted to gadolinium concentration, [Gd], in order to linearize the relationship of signal and [Gd] and to be able to have a common scaling between the AIF and myocardial signals which are acquired using different protocols. Cernicanu, et al.  proposed a method for conversion of normalized signal intensities using an analytic expression for readout using a gradient recalled echo (GRE) protocol. This formulation was extended to the dual sequence using a numerical Bloch calculation  which permitted application to b-SSFP readout of myocardium and FLASH readout of the AIF.
Proton density weighted images were acquired at the start of the scan for both AIF and myocardial image slices using a FLASH sequence without the SR preparation. The timing of the PD images matched the SR prepared images such that the images were acquired at the same cardiac phase. The images were used to correct the surface coil variation and were used as the signal normalization for look-up table (LUT) linearization. LUT calculations assumed that the native tissue T1 for the PD is prior to contrast (T10), therefore it was important that the PD signal intensity was sufficiently independent of the actual T1 since acquisition of rest perfusion scans typically follow the stress scans after only several minutes at which time the actual T1 is not fully recovered, i.e., the actual gadolinium concentration [Gd] > 0. The dependence of PD signal amplitude versus [Gd] for myocardium and blood tissue was calculated through simulation. This led to a selection of readout FA = 5°. Note that coil sensitivity maps, used both in parallel imaging reconstruction and adaptive coil combination, consisted of the average of all time frames including PD weighted images, therefore the number of PD frames acquired was set equal to the parallel imaging acceleration factor of the myocardial imaging (in this case, R = 3). The first PD image was used for normalization to avoid signal loss caused by previous heartbeat images.
Arterial input function
The AIF was acquired immediately after the R-wave trigger and was selected as the most basal of the slices prescribed in the first slice group. The AIF used the 6-pulse sequel for saturation preparation as described above followed by a dual echo low FA FLASH readout. The protocol parameters are listed in Table 1. A short readout (64 point) with wide bandwidth (3900 Hz/pixel) and short duration RF pulses (250 μs, time-bandwidth product = 2.0) were used to achieve low T2* losses (TE1 = 0.76 ms). T2* dephasing loss has been a known concern in estimating AIF and conversion to [Gd] and approaches to this problem have focused either on minimizing the loss by choosing adequately short echo time (TE)  or on correcting for T2* loss based on modeling the relationship between T1 and T2* [11, 12]. In this work, the dual sequence approach was modified to incorporate a 2 echo acquisition for measurement of T2* during the bolus passage. A dual echo acquisition with monopolar readout was used to acquire a second echo (TE2 = 1.76 μs) which was used for direct estimation of T2* during the first pass. The ratio of the 2 echo signals S1/S2 = S0 exp((T2-T1)/T2*) was used to calculate the signal amplitude S0 without T2* loss. This was performed for a blood pool region signal after left ventricle (LV) blood pool segmentation.
Blood pool segmentation was performed on the motion corrected low resolution AIF image series to extract arterial input function intensities signals for both echoes. First, the AIF PD image is used to detect the noise background. Since the noise standard deviation (SD) is unity after the SNR unit reconstruction , a simple threshold of 3 SD’s was used. For all foreground pixels as determined by the noise mask, the time intensity curves are analyzed using a scale-space based detector . Pixels with top 10% upslope and AUC values are picked as the candidates for LV blood pool mask. A connected component analysis is then used to separate RV and LV pixels based on the time to peak enhancement. The final LV blood pool mask is calculated using a further erosion step which seeks to drop border pixels which are a mixture of blood and myocardial tissue.
In order to shorten the imaging duration, 2-fold acceleration was achieved using parallel imaging with temporal generalized autocalibrating partially parallel acquisitions (TGRAPPA) . In order to minimize the non-linear response due to saturation at high gadolinium concentration a short saturation delay (TS) is desired, where TS is defined as the time from saturation to the k-space center. There is a tradeoff between the image signal-to-noise ratio (SNR), which is reduced at low TS, and linearity, which is improved at low TS. Early designs used a centric readout ordering to minimize TS [16, 17]. In centric ordering the saturation delay is nearly the trigger delay (TD), which may be as short as the gradient spoiler following the RF saturation (4 ms). The difficulty with centric order is 2-fold: 1) the SNR of the baseline images is too low, and 2) the k-space weighting during the saturation recovery leads to a strong high pass spatial filter that enhances the edge of the blood pool. The high pass spatial filter is problematic for accurate calculation of gadolinium concentration since the effective TS is a function of spatial frequency and will vary depending on how the blood region is segmented. For this reason, a linear acquisition ordering was chosen. Although not perfectly linear at high [Gd], it was sufficiently linear to enable LUT correction.
Using a low resolution image and parallel imaging factor 2, the k-space center was typically at N = 9 pulses. The sensitivity of the blood signal to in-flowing spins was estimated by comparing the signal after 9-pulses compared to the signal assuming all spins were new.
Conversion of the signal to gadolinium concentration, [Gd], was performed by LUT based on Bloch signal calculations. In this way, the LUT corrected signal was directly proportional to [Gd] and importantly was in the same units as the LUT corrected myocardial signal which was acquired with a different imaging protocol. The LUT was applied to the normalized signal SR/PD where SR and PD were the saturation recovery and proton density weighted images, respectively. It was important that the normalized signal SR/PD was not strongly dependent on the actual transmitted FA. The sensitivity of the LUT to transmit FA was calculated through simulation.
The readout of the PD image for the AIF may influence the initial magnetization of the 1st myocardial image which is at the same slice location since there is no SR preparation for the PD image. For this reason, a low PD FA is used (5°) which minimizes this effect. After 17 RF pulses (time bandwidth 2.0), the magnetization is reduced approximately 2%, as calculated by Bloch simulation.
The 2D multi-slice myocardial imaging sequence used the same 6-pulse sequel saturation preparation, followed by a trigger delay (TD) and single shot readout. The single shot readout was either FLASH or b-SSFP. Protocol parameters at 1.5 T are listed in Table 2 and may vary slightly at 3 T. The readout used parallel imaging with 3-fold acceleration using TGRAPPA, and there were 3 PD frames without SR using a FLASH readout at the start of the sequence. An optional chemical shift fat saturation may be used without any penalty in the timing since the TD accommodated the fat saturation RF pulse. Although the gadolinium concentration in the myocardium is typically < 1 mmol/L, the signal response is still somewhat non-linear and therefore, the normalized signal SR/PD was corrected by a LUT which converts the myocardial signal to gadolinium concentration units, [Gd]. The sensitivity of the LUT correction to the actual transmitted FA was calculated by simulation. The b-SSFP FA was limited to 50° in order to reduce sensitivity of LUT correction to variations in actual transmitted FA and well as to reduce the average SAR.
The duration of actual signal shot image was 70 ms for SSFP protocol using factor 3 acceleration and Partial Fourier factor of ¾ with the latter part of k-space omitted. There is a trade-off between contrast-to-noise ratio (CNR), spatial coverage (number of slices per RR), and linearity as illustrated in Fig. 3. The protocol was designed to work at a heart rate of 120 bpm which is commonly seen for patients under adenosine stress. It was possible to acquire the AIF plus 3 slices at 120 bpm with the proposed protocol using TS = 95 ms, or AIF plus 2 slices using TI = 160 with increased CNR. It is also possible to prescribe 2x the number of slices by acquiring slices at 2RR intervals. Although there is a gain CNR with longer TS, there is also a loss in performance when using 2RR sampling since there will be fewer samples of the myocardial signal during the first pass measurement.
Image reconstruction & flow estimation
Image reconstruction and processing steps are diagrammed in Fig. 4. Parallel imaging was used to accelerate the acquisition of both AIF and myocardial images. This helped to minimize cardiac motion contribution to dark rim artifacts (DRA) [18, 19] that appear as false perfusion defects, and to achieve adequate spatial coverage. Parallel imaging used TGRAPPA reconstruction with coil maps estimated by integrating the complete dataset . All acquisitions were with normal free breathing and parallel imaging auto-calibration was performed on the complete dataset resulted in images free of aliasing artifacts. Individual coil images are adaptively combined to minimize noise bias prior to magnitude detection.
Raw filtering was used to reduce edge ringing  (Gibb’s ringing) and to mitigate contamination of the measurements due to fat. In the low resolution AIF, a true Hanning window (i.e., without modification as in Tukey windows) was used to mitigate the influence of fat on the blood pool. For the low resolution AIF images the chest wall fat may be only 10–20 pixels from the LV blood pool. The Hanning window point spread function is <0.02% at a distance greater than 10 pixels whereas an un-windowed reconstruction would be <3.5% at this distance. The loss in resolution at full width half maximum is approx. 60%. In the higher resolution myocardial images, a truncated Gaussian filter was used, truncated to a width of 1.5 standard deviations. The first sidelobe was reduced to 36% compared to unweighted with a mainlobe broadening of 16% at full width half maximum.
Reconstruction and processing were implemented within the Gadgetron software framework  and were in-line and fully automatic. All images were respiratory motion corrected using a non-rigid image registration . All images were reconstructed in SNR units  to facilitate image scaling, SNR measurements, and calculation of fixed threshold noise masks.
Myocardial blood flow was calculated multiple tissue models: 1) a Fermi model , and 2) a blood tissue exchange (BTEX) model originally developed by Bassingthwaighte  which is a distributed model described by the partial differential equations (PDE):
where subscripts p and isf correspond to plasma and interstitial fluid, respectively, C the contrast agent concentration, F is the blood flow, PS is the permeability surface area product for the capillaries, Vp and Visf are the intracapillary plasma and interstitial fluid volume, respectively, D is the axial diffusion coefficient, L is the capillary length, and x is the distance along the capillary. These equations follow Bassingthwaighte  Eqs. (1) and (3) with the term for regional consumption ignored, assuming a single term for capillary leakage dominated by the gaps in the capillary wall, and that the gadolinium based contrast agent is extracellular. The BTEX implementation solved for 4 unknown parameters: myocardial blood flow, interstitial volume, plasma volume, and the permeability surface area product that governs the extraction efficiency, with fixed values for other parameters. The PDE was applied to the AIF to calculate the myocardial response for each set of model parameters, and the parameters with the minimum mean squared error were used as the estimate. The Fermi model was fit to the first pass only and the BTEX model was fit to the entire measurement. The influence of T2* correction of the AIF on myocardial blood flow was analyzed.
The sequence was simulated by Bloch equations to calculate the transverse magnetization as a function of all the protocol and tissue parameters in order to construct the LUT corrections for both the AIF and myocardial imaging protocols. Input to the LUT was the normalized signal SR/PD, where PD used the FLASH protocol and SR used either a b-SSFP or FLASH protocol. LUT were validated by phantom measurement by comparing the estimates of [Gd] after LUT correction with the known [Gd] using least squares fitting. A set of gadolinium doped saline phantoms were constructed at concentrations up to 10 mmol/L using both Gadoterate meglumine (Dotarem, Guerbet LLC) and Gadobutrol (Gadavist, Bayer Healthcare). LUT estimates of [Gd] vs known [Gd] were calculated with and without T2* correction. The phantom T1 values were measured using an inversion recovery GRE sequence at multiple inversion times (TI) with TR = 10s such that the longitudinal magnetization was fully relaxed after each RF excitation, and T1 was estimated by 3-parameter fitting to the mono-exponential inversion recovery S = A-Bexp(-TI/T1). The phantom T2 values were measured using a spin echo sequence (TR = 10s) with varying echo times (TE) and using T2 estimates from a 2-parameter fit to the mono-exponential decay curve, S = Aexp(-TE/T2). The coefficients for relaxivity rates (r1 and r2) were calculated from the T1 and T2 measurements vs known [Gd] using linear fitting, i.e., R1 = R10 + r1[Gd] with R1 = 1/T1.
In-vivo data measurements
The proposed sequence and in-line flow mapping was performed at stress and rest on 29 healthy normal volunteers (11 men and 18 women, mean age 25.4 ± 5.7 years) at the Karolinska University Hospital, Stockholm, Sweden. Studies were approved by the local Ethics Committee. Anonymized data was analyzed at NIH with approval by the NIH Office of Human Subjects Research OHSR (Exemption #13156). All imaging was performed at 1.5 T (Magnetom AERA, Siemens, software version VE11A). Gadolinium (Gd) contrast agent (Gadobutrol) was administered as a bolus with 0.5 dose (0.05 mmol/kg) at 4 mL/s with 20 mL saline flush. One cannula was used for administration of adenosine and another cannula for the administration of contrast agent. Adenosine was administered by continuous infusion for approximately 8 min at a dose of 140 μg/kg/min to allow for additional research scans at stress just prior to contrast administration. The SSFP protocol was used in this study with fat saturation enabled.
In-vivo studies were performed to test the sequence and LUT conversion of signal intensities. Peak [Gd] was measured for the AIF blood pool signal and myocardium, as well as peak SNR in the myocardium from SNR scaled signal intensities. Blood pool T2* values at peak [Gd] were measured as well as the influence of T2* correction on estimates of myocardial blood flow. Duration of the bolus first pass was measured automatically from the AIF signal from the foot of the curve on the upslope of the AIF to the foot of the downslope. The improvement in linearity of the AIF after conversion to gadolinium concentration was measured by the ratio of the AIF peak to valley following the peak, for the raw signal intensities and for the LUT corrected [Gd].
Simulations and Look-up table calculations
Performance of the 6-pulse saturation recovery preparation is shown in Fig. 5. Over the target design range of ±150 Hz and 17.5–29.6 μT (65–110% of effective FA) (dotted white box), the residual magnetization was < 0.5% in the myocardium for [Gd] up to 1 mmol/L, and for blood was < 0.5% up to 2.5 mmol/L, < 1% up to 5 mmol/L, and < 2.5% up to 10 mmol/L.
LUT’s relating the normalized signal (SR/PD) and [Gd] were calculated for ±20% variation in transmitted B1 from nominal (Fig. 6). The error in [Gd] due to the LUT using the assumed specified FA rather than actual was <10% over this range up to 1 mmol/L and <3% in the blood up to 10 mmol/L.
The T1 weighting of the PD images had < 1% variation in signal intensity over 0-2 mmol/L in the myocardium, and <1% over 0-1 mmol/L in the blood. The [Gd] in blood was < 1 mmol/L for a rest study following stress by several minutes, therefore using a fixed value of native T10 had <1% effect on the LUT. The effect of in-flow on the AIF LUT was calculated to be <4% assuming all spins were refreshed every RF pulse.
Measurement of r1 and r2 relaxivities was made for both Gadobutrol (Gadovist®) and Gadoterate meglumine (Dotarem®) doped saline phantoms. For Gadobutrol, the measured values for r1 and r2 were 5.5 L/mmol/s and 6.8 L/mmol/s, respectively and for Gadoterate meglumine the measured values were 4.6 L/mmol/s and 5.7 L/mmol/s, respectively. The measured [Gd] versus actual [Gd] is shown with and without T2* correction (Fig. 7). The linear fit for Gadobutrol was [Gd]estimate = 0.99 [Gd] + 0.0002 with T2* correction and was [Gd]estimate = 0.90 [Gd] + 0.08, without T2* correction. The linear fit for Gadoterate meglumine was [Gd]estimate = 1.02 [Gd] + 0.07 with T2* correction and was [Gd]estimate = 0.94 [Gd] + 0.12, without T2* correction. Measurements of [Gd] for the myocardial signal protocol are shown in Fig. 8 for TS = 95 ms. For TS = 95 ms the fits were 1.004[Gd] + 0.005 and 1.04[Gd] + 0.01 for SSFP protocol with Gadobutrol and Gadoterate meglumine, respectively, and were 0.96[Gd] + 0.01 and 0.96[Gd] + 0.02 for the FLASH protocol with Gadobutrol and Gadoterate meglumine, respectively. The measurements were made for TS = 65 to 125 ms in steps of 10 ms. For SSFP, the slopes of the fits were within 4% of unity slope for all TS values and both agents and for FLASH were within 5%.
Invivo AIF data
Adenosine stress studies were conducted on 29 normal healthy volunteers. The peak [Gd] was 3.5 ± 1.2 mmol/L (m ± SD) for stress and 4.4 ± 1.2 mmol/L for rest as measured in the AIF. The T2* in the LV blood pool at peak AIF [Gd] was 10.0 ± 2.4 ms at stress and 9.9 ± 1.7 ms at rest. The duration of the 1st pass was 10.3 ± 2.1 s at stress and 14.7 ± 3.2 s at rest. Example AIF images for echo 1 and echo 2 are shown in Fig. 9 for stress, and normalized AIF signal curves are shown in Fig. 10. In this example, the T2* corrected signal was approximately 8% higher than the echo 1 peak signal, which led to approximately 10% greater [Gd] after LUT correction. For 29 subjects, the peak to valley ratio was 5.6 for the raw signal intensities without correction, and was 8.3 for the LUT corrected AIF in gadolinium concentration units. The valley is indicated by the arrow in Fig. 10 stress signal intensity plot. This represents approximately 48% improvement in linearity.
Invivo myocardium data
The SNR of the myocardial time intensity was measured on a pixel-wise basis using SNR scaled reconstruction and measured at peak myocardial enhancement. Peak myocardial SNR was 21.8 ± 7.6 at stress and the peak [Gd] was 0.49 ± 0.15 mmol/L. Example of myocardial stress perfusion images are shown before and after normalization (Fig. 11) and myocardial blood flow maps are shown in Fig. 12. Images are well saturated as observed at baseline. Influence of T2* correction on flow comparing myocardial blood flow estimates with and without T2* correction is shown in Fig. 13. Without T2* correction the myocardial perfusion estimates of blood flow are overestimated by 10%. Estimates of perfusion stress flow using the BTEX model was 3.93 ± 0.38 and rest flow was 1.03 ± 0.19 ml/min/g (N = 29). Estimates for extraction fraction were 0.5 ± 0.04 and 0.85 ± 0.03, at stress and rest, respectively. Estimates of the permeability surface area product (PS) were 1.55 ± 0.2 and 1.33 ± 0.21 (ml/min/g), at stress and rest, respectively. Estimates for the interstitial volume fraction (%) were 27.4 ± 5.9 and 24.8 ± 5.9, at stress and rest, respectively. Estimates for the blood volume fraction (ml/g) were 13.0 ± 0.85 and 9.2 ± 0.76, at stress and rest, respectively. Estimates of perfusion flow using the Fermi model fit over the 1st pass were 3.4 ± 0.39 and 0.95 ± 0.16 ml/min/g, at stress and rest, respectively.
Raw images for a typical case corresponding to the example in Figs. 10 and 12 are provided as supplemental data as movies to include raw AIF images at stress and rest for both echo times before [see Additional file 1] and after respiratory MOCO [see Additional file 2], and the multislice stress and rest myocardium images including raw images [see Additional file 3], MOCO images [see Additional file 4], and MOCO images including surface coil correction [see Additional file 5].
The dual sequence approach was chosen for automated, in-line perfusion mapping since it is readily integrated into a clinical workflow. Unlike the dual bolus method, simultaneous measurement of the AIF and myocardial signals avoids physiological variation between bolus injections such as those due to differences in respiration. Another benefit of the dual sequence is that it decouples the measurement of the AIF from the myocardial imaging protocol so that they may be independently optimized. The inherent non-linear response of SR on the AIF was minimized by design of the protocol and post-processing. Earlier dual sequence AIF protocols used centric ordered acquisition to minimize the TS for improved linearity, but this leads to a high pass spatial filtering of the blood pool signal which becomes dependent on gadolinium concentration  and creates a dependence on the AIF and how the blood pool is segmented, i.e., the edges of the blood pool will have a longer effective saturation delay. Use of a linear ordering leads to a more homogeneous blood pool image. The saturation delay was minimized by use of parallel imaging acceleration to reduce the number of phase encodes lines actually acquired. The AIF signal with linear phase encoder order was slightly more non-linear than centric but could be corrected by a look-up-table approach based on Bloch signal calculations of the normalized signal (SR/PD). In this way, the blood pool signal could be automatically segmented and the AIF could be reliably estimated. Parallel imaging was used to accelerate the acquisition of myocardial perfusion images in order to reduce the single shot duration and thereby mitigate dark rim artifacts to some extent. Gibb’s ringing was suppressed by raw filtering .
It is difficult to make direct comparison between the T2* values reported here and many of the previously reported measurements since many of the previous publications used different Gd contrast agents (e.g., Magnevist), administered different concentrations of [Gd] (0.05-0.1 mmol/kg), and used different infusion rates (3-7 mL/s). Additionally, several publications measured the signal loss with specific sequence parameters and not actual measurement of T2* with multiple echo times. Thus effects of in-flow due to FA, phase encode acquisition order, slice thickness, will vary between sequences. Finally, not all work was done at 1.5 T. These factors make the direct comparison difficult. The paper by de Bazelaire  predicts a value of T2* of 7.5 ms at a peak concentration of [Gd] = 4 mmol at 3 T, and does not give values for 1.5 T which are expected to be longer. By comparison the T2* reported here at 1.5 T for estimated concentration of approx. 4 mmol was approx. 10 ms. Prior work with Magnevist at 1.5 T, 0.1 mmol/kg dose, 5 mmol/s infusion  reported T2* values of approx. 9 ms at peak [Gd] concentration. It is also important to adequately shim the volume to minimize intravoxel de-phasing such that the measured value represents the intrinsic T2*.
The correction of T2* in the AIF avoid underestimating the peak [Gd]. Underestimation the input function [Gd] will result in an overestimation of the flow. This is true in general for all perfusion models (BTEX, Fermi, exponential) that are implemented based on the absolute [Gd] signals.
Conversion to [Gd]
The conversion to gadolinium concentration units facilitated the use of the dual sequence which used a FLASH readout for the AIF and could support either FLASH or SSFP for myocardial perfusion imaging. The relaxivities (r1 and r2) were measured in Gd doped saline phantoms. Although the estimated concentration is dependent on the accuracy of these values, which may be slightly different in blood plasma, the quantified myocardial blood flow has been found through simulation to be quite insensitive since they affect the scale of both blood and myocardium. In-vivo values of blood T2* were significantly lower than for saline phantoms at a given [Gd] as previously reported . The average blood T2* at 1.5 T was approximately 10 ms at 4 mmol/L peak [Gd], whereas the T2* was approximately 25 ms in saline at the same concentration. The T2* decreases with field strength  and values as low as 4 ms are likely to be encountered at 3 T using a dose of 0.05 mmol/kg leading to higher signal losses and greater importance for T2* correction. The conversion to [Gd] in the presence of error in FA due to unknown B1 was analyzed through simulation to cause small errors in [Gd] which in turn will result in an error in estimated flow of < 15% for ±20% variation in B1.
The in-line automated perfusion mapping software has the capability of calculating perfusion estimates based on several different widely used models with differing complexity including the Fermi model  and distributed BTEX model . The distributed BTEX model explicitly estimates the permeability surface area (PS) product which is used to calculate an extraction fraction to account for the flow dependent leakage of Gd from the vascular space into the interstitium. The BTEX model has been used in 13N-Ammonia PET and validated against microspheres . There is growing interest in distributed tissue models for perfusion [28–30]. A more comprehensive comparison of models will be the subject of a detailed study. Despite significant differences in models, the mean values for perfusion estimates compare reasonably well with other reported values in normal subjects. In a study by Broadbent, et al. , perfusion was estimated using both Fermi and a distributed parameter (DP) model. Values for stress and rest perfusion estimates and myocardial flow reserve in that study were 3.8/1.5/2.7 for Fermi and 3.5/1.5/2.5 for DP, as compared to 3.4/0.95/3.6 for Fermi and 3.9/1.03/3.8 for BTEX in the present study. In a study of 10 normal subjects, Hsu, et al.  reported values of 3.39/1.02/3.3 using a Fermi tissue model.
The 6-pulse saturation preparation achieved > 99% saturation over a wide range of off-resonance, effective transmitter FA, and gadolinium concentration. Excellent saturation mitigates slice to slice cross-talk and allows the user to prescribe a mixture of short and long axis slices. There is possibility of slight cross talk between intersecting slices for the initial non-saturated PD frames, but this is quite small due to the low FA of PD weighted acquisition protocol. The 6-pulse is designed conservatively, and has recently been reduced to a 5-pulse design to reduce the overall duration and SAR. The newly design preparation has recently been evaluated and perform with >99% saturation over a design range of 50–110% nominal B1+. This newly designed SR preparation is 14.7 ms including crushers which translated to a total AIF duration of 57 ms. Using this 5-pulse design allows for slight increase of the saturation delay from 95 to 105 ms, thereby increasing the SNR while maintaining 3-slices up to 120 bpm.
A dual sequence for myocardial perfusion CMR and arterial input function measurement has been optimized for quantification of myocardial blood flow. A validation in phantoms was performed to confirm that the signal conversion to gadolinium concentration was linear. The proposed sequence was integrated with a fully automatic in-line solution for pixel-wise mapping of myocardial blood flow and evaluated in adenosine stress and rest studies on N = 29 normal healthy subjects. Reliable perfusion mapping was demonstrated and produced estimates with low variability.
Arterial input function
Cardiovascular magnetic resonance
Dark rim artifact
Fast low angle shot
Late gadolinium enhancement
Specific absorption ratio
Steady state free precession
Wilke N, Simm C, Zhang J, Ellermann J, Ya X, Merkle H, et al. Contrast-enhanced first pass myocardial perfusion imaging: correlation between myocardial blood flow in dogs at rest and during hyperemia. Magn Reson Med. 1993;29:485–97.
Jerosch-Herold M, Wilke N, Stillman AE, Wilson RF. Magnetic resonance quantification of the myocardial perfusion reserve with a Fermi function model for constrained deconvolution. Med Phys. 1998;25:73.
Christian TF, Aletras AH, Arai AE. Estimation of absolute myocardial blood flow during first-pass MR perfusion imaging using a dual-bolus injection technique: comparison to single-bolus injection method. J Magn Reson Imaging. 2008;27:1271–7.
Gatehouse PD, Elkington AG, Ablitt NA, Yang G-Z, Pennell DJ, Firmin DN. Accurate assessment of the arterial input function during high-dose myocardial perfusion cardiovascular magnetic resonance. J Magn Reson Imaging. 2004;20:39–45.
Xue H, Hansen MS, Nielles-vallespin S, Arai AE, Kellman P. Inline quantitative myocardial perfusion flow mapping. JCMR/ISMRM Work. 2016;18:4–6.
Chow K, Kellman P, Spottiswoode BS, Nielles-Vallespin S, Arai AE, Salerno M, et al. Saturation pulse design for quantitative myocardial T1 mapping. J Cardiovasc Magn Reson. 2015;17:84.
Nielles-Vallespin S, Kellman P, Hsu L-Y, Arai AE. FLASH proton density imaging for improved surface coil intensity correction in quantitative and semi-quantitative SSFP perfusion cardiovascular magnetic resonance. J Cardiovasc Magn Reson. 2015;17:16.
Cernicanu A, Axel L. Theory-based signal calibration with single-point T1 measurements for first-pass quantitative perfusion MRI studies. Acad Radiol. 2006;13:686–93.
Sekihara K. Steady-state magnetizations in angles and short repetition intervals. IEEE Trans Med Imaging. 1987;MI-6:157–64.
Gatehouse P, Lyne J, Smith G, Pennell D, Firmin D. T2* effects in the dual-sequence method for high-dose first-pass myocardial perfusion. J Magn Reson Imaging. 2006;24:1168–71.
Sánchez-González J, Fernandez-Jiménez R, Nothnagel ND, López-Martín G, Fuster V, Ibañez B. Optimization of dual-saturation single bolus acquisition for quantitative cardiac perfusion and myocardial blood flow maps. J Cardiovasc Magn Reson. 2015;17:21.
de Bazelaire C, Rofsky NM, Duhamel G, Zhang J, Michaelson MD, George D, et al. Combined T2* and T1 measurements for improved perfusion and permeability studies in high field using dynamic contrast enhancement. Eur Radiol. 2006;16:2083–91.
Kellman P, McVeigh ER. Image reconstruction in SNR units: a general method for SNR measurement. Magn Reson Med. 2005;54:1439–47.
Xue H, Zuehlsdorff S, Kellman P, Arai A, Nielles-Vallespin S, Chefdhotel C, et al. Unsupervised inline analysis of cardiac perfusion MRI. Lect. Notes Comput. Sci. (including Subser. Lect. Notes Artif. Intell. Lect. Notes Bioinformatics). 2009;5762 LNCS:741–9.
Breuer FA, Kellman P, Griswold MA, Jakob PM. Dynamic autocalibrated parallel imaging using temporal GRAPPA (TGRAPPA). Magn Reson Med. 2005;53:981–5.
Kim D. Influence of the k-space trajectory on the dynamic T1-weighted signal in quantitative first-pass cardiac perfusion MRI at 3 T. Magn Reson Med. 2008;59:202–8.
Elkington AG, He T, Gatehouse PD, Prasad SK, Firmin DN, Pennell DJ. Optimization of the arterial input function for myocardial perfusion cardiovascular magnetic resonance. J Magn Reson Imaging. 2005;21:354–9.
Zhao L, Salerno M, Kramer CM, Meyer CH. The contribution of cardiac motion to dark rim artifacts in myocardial perfusion scans. Proc ISMRM. 2010;18:3626.
Storey P, Chen Q, Li W, Edelman RR, Prasad PV. Band artifacts due to bulk motion. Magn Reson Med. 2002;48:1028–36.
Di Bella EVR, Parker DL, Sinusas AJ. On the dark rim artifact in dynamic contrast-enhanced MRI myocardial perfusion studies. Magn Reson Med. 2005;54:1295–9.
Hansen MS, Sørensen TS. Gadgetron: an open source framework for medical image reconstruction. Magn Reson Med. 2013;69:1768–76.
Chefd’Hotel C, Hermosillo G, Faugeras O. Flows of diffeomorphisms for multimodal image registration. IEEE Int Symp Biomed Imaging. 2002;753–756.
Bassingthwaighte JB, Wang CY, Chan IS. Blood-tissue exchange via transport and transformation by capillary endothelial cells. Circ Res. 1989;65:997–1020.
Zhou Z, Bi X, Wei J, Yang H-J, Dharmakumar R, Arsanjani R, et al. First-pass myocardial perfusion MRI with reduced subendocardial dark-rim artifact using optimized Cartesian sampling. J Magn Reson Imaging. 2017;45(2):542–55. doi:10.1002/jmri.25400.
Kellman P, Aletras AH, Hsu L-Y, McVeigh ER, Arai AE. T2* measurement during first-pass contrast-enhanced cardiac perfusion imaging. Magn Reson Med. 2006;56:1132–4.
Kalavagunta C, Michaeli S, Metzger GJ. In vitro Gd-DTPA relaxometry studies in oxygenated venous human blood and aqueous solution at 3 and 7 T. Contrast Media Mol Imaging. 2014;9:169–76.
Alessio A, Bassingthwaighte JB, Glenny R, Caldwell J. Validation of an axially distributed model for quantification of myocardial blood flow using 13 N-ammonia PET. J Nucl Cardiol. 2013;20:64–75.
Broadbent DA, Biglands JD, Larghat A, Sourbron SP, Radjenovic A, Greenwood JP, et al. Myocardial blood flow at rest and stress measured with dynamic contrast-enhanced MRI: comparison of a distributed parameter model with a Fermi function model. Magn Reson Med. 2013;70:1591–7.
Chung S, Shah B, Storey P, Iqbal S, Slater J, Axel L. Quantitative Perfusion Analysis of First-Pass Contrast Enhancement Kinetics: Application to MRI of Myocardial Perfusion in Coronary Artery Disease. PLoS One [Internet]. 2016;11:e0162067. Available from: http://dx.plos.org/10.1371/journal.pone.0162067
Kunze KP, Rischpler C, Hayes C, Ibrahim T, Laugwitz K-L, Haase A, et al. Measurement of extracellular volume and transit time heterogeneity using contrast-enhanced myocardial perfusion MRI in patients after acute myocardial infarction. Magn Reson Med. 2016. doi:10.1002/mrm.26320. [Epub ahead of print]
Hsu L-Y, Rhoads KL, Holly JE, Kellman P, Aletras AH, Arai AE. Quantitative myocardial perfusion analysis with a dual-bolus contrast-enhanced first-pass MRI technique in humans. J Magn Reson Imaging. 2006;23:315–22.
We acknowledge Kelvin Chow for his design of the saturation recovery preparation.
Supported by the National Heart, Lung and Blood Institute, National Institutes of Health by the Division of Intramural Research.
Availability of data and materials
The raw data that support the findings of this study are available from the corresponding author upon reasonable request subject to restriction on use by the Office of Human Subjects Research. Raw data require reconstruction processing.
PK and HX conceived of the study, developed the sequence and Bloch simulations, performed processing and analysis, and drafted the manuscript. SNV contributed to sequence development, and MSH contributed to the image reconstruction framework and data analysis. JN, RT, and MU acquired normal volunteer data. All authors participated in revising the manuscript and read and approved the final manuscript.
The authors declare that they have no competing interests
Consent for publication
Written informed consent was obtained from patients for publication of their individual details and accompanying images in this manuscript. The consent form is held in the patients’ clinical notes and is available for review by the Editor-in-Chief.
Ethical approval and consent to participate
The study was approved by Regional Human Subjects Research Ethics Committee in Stockholm, Sweden, and all subjects provided written informed consent. Imaging was performed at the Karolinska University Hospital, Stockholm, Sweden. Anonymized data was analysed at NIH with approval by the NIH Office of Human Subjects Research OHSR (Exemption #13156).
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Kellman, P., Hansen, M.S., Nielles-Vallespin, S. et al. Myocardial perfusion cardiovascular magnetic resonance: optimized dual sequence and reconstruction for quantification. J Cardiovasc Magn Reson 19, 43 (2017) doi:10.1186/s12968-017-0355-5
- Myocardial perfusion
- Perfusion quantification
- Arterial input function