In the following, we give an overview of the implementation of the proposed open-source cMRF sequence and the reference spin-echo sequences. Image reconstruction and parameter estimation are explained and the experiments and analysis are described.

The design of the open-source cMRF sequence follows the design presented in Ref. [11]. An overview of the sequence is given in Fig. 1. 2D data acquisition is carried out over 15 heartbeats in a pre-defined mid-diastolic phase of the cardiac cycle. In three cardiac cycles, the data acquisition is preceded by an adiabatic inversion pulse. In nine cardiac cycles, a \(T_2\)-preparation pulse is applied before the data acquisition with three different echo times. Data acquisition is carried out with a variable density spiral.

Overview of the cMRF sequence. a In the first RR-interval an inversion pulse is applied, in the next RR-interval no preparation pulse is applied, and in the subsequent three RR-intervals \(T_2\)-preparation pulses with different echo times are applied (e.g., 30 ms, 50 ms and 100 ms). This five heartbeat pattern is repeated three times, leading to a total scan duration of 15 RR-intervals. b The flip angle (FA) is varied during the data acquisition between \(4^\circ\) and \(25^\circ\). c Data acquisition is carried out with a variable density spiral trajectory. Data is sampled during the spiral-out part of the trajectory (red) which is followed by a gradient rewinder (black). d The MLEV-4 type \(T_2\)-preparation blocks consist of a (90x, 180y, 180y, −180y, −180y, 270x, −360x) pulse pattern, where the ±180y pulses are realized using (±90x, ±180y, ±90x) composite pulses [20]. The \(T_2\)-preparation block shown here corresponds to an echo time of 30 ms

The code to generate the variable density spiral trajectories is based on the MATLAB code of Brian Hargreaves [21]. It was translated into Python and adapted to provide a suitable interface for application in PyPulseq [22]. The code allows k-space spacing (1/FOV) to be varied as a function of the k-space radius to increase the sampling efficiency by undersampling the outer portions of k-space. The maximum allowed gradient amplitude and slew rate values can directly be extracted from PyPulseq system limitation objects and are taken into account during the trajectory generation. To calculate the most accurate and time-efficient rewinder gradients while accounting for the gradient raster time, the PyPulseq function make_extended_trapezoid_area was used. This function allows the desired start and end amplitudes, as well as the desired gradient area (0th gradient moment), to be specified and returns the shortest possible extended trapezoid gradient that fulfills these parameters.

The spiral arms were rotated with respect to each other by an approximation of the golden angle [23] as suggested by Hamilton et al. [24]. In total, 48 golden angle spiral arms were pre-calculated with a constant angular increment of \(2\pi /48\), which served as a dictionary for the available trajectories. The starting angle for the first readout was \(0^\circ\). For the following readouts, the starting angle was calculated by adding a golden angle increment. A trajectory was then selected from the trajectory dictionary, which had the starting angle closest to the current starting angle calculated based on the Golden angle increment. As an example, the starting angle for the spiral arm with index 100 is \(100 \times 2\pi \times (1 - 2/(1 + \sqrt))\) which is \(71^\circ\). Therefore, the trajectory with the starting angle of \(67.5^\circ\) was selected for the readout with index 100 which is the closest entry to \(71^\circ\).

The MLEV-4 type \(T_2\)-preparation block shown in Fig. 1 features a (\(90x, 180y, 180y, -180y, -180y, 270x, -360x\)) flip angle pattern, with the \(180^\circ\) refocusing pulses (±180y) implemented via a (\(\pm 90x, \pm 180y, \pm 90x\)) composite pulse design [20]. This entire \(T_2\)-preparation block is modular, allowing for adjustments to RF pulse duration and echo time. Similarly, the \(T_1\)-preparation/inversion block is modular, enabling customization of inversion time, RF pulse duration, and spoiler gradient duration.

To support image reconstruction, relevant sequence details and the spiral trajectory for each acquisition are saved in an ISMRM raw data (*.MRD) file [25], created in parallel with the Pulseq sequence. After data acquisition on any MRI scanner, the raw data can be exported and converted to the MRD format. Using the Python code provided in our GitHub repository, both MRD files can then be merged into a single file containing all data and header information, as well as trajectory information required for image reconstruction.

