Introduction

Cardiac magnetic resonance imaging (cMRI) has become an increasingly important imaging technique, which enables various non-invasive diagnostic options in clinical cardiology and basic cardiologic research. In addition to informative morphological and functional investigations, cMRI offers the possibility of enhanced tissue characterization to investigate tissue defects like edema, inflammation or fibrosis [1,2,3]. There are several methods facilitating magnetic resonance imaging to enable tissue characterization, such as late gadolinium enhancement (LGE), T1- and T2 mapping, and the quantification of the extracellular volume [4,5,6,7,8]. However, some of these methods require the administration of gadolinium contrast agent, which is contraindicated in certain patients, especially those with reduced renal function.

Native T1ρ quantification has been emerged to be a promising alternative. The T1ρ relaxation mechanism, also called spin lattice relaxation in the rotating frame, uses an on-resonant radiofrequency pulse that locks the magnetization and prohibits free relaxation of the spin ensemble [9,10,11]. This applied spin lock (SL) pulse causes a high sensitivity to low frequency processes at the molecular and cellular level (e.g. of biologic macromolecules) [12,13,14]. In the field of cMRI, the SL method, using a moderate locking amplitude, has been shown to suppress low frequency relaxation mechanisms that obscure endogenous contrast [15]. In chronic infarcts in the swine model, T1ρ revealed a significant increase compared to healthy remote myocardium even without the use of contrast agents [15]. Compared to other native quantification techniques, T1ρ improves the contrast ratio between diseased and healthy myocardial tissue [16,17,18]. In experimental models as well as in vivo, T1ρ has already been shown to be a sensitive marker for the detection of several tissue damages like edema, ischemia, fibrosis and myocardial infarction [19,20,21,22]. Moreover, T1ρ gives the possibility for dispersion measurements by manipulation of the effective SL pulse amplitude, enabling additional contrast mechanisms and further improved tissue characterization [11]. In this context it could already be shown that T1ρ dispersion quantification can be used as a myocardial biomarker for the classification of different stages of fibrosis [19]. This dispersion behavior is of great interest and might be highly beneficial not only in basic research in various animal models, but also in various scenarios in clinical practice [23,24,25]. However, fast and accurate T1ρ dispersion quantification in cMRI is still a major imaging challenge, since the preparation of T1ρ must be carried out immediately before each acquisition [26, 27] and the acquisition time is strongly limited by physiological parameters such as the breathing cycle and heart rate. For the calculation of a single T1ρ map, several images with different T1ρ weightings must be acquired. To quantify the dispersion, several maps must be measured using different SL amplitudes. This procedure can lead to excessive measurement times, especially in small animal studies.

Despite the many promising possibilities of cardiac T1ρ dispersion quantification, only minor attention has been paid to the technique itself so far. To date, there are only a few publications that deal with myocardial T1ρ dispersion analysis. In the work presented by Witschey et al. [15] myocardial T1ρ mapping was performed in vivo using an infarction model in pigs. Here, the infarction area could be clearly visualized using cardiac T1ρ mapping. However, T1ρ dispersion measurements were only performed ex vivo, revealing a significant difference in the dispersion behavior between infarcted, border zone and remote myocardium. Musthafa et al. [20] presented a study in a mouse infarction model at which T1ρ mapping was performed at several time points after infarction, accounting a significant increase of T1ρ at day 7 after infarction. The sequence used for high-resolution T1ρ mapping was based on respiratory and heartbeat triggered spin locking followed by a single Cartesian spin echo readout. This acquisition procedure ensured a high signal-to-noise ratio (SNR) but was relatively slow, requiring a total scan time of approximately 20 min for a single-slice T1ρ map, although only four different T1ρ weighted images were acquired for mapping. This work also introduced an accelerated proof of concept T1ρ dispersion measurement with reduced spatial resolution. However, these measurements only seem to roughly estimate the dispersion behavior and do not allow an exact pixel-by-pixel quantification. In the work published by Yin et al. [19] T1ρ dispersion analysis has been performed in a fibrosis model in several dogs in vivo. The authors introduced a dispersion dependent myocardial biomarker that enables the distinction between different grades of fibrosis. A drawback of the study is that only two relaxation maps were used for dispersion quantification and one of these maps did not represent a true T1ρ map but rather a T2 map. Over all, the main reason why there is no suitable myocardial T1ρ dispersion quantification method so far is the fact that an accurate measurement can take several hours.

The aim of our work is to introduce a very fast and accurate T1ρ quantification technique that enables high-resolution myocardial T1ρ mapping and T1ρ dispersion quantification even within the limited time window of a small animal in vivo study. For this, a Bloch simulation-optimized imaging sequence using high flip angles and a radial view-sharing method has been developed ensuring high SNR and efficient data sampling. In our study, the new radial sampling technique was compared in phantom measurements with a conventional radial acquisition scheme and a reference gold standard technique. In addition, measurements were carried out in N = 10 healthy mice to test the in vivo applicability and to carry out, for the first time, quantitative dispersion imaging in the myocardium.

