Electrocardiogram-less, free-breathing myocardial extracellular volume fraction mapping in small animals at high heart rates using motion-resolved cardiovascular magnetic reesonance multitasking: a feasibility study in a heart failure with preserved ejection fraction rat model

Background Extracellular volume fraction (ECV) quantification with cardiovascular magnetic resonance (CMR) T1 mapping is a powerful tool for the characterization of focal or diffuse myocardial fibrosis. However, it is technically challenging to acquire high-quality T1 and ECV maps in small animals for preclinical research because of high heart rates and high respiration rates. In this work, we developed an electrocardiogram (ECG)-less, free-breathing ECV mapping method using motion-resolved CMR Multitasking on a 9.4 T small animal CMR system. The feasibility of characterizing diffuse myocardial fibrosis was tested in a rat heart failure model with preserved ejection fraction (HFpEF). Methods High-salt fed rats diagnosed with HFpEF (n = 9) and control rats (n = 9) were imaged with the proposed ECV Multitasking technique. A 25-min exam, including two 4-min T1 Multitasking scans before and after gadolinium injection, were performed on each rat. It allows a cardiac temporal resolution of 20 ms for a heart rate of ~ 300 bpm. Myocardial ECV was calculated from the hematocrit (HCT) and fitted T1 values of the myocardium and the blood pool. Masson's trichrome stain was used to measure the extent of fibrosis. Welch’s t-test was performed between control and HFpEF groups. Results ECV was significantly higher in the HFpEF group (22.4% ± 2.5% vs. 18.0% ± 2.1%, P = 0.0010). A moderate correlation between the ECV and the extent of fibrosis was found (R = 0.59, P = 0.0098). Conclusions Motion-resolved ECV Multitasking CMR can quantify ECV in the rat myocardium at high heart rates without ECG triggering or respiratory gating. Elevated ECV found in the HFpEF group is consistent with previous human studies and well correlated with histological data. This technique has the potential to be a viable imaging tool for myocardial tissue characterization in small animal models.

Rodent models are widely used in preclinical studies of myocardial diseases because of the short development period, availability of genetically modified disease models, and low cost. However, it is technically challenging to acquire high-quality T 1 and ECV maps in small animals because of their high heart rates (often > 300 bpm) and high respiration rates (around 60 cpm). Therefore, CMR quantification of ECV is still not well established despite the important unmet needs in research. Several studies have been done to improve CMR T 1 mapping and/or ECV measurement in small animals. Coolen et al. proposed a 3D T 1 mapping method of the mouse heart using variable flip angle (VFA) analysis [7]. The method has the ability to detect regional differences in myocardium with excellent repeatability, but VFA-based methods have inherent problems with B 1 inhomogeneity, and the long scan time (more than 20 min) makes it impractical to be used in pre-and post-Gd studies. Messroghli et al. acquired myocardial T 1 and ECV maps in rats from a single (unsegmented) dataset using a small animal Look-Locker inversion recovery (SALLI) method [8,9]. It was able to reconstruct both cine CMR images and T 1 maps, and showed the feasibility to detect diffuse myocardial fibrosis, while the spatial resolution was limited due to signal-to-noise ratio (SNR) consideration. Segmented multi-shot FLASH methods were then proposed to acquired images with higher resolution [10,11].
In all previous studies with a Look-Locker scheme [12,13], electrocardiogram (ECG) triggering was used to monitor cardiac motion [8][9][10][11]. However, ECG triggering at high field strengths is unreliable due to elevated magnetohydrodynamic effects, which can introduce trigger-related motion artifacts and blurring effects [14,15]. High heart rates may make the situation even worse. Respiratory navigation was used in some studies [10], however, most studies simply used signal averaging, resulting in image blurring [8,11]. Furthermore, ECG and respiratory gating setup leads to complicated workflow.
In this study, we developed a ECG-less, free-breathing ECV mapping method using CMR Multitasking [16,17] (hereinafter abbreviated as ECV Multitasking) on a 9.4 T small animal CMR system. It allows continuous acquisition without ECG triggering or respiratory gating. The pre-and post-Gd data were acquired separately but reconstructed jointly, to allow image co-registration and direct ECV mapping. The feasibility of characterizing diffuse myocardial fibrosis was tested in a rat hypertensive heart failure model with preserved ejection fraction (HFpEF), which has been shown to have increased left ventricular (LV) interstitial fibrosis [18][19][20], and has been recapitulated in recent studies [21][22][23].