For the reconstruction of the different qualitative images, we used a sliding window approach [26]. For each dynamic image, multiple spiral arms are combined and spiral arms are shared between neighboring images (sliding window). An extended phase graph (EPG) approach is used to calculate the signal evolution during the cMRF acquisition [27]. The cMRF sequence is triggered to the cardiac cycle and therefore the time between data acquisition and preparation blocks in the sequence depend on the heart rate. A separate dictionary was calculated for each subject. The signal for different \(T_1\) and \(T_2\) times is estimated for each RF pulse with \(T_1\) values selected as [50:10:2000 ms, 2020:20:3000 ms, 3050:50:5000 ms] and \(T_2\) values selected as [6:2:100 ms, 105:5:200 ms, 220:20:500 ms]. To take the sliding window reconstruction into account, the signal in each sliding window range is then averaged to obtain the final signal curves. This signal dictionary is then matched to the reconstructed images on a pixel-by-pixel basis. The \(T_1\) and \(T_2\) times corresponding to the highest dot product between image and dictionary signal curves are selected for the final quantitative maps.

The reference \(T_1\) and \(T_2\) values were estimated with spin-echo sequences using a Cartesian sampling pattern, which obtain a single readout line after the \(90^\circ\)–\(180^\circ\)-RF-pulses. Data are acquired at different time points after an adiabatic inversion pulse for the \(T_1\) estimation and with different echo-times for the \(T_2\) estimation, respectively. The quantitative maps are obtained using dictionary mapping with a dictionary containing \(T_1\) values in [50:10:2000 ms, 2020:20:3000 ms, 3050:50:5000 ms] for the \(T_1\) reference scan and \(T_2\) values in [6:2:100 ms, 105:5:200 ms, 220:20:500 ms] for the \(T_2\) reference scan.

In this study, we carried out three different sets of experiments to evaluate the reproducibility and accuracy of the open-source cMRF sequence across different scanners and to evaluate the sequence in in vivo experiments.

To demonstrate that the proposed open-source cMRF sequence provides accurate \(T_1\) and \(T_2\) quantification on different scanner platforms, we carried out phantom experiments using a T1MES phantom [19]. We used the same phantom on four 3 T Siemens scanners (Siemens Healthineers, Erlangen, Germany) with different software versions (Table 1).

We selected a lower maximum gradient strength and gradient slew rate (leading to longer gradient durations) to what would have been possible on most of the scanners to ensure that the sequence can be run with the exact same settings on each scanner. Values of 30 mT/m for the maximum gradient amplitude and 100 T/m/s for the maximum slew rate were used for all of these phantom experiments. Data acquisition was carried out with a field-of-view (FOV) of 128 \(\times\) 128 \(\hbox ^2\) with a pixel size of 1 \(\hbox ^2\) and a slice thickness of 8 mm. For the open-source cMRF sequence, the echo times of the \(T_2\)-preparation pulses were selected as 30 ms, 50 ms and 100 ms. An ECG-signal was simulated for a heart rate of 60 bpm. In each cardiac cycle, 47 spiral arms were acquired with 355 sampling points along each spiral arm. The TE and TR were set to 1.6 ms and 10 ms, respectively.

In addition to the open-source cMRF sequence, we also acquired data with the \(T_1\) and \(T_2\) spin-echo reference sequences. The FOV, voxel size and slice thickness were set to be the same as those used in the cMRF sequence. For the \(T_1\) spin-echo reference sequence, data at inversion times of 25 ms, 50 ms, 300 ms, 600 ms, 1200 ms, 2400 ms and 4800 ms were obtained with TE of 20 ms which took 1 h 59 min. For the \(T_2\) spin-echo reference sequence, echo times of 24 ms, 50 ms, 100 ms, 200 ms and 400 ms were obtained which took 1 h 25 min. The TR for both sequences was 8 s. Although the same phantom was used for all phantom scans, the \(T_1\) and \(T_2\) reference sequences were carried out at each scanner to compensate for potential relaxation time variations, such as those caused by temperature differences.

A direct reconstruction was used for image reconstruction, where each image was reconstructed by applying a density compensation function to the spiral k-space data followed by gridding on a Cartesian grid and Fast Fourier transformation. For the Cartesian spin-echo reference scans only the Fast Fourier transformation was needed. The data from the different receiver coils were combined by a weighted sum using coil sensitivity maps.