Methods

Sequence design

The optimized T1ρ mapping sequence has been designed to acquire a series of T1ρ weighted images with different SL times tSL within a single measurement. Therefore, a signal intensity dependent radial sampling pattern and a k-space view-sharing method has been used. T1ρ preparation was performed by a balanced spin locking module, which includes two adiabatic half-passage (AHP) excitation pulses, three continuous wave SL pulses using alternating phases and two opposite 180° refocusing pulses for improved B0 and B1 insensitivity [28, 29]. For data acquisition, a golden angle radial gradient echo readout was used [30], acquiring four radial spokes after each SL preparation (Fig. 1). The acquisition window was positioned in end diastole using a dynamic trigger delay (depending on tSL) after the R-wave of the ECG signal. Each preparation experiment was separated by a waiting time trec for magnetization recovery, which is dependent on the respiratory cycle rate. Data sampling was segmented into 13 preparation experiments with identical preparation characteristics (identical SL time tSL and SL amplitude fSL). These 13 identical preparation experiments were carried out with eight different SL times, enabling the reconstruction of eight images with different T1ρ weightings using a k-space weighted image contrast (KWIC) filtered view-sharing method [31]. This leads to 13 × 8 = 104 consecutive preparation experiments and consequently 104 × 4 = 416 radial acquisitions in total for the generation of a single T1ρ map. The acquisition time at a respiratory rate of 1 Hz is ≈1.7 min.

Fig. 1
figure 1

Sequence design for myocardial T1ρ mapping in mice. After prospective respiratory gating a trigger on the R-wave has been used followed by a dynamic trigger delay and a balanced SL preparation module [28]. Thereafter, four radial spokes of the k-space were acquired in end diastole during one heartbeat. This procedure was repeated in each respiratory cycle separated by a waiting time trec for magnetization recovery. For each T1ρ weighting, 13 repetitions were used

Full size image

Concept of Bloch sorting

The fundamental idea of our fast T1ρ mapping technique is based on three key points. First, the increase of SNR using high flip angles in the gradient echo readouts. Second, the algorithmic search for an optimal radial sampling pattern based on Bloch simulations. Third, the image reconstruction of undersampled data using a KWIC-filtered view-sharing method.

Due to the very high heart rate in mice (≈ 450 bpm), only a few readouts after the spin lock preparation are possible in the acquisition window of the diastole (≈ 20 ms). Hence, short repetition times TR (≈ 5 ms) and only four readouts per preparation were used. The waiting time for magnetization recovery is set by the respiratory cycle (breath gating, trec ≈ 1500 ms) and is therefore relatively long. Since this means that only a small number of acquisitions are carried out over time, the signal level must be maximized. In this specific case, the maximum of the averaged signal intensity \(\overline{S }\) is achieved if appropriately high flip angles are used in the RF pulse train. Without taking relaxation effects into account, the optimal flip angle can be calculated as follows:

$${\alpha }_{\mathrm{opt}}=\underset{\alpha }{\mathrm{arg max}} \left[\overline{S }\left(\alpha \right)\right]=\underset{\alpha }{\mathrm{arg max}}\left[{\sum }_{k=1}^{4}\mathrm{sin}\left(\alpha \right) {\mathrm{cos}\left(\alpha \right)}^{ k-1}\right]={43.51}^{^\circ }$$
(1)

By considering typical relaxation times of myocardial tissue (T1ρ = 40 ms [20], T1 = 1400 ms [32]) and the sequence timings (TR = 5 ms, trec = 1500 ms) described above, the signal maximum is reached at a slightly different value. In Fig. 2a the signal of the NR = 4 readouts was simulated for different SL times (tSL = 4, 12, 20, 28, 36, 44, 52, 60 ms) by solving the Bloch equations [33]. In Fig. 2b, the simulation was carried out for different flip angles and plotted averaged over all SL times and readouts. The signal maximum is reached at αopt = 39.35°. Based on this estimate, we consistently used α = 40° for all measurements and simulations in this study. Since the T1 recovery is not complete between sequence repetitions, we carried out two dummy cycles prior to the first preparation experiment to prevent increased signals in the first acquisitions. The high flip angles are advantageous here because the steady-state is reached quickly.

Fig. 2
figure 2

Optimization of the flip angle for the radial gradient echo readout. a Results of the Bloch simulation for two exemplary flip angles. At low flip angles, the signal intensities of the different readouts differ only slightly. At higher flip angles, a variation in the signal intensities is observed. b The plot shows the simulated signal averaged over all tSL for different flip angles. For the sequence parameters TR = 5 ms, trec = 1500 ms, tSL = 4…60 ms and the relaxation times T1 = 1400 ms and T1ρ = 40 ms (typical for myocardial tissue), the mean signal level \(\overline{S }\) is maximized for αopt = 39.35°. This value strongly depends on the number of readouts after preparation (NR = 4)

Full size image

The signal maximization generated using high flip angles for the readout (α = 40°) has the disadvantage of strong signal variations after each T1ρ preparation (Fig. 3a). This would ultimately lead to incorrect T1ρ weighting if all projections were used equally for the sampling of k-space. A fundamental idea of our method is therefore to generate a smooth variation of signal intensities across k-space using an optimized sorting of golden angles. For this a Bloch simulation-optimized sampling scheme was developed. The signal intensity for every acquisition window was roughly predicted prior to the measurement using the known sequence parameters (TR, trec and α) and estimated T1 and T1ρ values of the probe under investigation by solving the Bloch equations [33]. The results were used to create a sorting (Bloch sorting) of the predicted signal levels for the corresponding acquisition windows (Fig. 3b). The projection angle \(\phi\) was calculated for each radial readout by linking neighboring golden angles (111.25°, 222.50°, 333.75°, …) with the corresponding Bloch sorting position \({N}_{\mathrm{BS}}\) (Fig. 3c).

Fig. 3
figure 3

Concept of Bloch sorting (a–c) and the image reconstruction technique (d). a The expected signal intensity is plotted versus the chronological readout number of the acquired spokes. For Bloch simulation the values T1 = 1400 ms, T1ρ = 60 ms, TR = 5 ms, trec = 1000 ms, α = 40° and tSL = 4…102 ms (linear spacing) were used in the phantom experiments. The colors of the plateaus represent the different acquisition windows after SL preparation (see Figs. 1 and 2). The roman numerals indicate the sorting position of the plateaus. The numbers in gray represent the indices of the golden angles numbers. b Signal intensity of subsequent golden angles after Bloch sorting. In the ideal case, the curve is monotonically decreasing. The change in contrast between subsequent spokes is minimized. c Bloch simulation-optimized sampling pattern. A projection angle ϕ is calculated for each readout by linking subsequent golden angles to the sorting position. d KWIC profile for the reconstruction of a T1ρ weighted image and the corresponding k-space sampling trajectory. The 13 readouts directly after T1ρ preparation are used exclusively to sample the corresponding k-space center. The k-space periphery is also filled with other weightings

Full size image
$$\phi ={N}_{\mathrm{BS}}\cdot \frac{360^\circ }{1+\sqrt{5}}\approx {N}_{\mathrm{BS}}\cdot 111.25^\circ$$
(2)

The image reconstruction was performed using an adapted KWIC-filtered technique [31] (Fig. 3d). Here, the k-space center for a desired T1ρ weighting (with desired tSL) is exclusively selected from the first readout after preparation. This ensures that the main image contrast of the specific T1ρ weighted image is provided by the correct SL time. Since we have 13 (Fibonacci number) radial projections available for each preparation experiment, the sampling density in the azimuthal direction is homogeneous and only two different azimuthal gaps exist (Fig. 3d) [30]. Due to the Bloch sorting technique, the change in contrast between neighboring golden angles used for the k-space filling (based on the KWIC filter) is minimized for the k-space center. The k-space peripherals were also chosen from other acquisition windows as well as other T1ρ weightings. Here, the k-space was filled in segments according to the Bloch sorting principle. In each segment, the total number of projections \({n}_{\phi }\) agrees to a Fibonacci number (13, 21, 34, 55, …) to ensure a homogenous sampling density. We further ensured that the Nyquist criterion was met in each region of the segments (detailed explanation in the supplementary material, Online Fig. 1). For \({n}_{\phi }\) projections in a segment, the number of samples in radial direction \({n}_{r}\) (kr direction, Fig. 3d) was linked to the condition \({n}_{\phi }>\pi \cdot { n}_{\mathrm{r}}\). The reconstruction of the radial k-space data was done using a nonuniform fast Fourier transform (NUFFT) with an open source software toolbox [34, 35].

Phantom measurements

All measurements were performed on a 7.0 T small animal imaging system (Bruker BioSpec 70/30, Bruker BioSpin MRI GmbH, Ettlingen, Germany) with a maximum gradient field strength of 470 mT/m. A 35 mm homebuilt quadrature transmit-receive birdcage was used for signal detection. The phantom used consisted of four cylindrical sample tubes with a diameter of 17 mm and a length of 120 mm. The tubes were filled with different concentrations (10, 15, 20, and 25%) of BSA (Bovine Serum Albumin, Sigma-Aldrich, St. Louis, MO, USA) resulting in different T1ρ values in the typical range of biological tissue. The sample tubes were arranged in a quadratic array, which was placed in the isocenter of the magnet. In the phantom measurements, ECG triggering of the sequence was deactivated and the recovery time was fixed to a constant value trec = 5000 ms.

To demonstrate the advantages of the Bloch simulation-based sampling scheme, the image quality and the T1ρ quantification accuracy were compared with a conventional sampling scheme using a serial sorting of golden angles for subsequent readouts (Fig. 3c). In this setup, the projection angles were not optimized for the expected signal intensity. The remaining sequence parameters for the acquired T1ρ weighted images were adjusted identical for both sampling schemes: field of view (FOV) = 38.4 × 38.4 mm2, slice thickness = 1.5 mm, acquired/reconstructed resolution = 128 × 128 pixels, repetition time (TR) = 5.0 ms, echo time (TE) = 2.0 ms, bandwidth = 75 kHz, flip angle α = 40°, acquired spokes after SL preparation = 4, tSL = 4…102 ms (8 different, linear spacing), fSL = 1500 Hz. The calculation of a single T1ρ map required 104 (13 × 8) SL preparation experiments.

An artifact/SNR analysis was performed to determine the sampling scheme that achieves the best image quality with least artifacts, due to the undersampling of the k-space and the KWIC filter based image reconstruction. Therefore, the SNR for every T1ρ weighted image has been calculated and averaged over all reconstructed images. The individual SNR values were calculated from the magnitude images using four signal masks and a noise mask (Fig. 5) avoiding the edges of the phantoms. The ratio was built from the mean value within the signal masks and the standard deviation within the noise mask. In addition, the coefficient of determination R2, which represents a measure of the agreement with the mono-exponential model function, has been calculated pixel-wise for each fitting process of the T1ρ exponential decay and was finally averaged to obtain a global R2 indicator for both sampling schemes.

For validation of the T1ρ quantification accuracy pixel-wise comparisons with a Cartesian SL prepared turbo spin echo (TSE) sequence were performed, which serves as the gold standard reference. The TSE sequence parameters were chosen similar to the radial sequence and the identical balanced spin lock preparation was used: FOV = 38.4 × 38.4 mm2, slice thickness = 1.5 mm, resolution = 128 × 128 pixels, TR = 5031.4 ms (including trec), TE = 7.0 ms, bandwidth = 50 kHz, turbo factor = 4, trec = 5000 ms, tSL = 4…102 ms (eight different, linear spacing), fSL = 1500 Hz. T1ρ mapping using the reference TSE sequence required 256 (32 × 8) SL preparation experiments. Hence, our new radial T1ρ mapping sequence is about 2.5 times (32:13) faster than the TSE reference measurement.

Furthermore, the accuracy of the T1ρ dispersion quantification was examined by comparing the novel radial Bloch sorting technique and the TSE reference. Here, T1ρ mapping was performed for eight different SL amplitudes (fSL = 750…2500 Hz, linear spacing). Thus 8 × 8 × 13 = 832 preparation experiments were required for the radial acquisitions and 8 × 8 × 32 = 2048 preparations for the TSE sequence. For data analysis, circular regions of interest (ROIs) were drawn at the positions of the four T1ρ sample probes. In these ROIs, the T1ρ mapping results were compared pixel-wise with the corresponding values obtained from the TSE reference sequence. The mean quantification accuracy and its variance were determined based on the eight mapping experiments. In addition, the T1ρ dispersion was analyzed using a linear dispersion model to the eight acquired T1ρ maps with various SL amplitudes fSL.

$${T}_{1\uprho }\left({f}_{\mathrm{SL}}\right)={T}_{1\uprho }^{0}+{m}_{1\uprho }\cdot {f}_{\mathrm{SL}}$$
(3)

The quantification accuracy of the dispersion offset \({{T}_{1\uprho }\left({f}_{\mathrm{SL}}=0\right)=T}_{1\uprho }^{0}\) and the dispersion slope \({m}_{1\uprho }\) was also investigated by a pixel-wise examination for each sample tube, identifying its deviation from the corresponding pixel in the TSE reference.

In vivo measurements

To check the applicability of our new method in vivo, measurements were carried out in healthy mice. Here the Bloch sorting scheme was used to quantify the T1ρ relaxation time and T1ρ dispersion, since—with respect of the results of the phantom measurements—the Serial sorting scheme did not allow a reasonable application in vivo.

The mice (Naval Medical Research Institute, Charles River Laboratories MA) were imaged in prone position. The animals were anesthetized with isoflurane inhalation (1.5–2 Vol. % in oxygen) and were kept at a constant body temperature of 37 °C. For signal detection the same 35 mm quadrature transmit-receive birdcage as for the phantom measurements was used. Two ECG electrodes attached to the forepaws of the mice were used for ECG triggering and a pressure sensitive balloon placed on the abdominal wall for breath gating. All experimental procedures were in accordance with institutional guidelines and were approved by local authorities.

Myocardial T1ρ mapping was performed in N = 10 mice in a single slice for fSL = 1500 Hz. After a standard planning procedure using high-resolution cine images, a midventricular short-axis imaging slice has been selected for the T1ρ measurements. The acquisition window was positioned in end diastole using a variable trigger delay dependent on tSL. The recovery time trec was chosen dependent on the respiratory cycle, which was ≈ 1500 ms. The further sequence parameters were adjusted similar to the phantom measurements: FOV = 32 × 32 mm2, slice thickness = 1.5 mm, TR = 4.7 ms, TE = 1.9 ms, bandwidth = 75 kHz, α = 40°, acquired spokes after SL preparation = 4, acquired/reconstructed resolution = 128 × 128 pixels, tSL = 4…60 ms (eight different, linear spacing). The in vivo measurement time for each T1ρ map was approximately 2.5 min.

For analysis of the best case reproducibility of the radial T1ρ mapping sequence, a series of ten T1ρ maps with identical SL amplitudes was acquired in one animal in direct succession. A global left ventricular ROI has been selected and was then copied to every repeated T1ρ map. Subsequently, the mean T1ρ values as well as the standard deviations have been calculated at the position of the ROIs. In addition, a comparative measurement using an equivalent Cartesian gradient echo readout was carried out in one animal. Here, fully sampled data sets were acquired for all T1ρ weighted images, leading to a total measurement time of ≈ 6.2 min for the matrix size of 128 × 128. All T1ρ maps obtained were analyzed in a left ventricular segmentation model according to the American Heart Association (AHA).

Furthermore, a detailed analysis of the T1ρ dispersion was carried out in one animal. Therefore, eight T1ρ maps with different SL amplitudes fSL = 750…2500 Hz (eight different, linear spacing, 8 × 8 × 13 = 832 preparation experiments) were acquired. With the acquired data, eight myocardial T1ρ maps and a quantitative T1ρ dispersion slope map were calculated (Eq. 3). Finally, the dispersion behavior of myocardial tissue and the left ventricular blood pool was quantitatively analyzed. Here the measurement time for the acquisition of the complete T1ρ dispersion map was approximately 20 min.

Results

Phantom measurements

The results in Fig. 4 show the measured signal intensities over the projection angle of the respective readouts. Here, a good agreement with the predictions of the Bloch simulation (Fig. 3) can be observed. In the case of Bloch sorting, the intensities were sorted almost perfectly in a descending manner (98%). Only a few readouts were acquired in incorrect order. In the case of Serial sorting fast signal changes were observed due to the high flip angles used.

Fig. 4
figure 4

Comparison of the measured signal intensities for Serial sorting and Bloch sorting in the phantom experiment. a The k-space projections (absolute values) are visualized as heatmaps in order of the corresponding projection angles (golden angle number \({N}_{\upphi }\)). b The signal intensities of the measured k-space centers are plotted versus the golden angle number \({N}_{\upphi }\). In the case of Bloch sorting the change of contrast for subsequent angles was minimized. The mean deviation from ideal sorting (monotonically decreasing) was 2%. For Serial sorting fast signal changes were observed

Full size image

Figure 5a shows the results of the artifact/SNR analysis. The images reveal clearly less artifacts using the optimized Bloch sampling scheme compared to the standard Serial sorting method, which is a result of the reduced contrast differences between acquired subsequent golden angles. The SNR values averaged over all reconstructed T1ρ weighted images were SNRSerial = 17.9 and SNRBloch = 39.3. This represents a SNR increase of 120% using the Bloch sorting scheme. The R2 map (coefficient of determination, Fig. 5b) also shows a better agreement with the fitted exponential T1ρ decay, since the R2 values almost perfectly fit the value 1 using Bloch sorting. In the R2 map of the Serial sorting scheme, certain structural shapes are visible, which might be caused by the artifacts within the reconstructed T1ρ weighted images. Hence, the R2 maps also represents a measure of potential artifacts. The mean R2 values were determined to be \({R}_{\mathrm{Serial}}^{2}=0.973\) and \({R}_{\mathrm{Bloch}}^{2}=0.999\).

Fig. 5
figure 5

Artifact and SNR analysis of the novel Bloch simulation-based sampling scheme and the standard Serial sorting scheme. a Comparison of exemplary T1ρ-weighted images of the four BSA phantoms. Left column: the images are “normally” scaled (all images were acquired with an identical tSL and fSL combination). Right column: the scaling of the images is adapted to highlight artifacts. It is apparent, that Bloch sorting produces much less artifacts than the Serial sorting technique. The SNR values have been determined to SNRSerial = 16.7 and SNRBloch = 44.6 for the images shown. b Comparison of R2 maps. The R2 values of the Bloch sorting method almost perfectly fit the value 1, while the other method only reaches lower values (\({R}_{\mathrm{Serial}}^{2}=0.973\) and \({R}_{\mathrm{Bloch}}^{2}=0.999\)). Structural shapes within the R2 maps are an indication of artifacts within the reconstructed T1ρ weighted images and a measure of reduced image quality

Full size image

The resulting T1ρ maps for both sampling schemes are shown in Fig. 6a. An improved image quality for Bloch sorting is also visible here, since fewer streaking artifacts can be seen in the phantoms. The results of the T1ρ quantification accuracy and precision measurements are illustrated in Fig. 6b. The quantification errors compared to the TSE reference have been determined for the individual BSA tubes using a pixel-wise evaluation (Table 1). The relative error averaged over all phantoms was − 4.01 ± 5.57% (mean value ± standard deviation) for Serial sorting and − 0.17 ± 2.79% for Bloch sorting. The optimized Bloch sorting scheme achieves a considerable improvement in accuracy (+ 56%) and precision (+ 49%). However, for the individual tubes, Bloch sorting showed a tendency to slightly underestimate long T1ρ times and overestimate short times (Table 1).

Fig. 6
figure 6

T1ρ quantification of Serial Sorting and Bloch Sorting. a Calculated T1ρ maps of the four BSA phantoms. Streaking artifacts were clearly reduced with Bloch Sorting. b Comparison of quantification errors using the TSE reference method. For the individual BSA tubes, the relative error in relation to the TSE reference was pixel-wise evaluated. The accuracy and precision has been improved by Bloch sorting for all phantoms. However, a slight underestimation (low BSA concentrations) and overestimation (high BSA concentrations) was obtained

Full size image
Table 1 T1ρ quantification of Serial sorting and Bloch sorting
Full size table

The results of the T1ρ dispersion investigations are shown in Fig. 7a. The mean T1ρ values of the four sample probes are plotted versus the SL amplitude for our radial Bloch sorting scheme and the reference values obtained from the TSE measurement (Fig. 7b). The quantification error was averaged for all SL amplitudes and determined to be − 0.46 ± 1.84%. In Fig. 7c the calculated dispersion slope maps and dispersion offset maps using the linear dispersion model (Eq. 3) are shown. The obtained dispersion slopes \({m}_{1\rho }\) and offsets \({T}_{1\uprho }^{0}\) for the different sample probes and the comparison with the reference TSE values are listed in Table 2. The dispersion behavior of the four samples can be clearly differentiated using our new radial sequence. For the dispersion slope a mean deviation of − 1.3% and a maximum deviation of − 2.8% has been observed. Mean deviation of the dispersion offset was − 0.2% and the maximum deviation was − 3.4%. For the dispersion analysis, good agreement with a linear model could be determined in the range of the SL amplitudes used. The coefficient of determination of the linear regression averaged at least 0.993.

Fig. 7
figure 7

Comparison of the radial Bloch sorting method with the Cartesian TSE reference gold standard. a Results of T1ρ mapping using different SL amplitudes (fSL = 750…2500 Hz). Both methods clearly show T1ρ dispersion for all BSA probes. b Mean T1ρ values and standard deviation of the different sample probes plotted versus the SL amplitude. The calculated T1ρ values of our new method are in good agreement with the reference TSE values (shaded area, mean ± std). c Heat maps of the dispersion offset \({\mathrm{T}}_{1\uprho }^{0}\) and the dispersion slope \({\mathrm{m}}_{1\uprho }\) calculated using a linear dispersion model (Eq. 3). The results of the radial Bloch sorting scheme and the TSE reference measurements are in good agreement and show a comparable image quality. Only at the edges of the sample tubes, some deviations are visible using the Bloch sorting scheme. This might be caused by data undersampling and the KWIC-filtered view-sharing technique

Full size image
Table 2 Dispersion analysis of the phantom experiments
Full size table

In vivo measurements

The calculated T1ρ maps of the N = 10 mice are shown in Fig. 8. A segmentation of the left ventricle according to the AHA model was carried out for each map. The results are listed in Table 3 for all animals (physiological parameters in the supplementary material, Online Table 2). Here, the breathing cycle was 1460 ± 154 ms and the cardiac cycle length was 137.3 ± 5.9 ms averaged over all animals. The global left ventricular T1ρ was 39.5 ± 1.2 ms for fSL = 1500 Hz. Slightly increased T1ρ values could be found in the segments 3 (41.7 ± 2.5 ms) and 4 (41.1 ± 2.1 ms). The comparison of the Cartesian reference measurement (37.8 ± 3.7 ms) with the accelerated radial method (38.3 ± 3.3 ms) revealed a small deviation of + 1.3%.

Fig. 8
figure 8

Results of in vivo T1ρ mapping in N = 10 mice. a Images of the acquired T1ρ maps (short-axis view, isotropic resolution 250 μm) with identical SL amplitudes of 1500 Hz. All maps were segmented in the left ventricle according to the AHA model. In animal 2, a map was also calculated from a fully sampled Cartesian data set. b The results of the segmentation are shown in a box plot. In the segments 3 and 4 the values are slightly increased. Segment 1 tends to have the lowest values

Full size image
Table 3 Results of in vivo T1ρ mapping in mice
Full size table

In Fig. 9 the results of the in vivo reproducibility study are depicted. Figure 9a shows the ten repetitive T1ρ maps, indicating a very good image quality. There are hardly any visual differences in the images. All structures are at the same position and the quantified T1ρ values appear to be almost identical. A ROI-based analysis of the left ventricular myocardium shows a mean T1ρ of 38.52 ms and only a small variation in the individual T1ρ values with a standard deviation of ± 0.54 ms (Fig. 9b). The maximum deviation observed in the separate experiments was 2.1%.

Fig. 9
figure 9

Results of the in vivo reproducibility measurements. a Images of the ten repeatedly acquired T1ρ maps (short-axis view, isotropic resolution 250 μm) with identical SL amplitudes of 1500 Hz. All T1ρ maps show the same image quality without any motion, flow or reconstruction artifacts. The dashed lines represent the chosen ROI for data analysis. b Mean (blue dots) and standard deviation (black lines) of the myocardial T1ρ values of the ten repetitive measurements. Mean T1ρ has been determined to 38.52 ± 0.54 ms (dotted red line and light red area)

Full size image

Figure 10 shows the results of the in vivo T1ρ dispersion quantification measurements. The T1ρ maps show a very good image quality with minimal motion, flow or reconstruction artifacts. The T1ρ dispersion map reveals a diagnostic image quality with little blurring at the myocardial borders. This might be caused by slightly different heart phases within the underlying T1ρ maps. The T1ρ dispersion slope has been determined to 4.76 ± 0.23 ms/kHz for left ventricular myocardium and 21.57 ± 2.56 ms/kHz for the left ventricular blood pool. The dispersion offset could be determined to be 32.73 ± 0.36 and 48.72 ± 3.41 ms, respectively.

Fig. 10
figure 10

In vivo T1ρ mapping including T1ρ dispersion quantification (short-axis view, isotropic resolution 250 μm). a Exemplary T1ρ maps acquired with different SL amplitudes (fSL = 750 Hz…2500 Hz). The increase of T1ρ with increasing fSL is clearly visible. All T1ρ maps show diagnostic image quality without any motion, flow or reconstruction artifacts. b Myocardial dispersion slope map. The image quality is also good with a slight blurring at the myocardial borders. The dashed lines represent the location of the chosen ROIs for the T1ρ dispersion quantification of myocardium and the blood pool. c Analysis of the dispersion behavior of myocardial tissue and blood. The T1ρ dispersion slopes are 4.76 ± 0.23 ms/kHz for left ventricular myocardium and 21.57 ± 2.56 ms/kHz for the left ventricular blood pool

Full size image

Discussion

In the current study, we introduced a fast cardiac T1ρ quantification technique that enables high-resolution in vivo T1ρ mapping and even T1ρ dispersion mapping in small animals. The method consisted of three concepts; Signal maximization using high flip angles, a specially adapted radial data sampling (which was optimized using Bloch simulations) and a KWIC-filtered view-sharing technique for efficient image reconstruction. With these concepts, an appropriate data undersampling (factor 32:13 ≈ 2.5, 13 preparations for radial sampling, 32 = 128/4 preparations for Cartesian sampling) is possible requiring only a fraction of the commonly needed data for the acquisition of T1ρ weighted images.

The optimized radial sequence ensures a high signal-to-noise-ratio and is combined with a very efficient sampling strategy that accelerates data acquisition while reducing the occurrence of artifacts and preventing incorrect T1ρ weighting. The main reason for this is that the pre-calculation of the expected signal intensity using Bloch simulations avoids strong changes in the signal levels of consecutive spokes. This results in an improved image quality with minor artifacts and an SNR increase of 120% in contrast to a conventional Serial sorting scheme of golden angles. Besides, the optimized Bloch sampling strategy reveals a higher accuracy and precision of the calculated T1ρ values compared to Serial sorting. The T1ρ quantification accuracy of Bloch sorting has proven in phantom measurements to be very high with a mean deviation of − 0.46 ± 1.84% compared to a TSE reference measurement. We observed the largest quantification errors in the phantom with 10% BSA concentration. Compared to the other phantoms, this has a much longer T1ρ relaxation time. Here the influence of the KWIC filter could possibly affect the contrasts. As a result, systematically shorter relaxation times were measured (≈ − 3%). However, this shift did not affect the quantification of the dispersion slope and offset and the T1ρ range of this phantom is hardly relevant for the characterization of myocardial tissue.

As described by Song et al. [31], the variation of signal intensities between the radial profiles possibly results in different relative proportions of high and low spatial frequencies for different weighted images. This effect depends on the design of the KWIC filter and was examined in the supplementary material. If the Nyquist criterion is not met in all regions of k-space (Nyquist factor fnyq < 1), undersampling artifacts and reduced SNR occur (Online Fig. 2). Opposed, excessive sharing of different contrasts (fnyq > 2) results in edge blurring and edge sharpening (Online Fig. 3). As the results in the supplementary material show, the T1ρ quantification accuracy is also affected by the choice of the KWIC filter design (Online Figs.4 and 5). However, in the range chosen in this work (fnyq = 1.1) the accuracy is stable and only shows systematic errors for the phantom with the highest relaxation time. Using moderate Nyquist factors fnyq = 1.0…1.5, the influence of the KWIC filter design on the quantification is small (± 0.3%) and can be neglected in the context of in vivo experiments.

In this work constant high flip angles (40°) were used for the gradient echo readout. For small flip angles, Bloch sorting can be avoided. However, in a previous work we could show that the resulting image quality is significantly reduced using an echo number based sorting and 10° flip angles, for example [36]. Yet, it is also possible to use ramped flip angles after the preparation. This could prevent the formation of signal plateaus and, if the RF pulse accuracy is well calibrated, the Bloch sorting could be omitted. However, different T1ρ weightings require different ramps, whereby the first readout must always be carried out using the same flip angle in order not to impair T1ρ quantification.

In the in vivo experiments, a good comparability of the T1ρ values in the left ventricle could be achieved for different animals. The T1ρ values determined at fSL = 1500 Hz varied between 37.9 and 41.8 ms. Musthafa et al. obtained values in the range 32…38 ms for fSL = 1250 Hz and B0 = 9.4 T [20]. The deviation could be explained by the T1ρ dispersion, since, according to basic relaxation theory, smaller values are expected for decreasing fSL. In addition, a high goodness of the mono-exponential T1ρ fit could be observed in most of the measurements. The values in the left ventricle were R2 = 0.993 on average. The slightly increased T1ρ values in the segments 3 (+ 5.6%) and 4 (+ 4.1%) may have been caused by partial volume effects, since the wall thickness was also reduced here. However, streaking and blurring artifacts could also have an impact. The fully sampled Cartesian measurement shows a slightly improved image quality compared to the accelerated radial technique. However, the quantification in the respective segments is in good agreement (+ 1.9 ± 3.3%) and the deviation in the global left ventricular ROI was only + 1.3%. It also needs to be considered that the cardiac phases were not exactly identical. The cycle length for the radial acquisition was 131.5 ± 1.2 ms and for the Cartesian measurement 153.2 ± 0.7 ms. This can also be seen in the contrast of the blood pool.

The in vivo reproducibility study showed, that our method enables a high degree of comparability for successive measurements, which is also crucial for dispersion imaging. This is a consequence of the radial readout that is more robust against moving tissues and blood flow. However, the main reason might be the very short measurement time of around 2.5 min. Alternative methods that require longer measurement times can lead to problems in the stability of the animal under investigation due to variations in the cardiac and respiratory cycle. This would require additional compensation and/or correction mechanisms, which is quite complex and hampers the practicality of cardiac T1ρ and/or T1ρ dispersion mapping.

Furthermore, a new approach for quantifying the dispersion behavior was introduced in this work. Due to the high acceleration of the radial data acquisition, the measurement could be accomplished in 20 min. In this total measurement time, only relatively small variations in the physiological parameters were observed. The mean cardiac cycle length was 141.5 ± 3.9 ms while the breathing cycle was 1430 ± 63 ms. In comparison, a fully sampled Cartesian method would have taken 50 min, potentially leading to higher variations. The determination of the dispersion slope m1ρ and the offset \({T}_{1\uprho }^{0}\) was based on the approach of Yin et al. [19], in which the difference \({T}_{1\uprho }^{0}-{T}_{1\uprho }({f}_{\mathrm{SL}})\) is presented as a potential myocardial fibrosis index. The calculation of a slope has an advantage over a differential measurement because the value is thus normalized. In addition, the determination of the dispersion slope based on eight independent T1ρ maps increases the accuracy compared to the previously described method of Yin et al. [19]. Nevertheless, it is necessary to ensure that the linear range of T1ρ dispersion is not exceeded, since it is known from relaxation theory that T1ρ dispersion is formally more complex than a linear model [10, 11]. In [19] the range fSL = 0…510 Hz was used for dispersion quantification, although it was not explained why this range is particularly suitable. Future studies must therefore primarily investigate which ranges of SL amplitudes are most suitable for calculating a fibrosis index.

For the translation of the method to human studies, the SAR (specific absorption rate) is a limiting factor. In this study, relatively high SL amplitudes were used, which show a low susceptibility to B0 inhomogeneities [27, 28]. These amplitudes are hardly feasible using radiofrequency amplifiers available on clinical scanners, which ultimately limits routine clinical applications [29]. The concept of Bloch sorting can, however, be transferred, whereby far more readouts would be possible after the preparation due to the lower heart rate in humans. Here the optimal flip angle arises with significantly lower values, which could possibly be utilized for SAR reduction.

A drawback of our study is that we are not yet able to present results in fibrotic myocardial tissue. The quantification of the T1ρ dispersion is of particular interest for future studies. Measurements in the infarct model, in which the dispersion is to be examined in detail, are planned. Besides, work is in progress to develop an extension of the sequence design presented in this study. Here, not only one T1ρ map should be acquired in a single scan, but the whole data set required for dispersion mapping, since separate maps can ultimately lead to quantification errors due to misregistration and cardiac motion. For this purpose, an optimization of the radial sampling pattern is also considered based on the Bloch sorting principle.

Conclusion

In conclusion, our new T1ρ quantification technique represents a reasonable tool for cardiac T1ρ mapping and T1ρ dispersion imaging. The method is feasible for improved cardiac tissue characterization and possibly enables enhanced diagnostics for several diseases as it might identify diffuse fibrosis. Due to its very short measurement time, stability and robustness, this method might be included in a common in vivo measurement protocol for small animals without major difficulty.