Animal model
All animal experiments were approved by the Cedars-Sinai Institutional Animal Care and Use Committee. Dahl salt-sensitive (DSS) rats can develop hypertension followed by HFpEF on a high-salt diet [18,19]. In this model, male DSS rats (Charles River Laboratories, Wilmington, Massachusetts, USA) were normally fed (0.3% NaCl) until the age of 7 weeks. Rats were then randomly assigned to a high-salt (HS) diet group (8% NaCl) to induce HFpEF or a normal-salt (NS) diet group (0.3% NaCl) to serve as controls, until the age of 14 weeks [24]. HS rats with heart failure symptoms (including decreased activity, cachexia, labored breathing, and body edema) and echocardiographic evidence for diastolic dysfunctionwere diagnosed as HFpEF.
Control rats (n = 9; weight, 335 ± 44 g) and HS fed rats diagnosed with HFpEF (n = 9; weight, 286 ± 47 g) were imaged. Imaging experiments and all measurements were done between the age of 14 weeks and 15 weeks. After the imaging study, animals were euthanized, and the hearts were excised. Mid-ventricular heart tissues of the control rats (n = 9) and HFpEF rats (n = 9) were sectioned and stained with Masson's trichrome staining.

CMR protocol
All CMR data were acquired on a 9.4 T preclinical system (BioSpec 94/20 USR; Bruker Biospin, Billerica, Massachusetts, USA) using a single-channel volume coil. The T 1 Multitasking sequence was implemented in Paravision 5.1 by modifying the built-in FLASH sequence.
The rat was anesthetized before the scan, and the tail vein was cannulated for later Gd injection within the scan. During the scan, anesthesia was maintained by ventilation human studies and well correlated with histological data. This technique has the potential to be a viable imaging tool for myocardial tissue characterization in small animal models.
Keywords: Cardiovascular MR, T 1 mapping, Extracellular volume fraction, HFpEF with 1.5% isoflurane-oxygen. After the scan, hematocrit (HCT) level was measured for ECV calculation. Figure 1a shows the imaging workflow for one rat study. A self-gated (IntraGate) localizer was used to select a slice with a mid-cavity short-axis LV view. The T 1 Multitasking sequence [17] (Fig. 1b) was then performed. Gd contrast agent (Gadavist, 0.2 mmol/kg; Bayer Schering Pharma, Berlin-Wedding, Germany) was manually injected immediately afterwards. Fifteen minutes later, the T 1 Multitasking sequence was repeated on the same slice using identical imaging parameters. The total exam time was ~ 25 min.
Data acquired from pre-and post-Gd T 1 Multitasking sequence were reconstructed and analyzed jointly, as described in the following sections. The joint pre-and post-Gd reconstruction can improve respiratory and cardiac binning and promote image co-registration in native and post-Gd T 1 fitting.

Image reconstruction
Images from the proposed ECV Multitasking protocol are represented as a high-dimensional image A(x, t c , t r , t T 1 ) with 2 spatial dimensions and 3 temporal dimensions (cardiac phase t c , respiratory phase t r , and T 1 recovery time t T 1 ). The high-dimensional image can be discretized and viewed as a low-rank tensor A , and thus partially separable [25], i.e.
where columns of U x represents spatial basis functions; columns of U c , U r , and U T 1 contain cardiac, respiratory, and T 1 recovery temporal basis functions, respectively; ⊗ denotes the Kronecker product; and A (1) and C (1) are mode-1 matricization of the image tensor A and the core tensor C respectively [26].
T , in which and U x contain separate temporal and spatial bases respectively. Therefore, can be first recovered using only the training data d tr , which is frequently sampled in time, via Blochconstrained low-rank tensor completion followed by high-order singular value decomposition (HOSVD) [16,27]. With determined, the spatial basis U x can then be reconstructed from the imaging data d im by solving the following problem: where E is the signal encoding operator, including Fourier transformation and coil sensitivities weighting (optional), Imaging workflow and T 1 Multitasking sequence diagram. a Imaging workflow: A self-gated (IntraGate) localizer was used to select a mid-cavity short-axis slice of the left ventricle (LV). The CMR Multitasking T 1 mapping sequence (T 1 Multitasking) was then performed. Gd contrast agent (Gadavist, 0.2 mmol/kg) was manually administered. 15 min after injection, T 1 Multitasking was repeated. The total exam time was around 25 min. b Sequence design: T 1 Multitasking was performed using a continuous FLASH acquisition with repeated non-selective inversion recovery (IR) magnetization preparation pulses. Odd-numbered readouts followed randomized Gaussian-density sampling in the phase encoding ( k y ) direction (used as the imaging data), and even-numbered readouts collected the k-space center line ( k y = 0 , used as the subspace training data) is the undersampling operator, and R is a regularization functional which was chosen here as a wavelet sparsity penalty in order to additionally exploit compressed sensing [28].
A detailed description of the basic image model and reconstruction scheme used in CMR Multitasking can be found in previous work [16,17,29]. In this work, it is divided into the following steps:

Determination of T 1 recovery basis functions in U T 1 with Bloch simulation
The inversion recovery (IR)-prepared T 1 recovery process is modeled by Bloch simulations to determine U T 1 . In this work, pairs of pre-and post-Gd T 1 recovery curves were jointly simulated as shown in Fig. 2, for the purposes of joint physical modeling during respiratory and cardiac binning (Step 3) and multidimensional image reconstruction (Step 4). The detailed information of the T 1 recovery modeling is available in Appendix.
A dictionary was generated with this model using 101 T 1 values ranging from 100 to 3000 ms ( T 1,pre > T 1,post combinations only), 15 FLASH flip angles from 0.5 to 7.5 • , and 21 inversion efficiency values from -1 (complete inversion) to 0 (no inversion). The singular value decomposition (SVD) of this dictionary yielded the T 1 recovery basis functions U T 1 .

Real-time image reconstruction
A "real-time" (i.e., ungated before sorting into multiple time dimensions) image was first generated for motion identification: where rows of rt correspond to the real-time temporal basis functions, and columns of U x,rt correspond to the spatial basis functions. (3) The real-time temporal basis rt was estimated from the SVD of the training data d tr . Then the spatial coefficients U x,rt were recovered by solving the least-squares optimization problem:

Respiratory and cardiac binning
Different time points of X rt were assigned to multiple respiratory motion states ("respiratory bins") and cardiac motion states ("cardiac bins") based on the real-time images generated in Step 2. A modified k-means clustering method was used to automatically group the data into different bins, incorporating the predetermined low-rank T 1 recovery model in U T 1 to address the contrast change of X rt from T 1 recovery [16].

Tensor formation and multidimensional reconstruction
After cardiac and respiratory binning, each readout time point was assigned three temporal indices: cardiac phase, respiratory phase, and T 1 recovery index. A 4-way training data tensor D tr with one k-space readout dimension and three temporal dimensions was recovered from the subspace training data d tr by solving a Bloch-constrained low-rank tensor completion problem.
After D tr is completed, the cardiac basis functions U c , the respiratory basis functions U r , and the core tensor C can be extracted from the HOSVD of D tr , fully determin- T . The spatial coefficients U x can then be reconstructed by solving the problem in Eq. (2).

Parameter fitting and image analysis
One specific respiratory phase and cardiac phase were selected for T 1 and ECV map generation, corresponding Fig. 2 Illustration of joint pre-and post-Gd T 1 recovery modeling. Representative T 1 recovery curves from pre-and post-Gd acquisition are shown separately. Black line readouts represent image data, while green line readouts represent training data. A pair of single blocks from pre-and post-Gd acquisition are Bloch simulated together, with the same FLASH flip angle α and inversion efficiency B to end-expiration and end-diastole respectively. A joint pre-and post-Gd T 1 recovery model was also used in pixel-wise T 1 fitting. By fitting this joint model (see Appendix), we can get the pre-and post-Gd T 1 maps simultaneously.

Myocardial T 1 maps
Parameter fitting was done using the joint pre-and post-Gd T 1 recovery model in Eq. (11). T 1,pre , T 1,post , α, B, and M 0 z were fitted using the 832 inversion time images at the selected reparatory and cardiac phase. The pre-and post-Gd myocardial T 1 maps were directly generated from the T 1,pre and T 1,post fitting results.