For all phantom experiments, parameter maps obtained with the open-source cMRF sequence were compared to the \(T_1\) and \(T_2\) maps obtained with the spin-echo reference sequences for each of the scanners separately. In each tube of the T1MES phantom, a region-of-interest (ROI) was manually delineated. The mean value and standard deviation within this ROI was compared. A linear fit was carried out and the relative root-mean-squared-error (RMSE) over all nine tubes was calculated between the cMRF and the reference sequences. Figures showing \(T_1\) and \(T_2\) maps follow the recommendation provided by Fuderer et al. [28].

In the second set of experiments, we obtained \(T_1\) and \(T_2\) maps in a healthy volunteer (male, 34 years). We also obtained \(T_1\) and \(T_2\) maps with a vendor-specific cMRF sequence. The sequence follows the method outlined in Ref. [11] and is available as a prototype sequence for Siemens scanners. These experiments were only carried out on scanner 1 as the vendor-specific cMRF sequence was only available there. The FOV was 300 \(\times\) 300 \(\hbox ^2\) with a resolution of 1.56 \(\times\) 1.56 \(\hbox ^2\) and a slice thickness of 8 mm. We maximized the peak gradient strength (180 mT/m) and gradient slew rate (150 T/m/s) and adapted the sequence parameters of the open-source cMRF sequence as much as possible to the vendor-specific cMRF sequence. \(T_2\)-preparation pulses with echo times of 30 ms, 50 ms and 80 ms, a TE of 1.6 ms and a TR of 6.6 ms was used. 47 spiral arms were acquired in each cardiac cycle with 355 readout points along a spiral arm. For the phantom scans, a heart rate of 60 bpm was simulated. For the in vivo scan, data acquisition was carried out during mid-diastole. A cine acquisition was carried out to determine the rest-period of the heart during diastole. The parameters of the vendor-specific cMRF sequence were the same, except for TE and TR which were 1 ms and 5.5 ms, respectively. In one cardiac cycle 45 spiral arms were obtained, each with 1388 sampling points. Fat-suppression was used for this sequence.

For the in vivo experiments, the spin-echo reference sequences could not be used due to their long acquisition times. Clinically recommended sequences were used instead [29]. For \(T_1\) mapping, a 5(3)3 Modified look-locker inversion recovery (MOLLI) sequence was used [30]. \(T_2\) maps were obtained using a \(T_2\)-prepared single-shot gradient echo sequence (T2prep) [31]. The echo times of the \(T_2\)-preparation blocks were 0 ms, 35 ms and 55 ms. Both sequences were ECG-triggered with the same trigger time as the cMRF sequences. MOLLI and T2prep are part of a software package for cardiac imaging provided by the vendor. \(T_1\) and \(T_2\) estimation was carried out on the scanner using the vendor software. The accuracy of sequences used for in vivo imaging was also evaluated in phantom experiments by comparing the results to the spin echo reference sequences. In the following we refer to the MOLLI and T2prep sequence as clinical sequences.

In the third set of experiments, we carried out a traveling volunteer study, where we obtained \(T_1\) and \(T_2\) maps with the open-source cMRF sequence in three volunteers on two different scanners. The scans were conducted on two successive days.

All in vivo scans of the second and third set of experiments were carried out on a short-axis orientation and each scan was obtained in a single breath-hold. These in vivo experiments were approved by the ethic committee of the responsible institution and carried out in accordance with relevant guidelines and regulations. The subjects gave written informed consent to participate in this study.

For the in vivo scans, each image was reconstructed using an iterative SENSE approach due to the larger FOV and hence higher level of undersampling artifacts [32]. The coil sensitivity maps needed for both reconstructions were calculated from the average image utilizing all acquired data with the method proposed by Inati et al. [33]. Image reconstructions and parameter estimations for the open-source cMRF sequence were performed offline using MRpro [34], a PyTorch-based MR image reconstruction and processing package. The processing time for a single cMRF in vivo acquisition took less than 2 min (50 s for iterative image reconstruction, 50 s for dictionary calculation, 1 s for dictionary matching) on a PC with 12 CPUs and 192 GB RAM without GPU acceleration.