Blood T 1 fitting
First, the pixels to be used for blood T 1 fitting were automatically selected based on the thresholding of the fitted M 0 z map. Then, the average T 1 recovery curve from these pixels was calculated, weighted by the M 0 z map. Finally, parameters were fitted to the average blood T 1 recovery curves using the model in Eq. (11) without Look-Locker correction (i.e., with M n M 0 z , T 1 , α = 0, B ) to better model the inflow of unexcited blood spins into the imaging slice [17]. The fitted T 1,pre and T 1,post values were then used as the pre-and post-Gd blood T 1 values.

ECV calculation
The ECV map was generated according to the following equations: in which R 1,myo is the pixel-wise R 1 changes, and R 1,blo is the change of blood R 1 , and HCT is the hematocrit value. All statistical analyses were performed for septal ECV, i.e., where R 1,myo was calculated as the mean value within the septal myocardium. All image reconstruction and curve fitting was done in MATLAB 2018a (MathWorks, Natick, Massachusetts, USA).

Histological analysis
Masson's trichrome staining was used to measure the extent of fibrosis. Mid-ventricular heart tissues were sectioned and stained following the manufacturer's protocol (Sigma-Aldrich, St. Louis, Missouri, USA). Quantitative histological analysis was done with ImageJ 1.52a (National Institutes of Health; http://image j.nih.gov/ ij). In each section, the extent of myocardial fibrosis was quantified by the percentage of total fibrosis area, which was calculated as the number of blue-stained pixels divided by the total number of pixels in the ventricular area (excluded: intramural vascular structures, perivascular collagen, endocardium, and LV trabeculae) [20,30]. For each rat, the average percentage from five different sections was reported.

Statistical analysis
Welch's t-test was performed to compare ECV values between the control group and HFpEF group. The Pearson coefficient was measured to evaluate the correlation between ECV values and quantitative fibrosis percentages. A two-tailed value of P < 0.05 was considered to be statistically significant. Statistical graphs were generated using GraphPad Prism 8 (GraphPad Software, San Diego, California, USA) and Excel (Microsoft Corporation, Redmond, Washington, USA).

Results
The reference T 1 of the phantoms measured by RARE-VTR ranged from 333 to 2433 ms, with a corresponding R 1 range from 0.41 to 3.00 s −1 . Figure 3 shows the linear regression result of R 1 values measured by Multitasking and RARE-VTR, in which a strong linear relationship can be found (R = 0.9983, P < 0.0001). Five respiratory bins and ten cardiac bins were used for respiratory and cardiac binning, corresponding to a cardiac temporal resolution of 20 ms for a heart rate of ~ 300 bpm. Figure 4 shows the images of the heart from a representative healthy control rat right before the IR preparation pulse, generated from different cardiac and respiratory phases, i.e. A x, n c , n r , n T 1 n c =1:10,n r =1,n T 1 =N =416 (Fig. 4a) and A x, n c , n r , n T 1 n c =6,n r =1:5,n T 1 =N =416 (Fig. 4b). The diastole (Fig. 4a6) and the systole (Fig. 4a1), or the end-expiration phase (Fig. 4b1) and the end-inspiration phase (Fig. 4b5) can be clearly differentiated in the figure. Movies showing the dynamic images from (1) different respiratory bins, (2) different cardiac bins, (3) native T 1 recovery process, and (4) post-Gd T 1 recovery process are also available in Additional file 1. Figure 5a, b shows representative native and post-Gd T 1 maps from the control group and heart failure with preserved ejection fractino (HFpEF) group. T 1 maps were smooth and homogeneous within the left ventricular (LV) myocardium, except for areas in the inferior and lateral walls with imperfect fitting results (as indicated by white arrows in Fig. 5a1, a2). The native myocardial T 1 values (mean ± SD, in ms) were 1662 ± 152 in the HFpEF group (vs 1534 ± 151 in the control group; P = 0.09). Figure 5c shows the corresponding ECV maps. In this study, signal from the aorta rather than the LV was usually selected for blood T 1 fitting, because of the signal loss resulting from LV inflow in this single slice setup; and the septal myocardium areas selected for R 1,myo calculation were indicated by red arrows. Welch's t-test showed that ECV was significantly higher in the HFpEF group (22.4% ± 2.5% vs. 18.0% ± 2.1%; P = 0.0010, Fig. 5d).
In eight rats (4 control + 4 HFpEF), an additional post-Gd T 1 Multitasking sequence was run 10 min after the Gd injection, with identical imaging parameters. ECVs were also calculated using this separate post-Gd T 1 measurement. Figure 6 shows the Bland-Altman plot comparing the ECVs measured from the 15-min post-Gd T 1 (ECV 15 ) and the 10-min post-Gd T 1 (ECV 10 ). The interclass correlation coefficient (ICC) was 0.817. The root-mean-square within-subject standard deviation was 1.4, yielding a coefficient of variation of 6.7%. There was no significant difference between ECV 15 and ECV 10 (P = 0.66). Figure 7a, b shows representative Masson trichromestained sections of a HFpEF and control rat. The myocardial fibrosis can be clearly seen as the diffused dark blue areas in Fig. 7b. The extent of fibrosis significantly increased in HFpEF hearts (11.7% ± 1.9% vs 3.4% ± 0.8%, P < 0.001) Fig. 7c shows the relationship between the ECV value and the extent of fibrosis, in which a moderate correlation can be found (R = 0.59, P = 0.0098).

Discussion
We developed a novel ECV Multitasking protocol which can map ECV at high heart rates without ECG triggering    [18,19], this newly developed technique can easily be extended to other small animal applications. More importantly, the non-invasive quantitative imaging protocol not only provides a diagnostic tool, but also a method for longitudinal therapy monitoring of the same subject [31].
In the control group, native T 1 values were 1534 ± 151 ms and ECV values were 18.0% ± 2.1% (both calculated from the septal area). Previous rodent studies have reported variable myocardial T 1 values, depending on field strength and T 1 encoding schemes, both of which are known to affect T 1 estimates. At 7.0 T, T 1 values were reported as 1620 ms (Vandsburger et al.) and 1638 ms (Zhang et al.) [10,32]. At 9.4 T, reported T 1 values ranged from ~ 1200 ms to 1764 ms (Kim et al. Li et al. and Coolen et al.) [7, 33,34]. ECV, a physiological quantity in principle unaffected by field strength or T 1 encoding scheme, allows more direct comparison. In healthy rats, myocardial ECV was previously reported as 16% [35], 18% [36], 17.2% [9], and 15.5% [34], comparable with the 18.0% value measured in our study. A direct comparison to these previous methods on the same scanner and in the same subjects was not performed in the present study, but would be a useful future validation step.
One major technical novelty of this work compared with previous CMR Multitasking work is the joint reconstruction of pre-and post-Gd data and the paired modeling of the pre-and post-Gd T 1 recovery. This has two primary benefits. First, it allows the respiratory and cardiac binning to be done in the concatenated pre-and post-Gd data. In this way, the data with a similar motion state-no matter whether from pre-Gd or post-Gd-will be clustered into exactly the same respiratory or cardiac bin, removing the need for respiratory co-registration between pre-and post-Gd T 1 maps. Second, it can improve the robustness of T 1 fitting by constraining pre-and post-Gd images to share the same thermal equilibrium magnetization M 0 z , FLASH flip angle α and inversion efficiency B , exploiting internal a priori knowledge and reducing the total number of fitting parameters from 8 to 5.
The phantom study found that the Multitasking based T 1 fitting results showed a strong linear correlation with the RARE-VTR reference. If calibrated using the linear correlation ( ∼ R1 = LC(R 1 ) = 1.166R 1 − 0.331 ), the relative R 1 difference between Multitasking method and RARE-VTR reference would be 2.7% ± 3.6%. Note that ECV measurements are insensitive to linear transformation: the ECV value is calculated by the ratio of R 1 changes of myocardial and blood pool, so both the intercept (− 0.331) and slope (1.166) of the linear transformation LC will be canceled, leaving the final ECV quantification result unchanged.
Ten bins were used to separate different cardiac phases in the current protocol for CMR ECV Multitasking reconstruction, corresponding to a cardiac temporal resolution of 20 ms for a heart rate of 300 bpm. As shown in Fig. 2, the diastolic and systolic phases can be differentiated at this temporal resolution. However, for higher heart rates, such as those in mice (around 450 bpm), more cardiac bins may be required, which may have an impact on SNR and ECV homogeneity. A nearly-significant negative correlation (R = − 0.44, P = 0.06) was found between the septal SNR of the ECV map and the heart rate (the scatterplot can be found in Additional file 1). This may be an effect of increased intra-bin motion at higher heart rates. More cardiac bins can possibly be achieved by shortening the echo spacing, or by acquiring k-space training line and imaging lines in separate echoes after each FLASH flip angle [37].
In this work, blood T 1 was extracted by fitting the weighted average T 1 recovery curve of the selected blood pool. The FLASH flip angle α in Eq. (8) was fixed as 0 to remove Look-Locker correction and better model the inflow of unexcited blood spins. This approach is particular to 2D slice-selective excitation, for which blood spins flow through the slice quickly enough that Look-Locker correction is not required. For a volumetric 3D variation of our method, the blood spins would be driven towards steady-state, and the Look-Locker correction could be retained in blood.
T 1 inhomogeneity was present in the inferior and lateral walls due to potential B 1 issues at high B 0 field, as visible in Fig. 4. In diffuse myocardial disease, it does not affect the ECV calculation from septal area. However, it will affect the accuracy of the ECV map in the inferior and lateral walls. This may reduce the reliability of the method in the diagnosis of focal myocardial disease, such as focal myocardial infarction. Further sequence improvements, such as a dual-flip-angle acquisition [38], should be made to address this problem.
The method described here did not include bulk motion compensation, as anesthesia and additional immobilization setup did not result in any apparent bulk motion. However, bulk motion compensation is compatible with the CMR Multitasking framework, as described in [39], where a translational inter-bin registration to k-space data was applied prior to joint reconstruction. Similar motion compensation could be incorporated into the proposed method if used with different experimental setups more susceptible to bulk motion.
Currently, a total of 85 IR preparation pulses were applied in each T 1 Multitasking module, resulting in a single-slice scan of just over 4 min. However, the minimum required scan time was not systematically explored, so further shortening of the scan time may yet be possible.
Future work will include optimizing the number of cardiac bins, expanding spatial coverage with 3D volumetric imaging, as well as addressing B 1 and B 0 issues affecting T 1 homogeneity in the inferior and lateral walls. Further imaging tests will also include using ECV Multitasking for longitudinal therapy monitoring, such as the cardiosphere-derived cell (CDC) treatment of HFpEF [31].

Conclusions
We developed an ECG-less, free-breathing CMR Multitasking ECV mapping method at high heart rates. The pre-and post-Gd data were concatenated and reconstructed using a joint T 1 recovery model to generate ECV maps. Elevated ECV found in the HFpEF group agrees well with previous human studies and shows a moderate correlation with the histological data. This technique can serve as a viable CMR imaging tool for myocardial tissue characterization in small animal models.
Additional file 1: Page 1. The four videos show reconstructed results of (1) steady-state images from five respiratory bins, (2) steady-state images from ten cardiac bins, (3) native inversion recovery process, and (4) post-Gd inversion recovery, respectively. Page 2: The scatterplot between the SNR of the ECV map and the heart rate. The SNR was measured as mean/ std of the ECV map within the septal area. A nearly-significant negative correlation (R = -0.44, P = 0.06) can be found between the SNR and the heart rate. z-magnetization at true steady-state ( t → ∞ ); M N +1 is the actual final z-magnetization in each recovery period incorporating incomplete relaxation; and B ∈ [−1,0] represents the inversion efficiency of the IR pulse. Our pulse sequence used N = 416 , TR = 7 ms; α was prescribed as 5° but was not assumed to be perfectly homogenous during Bloch simulations and during fitting.
Combining Eqs. 7, 8, 9 yields the combined z-magnetization equation Assuming that the B 1 field does not change between preand post-contrast acquisitions, the unknown actual FLASH flip angle α and inversion efficiency B can be assumed to be constant for each pixel. Therefore the paired signal model concatenating the pre-and post-Gd signal evolution is: