Quantitative Magnetic Resonance Imaging by Nonlinear Inversion of the Bloch Equations

Purpose: Development of a generic model-based reconstruction framework for multi-parametric quantitative MRI that can be used with data from different pulse sequences. Methods: Generic nonlinear model-based reconstruction for quantitative MRI estimates parametric maps directly from the acquired k-space by numerical optimization. This requires numerically accurate and efficient methods to solve the Bloch equations and their partial derivatives. In this work, we combine direct sensitivity analysis and pre-computed state-transition matrices into a generic framework for calibrationless model-based reconstruction that can be applied to different pulse sequences. As a proof-of-concept, the method is implemented and validated for quantitative $T_1$ and $T_2$ mapping with single-shot inversion-recovery (IR) FLASH and IR bSSFP sequences in simulations, phantoms, and the human brain. Results: The direct sensitivity analysis enables a highly accurate and numerically stable calculation of the derivatives. The state-transition matrices efficiently exploit repeating patterns in pulse sequences, speeding up the calculation by a factor of 10 for the examples considered in this work, while preserving the accuracy of native ODE solvers. The generic model-based method reproduces quantitative results of previous model-based reconstructions based on the known analytical solutions for radial IR FLASH. For IR bSFFP it produces accurate $T_1$ and $T_2$ maps for the NIST phantom in numerical simulations and experiments. Feasibility is also shown for human brain, although results are affected by magnetization transfer effects. Conclusion: By developing efficient tools for numerical optimizations using the Bloch equations as forward model, this work enables generic model-based reconstruction for quantitative MRI.


Introduction
Conventional quantitative magnetic resonance imaging is based on a two-step process, where first intermediate images are reconstructed and then physical models are fitted pixel-wisely to obtain parameter maps.Acquiring a sufficient amount of high quality images with carefully designed contrasts is required for achieving a good fit.For this reasons, these methods are too slow for many clinical applications.In contrast, nonlinear model-based reconstruction methods formulate image reconstruction as a single inverse problem.They exploit a physical model of the measurement process and directly estimate quantitative parameter maps from k-space.Thus, they make optimal use of the available data and enable highly efficient parameter mapping from signals acquired with sequences that make use of transient magnetization dynamics [1,2,3,4,5].These techniques have two problems: They are computationally demanding and they need to be specially designed for each application.
Alternatively, fingerprinting [6] uses a lookup dictionary obtained by Bloch simulations to map the pixels of intermediate images computed directly from undersampled data to quantitative parameter maps.This enables multi-parametric mapping with high acceleration in a flexible and computationally efficient framework, but is not optimal due to its lack of a least-squares data consistency term.Subspace models can be exploited for a more efficient mapping by approximating the physical signal with a larger linear subspace.They reduce the computational demand of the reconstruction very efficiently [7,8,9,10,11], but are still not optimal because a linear subspace is used to approximate the manifold of possible signals.For complicated spin dynamics a larger number of subspace coefficients may be needed to accurately represent the signal, rendering subspace methods less efficient [5].
The aim of this work is to develop a generic framework for nonlinear modelbased reconstruction with accurate signal models for different MRI sequences even with complicated spin dynamics.The generalization of the forward model then allows the use of optimized sequences for which no analytical expression for their signal can be derived.Fundamentally, such a generic method requires efficient techniques to compute the partial derivatives of the Bloch equations.So far, two different methods were used in MRI.First, symbolic derivatives can be calculated for analytical solutions of the Bloch equations for special sequences [12,13,14] and this can be generalized to chains of small blocks using automatic differentiation [15,16].These methods require idealized assumptions such as hard pulse approximation and perfect inversion.In more complicated scenarios with long pulses or imperfect inversion they require high discretization rates or suffer from errors.Second, derivatives can be estimates using difference quotients [17,18].This method fully exploits the generality of full Bloch simulations, but is computational expensive and requires careful tuning to balance accuracy and noise amplification.
To overcome these limitations, this work uses a direct sensitivity analysis [19] to compute the derivatives of the Bloch equations for arbitrary sequences with high accuracy by solving an extended system of ordinary differential equations (ODE).The technique is validated using an analytical model of an IR bSSFP sequence and is compared to results obtained from difference quotients methods.To further improve computation speed, pre-computed state-transition matrices are applied to arbitrary initial conditions solving the required ODEs for all repeating parts of an MRI sequence efficiently.They are validated by comparing them to the direct application of a Runge-Kutta ODE solver.We further integrate both techniques in a nonlinear model-based reconstruction with integrated calibration-less parallel imaging.For IR FLASH, we show that the methods reproduce the results of an analytical model.In a numerical and measured phantom study with an IR bSSFP sequence we refine the flexible forward model of the generic model-based reconstruction to include realistic simulations with slice-selective excitations and hyperbolic secant inversion pulses.Thus, we show that the reconstruction quality benefits much from the more physically accurate modelling leading to accurate T 1 and T 2 parameter maps.Finally, we test the developed technique on in-vivo brain data from a healthy volunteer.

Theory
In the following, we briefly explain the concepts of a direct sensitivity analysis and its application to the Bloch equations (SAB).We then describe how state-transition matrices (STMs) can be used to accelerate the solution of the ODEs.Afterwards, both methods are integrated into a nonlinear model-based reconstruction method.

Sensitivity Analysis of the Bloch Equations
We consider the temporal evolution of a magnetization vector M (c, t) depending on a vector of parameters c and time t.The temporal behaviour of its components M i (c, t) is described by the Bloch equations as a system of ODEs where f defines the dynamics.The partial derivative of the component M i with respect to the parameter c j defines the (i, j)-th entry of the sensitivity matrix Z(t).
Using direct sensitivity analysis [19] one obtains Z(t) by solving an additional set of ODEs.Assuming that the partial and ordinary derivatives interchange, the time derivative of the (i, j)-th entry of Z is Substituting dMi(c,t) dt then yields With the chain rule, the resulting ODE becomes where ∂fi(M (c,t),c,t) ∂Mj describes the (i, j)-th element of the Jacobian J i,j .This can be written compactly for the sensitivity matrix Z(t) as If a direct sensitivity analysis is applied to the Bloch equations for the parameters R 1 , R 2 and B 1 , the ODE in Eqs.6 describing the temporal evolution of the sensitivities becomes depending on the x, y and z components of the magnetization M , the magnetic fields B z and B 1 as well as the RF pulse phase φ.Eq. 7 is solved jointly with the Bloch Eqs. 1 which provide the time-dependent solutions for M x , M y , and M z .

State-Transition Matrices
By embedding the magnetization vector into a four-dimensional space we obtain a formulation of the Bloch Eqs. 1 as a system of homogeneous ODEs with the system matrix at location r depending on time t, the z-gradient G z and the magnetic fields B x,y .The Bloch equations can be solved directly for time-dependent coefficients A(t) using standard ODE solvers.Here we describe the pre-computation of STMs as more efficient way to solve the equations for MRI sequences with repeating patterns.A STM S t1→t2 describes the evolution of an arbitrary starting magnetization M (c, t 1 ) for the time span from t 1 to t 2 including all effects from relaxation and time-dependent external RF fields.This compresses the temporal evolution to a single matrix multiplication The computation of S t1→t2 is based on derivation of Eq. 11 Using the Bloch Eqs. 9 to replace the time derivative on the left side and using dM (t1) dt2 = 0 for the right, we obtain By using Eqs.11 and switching both sides we obtain As this holds for arbitrary M (c, t 1 ) and by renaming t 2 as t a system of ODEs for the entries of the STM is derived.This ODEs 15 can be solved for estimating S t1→t2 column-wisely with an ODE solver [24,25] with initial conditions The solution of the state-transition ODE in Eqs. 15 can be formally defined as a time ordered exponential ≡ lim This links the proposed technique to approximation methods based on matrixexponentials computed using discretized sampling [26,27].This technique is not limited to the Bloch Equations, but can be extended to also include the sensitivity analysis for the three partial derivatives R 1 , R 2 and B 1 .This is further described in Appendix A.

Bloch Model-Based Reconstruction
In the following, we integrate the generic Bloch operator B into a nonlinear model-based reconstruction framework with non-Cartesian, calibrationsless, parallel imaging and compressed sensing as illustrated in Figure 1.The reconstruction method solves the nonlinear inverse problem for the maps x = (x p x c ) T with the physical parameters T and coil sensitivities x c = (c 1 . . .c N ) T by optimizing x = argmin x y − A(x) 2  2 + αQ(x c ) + βR(x p ) .
Eq. 19 includes the forward-operator A, the measured data y, the Sobolev norm Q with its regularization parameter α to enforce the smoothness of coil profiles [28] and B 1 maps.A joint sparsity constraint R is applied to the other parameter maps [29,30].The full forward operator is A = PF CB.It is solved by the Iteratively Regularized Gauss-Newton Method (IRGNM) with the Jacobian DA(x n ) and the regularization parameters α n = α 0 • q n and β n = β 0 • q n at the nth iteration step.Here, C is the nonlinear parallel imaging operator combining the signal with the coil profiles, F represents the Fourier operator, P the sampling pattern.The generic operator B takes information about the applied sequence and outputs the simulated signal based on the STM technique.The partial derivatives of B are calculated using the direct sensitivity analysis.The derivatives of A are described in Appendix B.

Implementation
All simulations and reconstructions are implemented in the Berkeley Advanced Reconstruction Toolbox (BART) using single-precision floating point arithmetic [31].The Bloch operator B is implemented in BARTs nonlinear operator framework [32].The calibrationless model-based reconstruction is based on an IRGNM-FISTA following Wang et al. [13] and we refer to this work for further details.
The Bloch operator includes a pixel-wise calculation of the signal evolution using STMs (Section 2.2) and of the partial derivatives with SAB (Section 2.1).ODEs are solved using the a Runge-Kutta algorithm (RK54) with adaptive step-sizes.
The error tolerances are chosen to be 10 −7 for the simulation comparisons in Section 2.1 and 2.2 as well as 10 −6 for further reduced computational costs in the Bloch model-based reconstructions.The Runge-Kutta solver exploits weights published by Dormand and Prince [33].For balancing the relative scaling of the partial derivatives during the optimization of Eq. 19 pre-conditioning following Wang et al. [13] is used.The initial wavelet regularization is set to α 0 =β 0 =1 and decreased by q=1/2 in each Newton iteration.The output of the Bloch model operator B is scaled according to Section D. As a globalized Newton method, the IRGNM does not require fine tuning of initial values.Here, the maps are initialized with the constants R 1 = 1 Hz, R 2 = 1 Hz, M 0 = 1, B 1 = 0 and the coil profiles are initialized with zero.
For comparison, we also implemented the reparameterized Look-Locker model from equation 29 in the same model-based reconstruction framework following Wang et al. [13].

Validation of Bloch Simulation
The accuracy of the SAB technique is validated with an IR bSSFP sequence for tissue with T 1 /T 2 = 1250/45 ms.An analytical solution for the IR bSSFP signal can be derived from the Bloch equations [34] assuming hard pulses, a perfect inversion and a perfect non-selective excitation.The symbolic partial derivatives with respect to the parameters R 1 , R 2 , M 0 and B 1 are calculated in Appendix F and are used as ground truth.For validation of the derivatives the ODE simulation parameters are chosen to be close to the assumptions of the analytical model.The ODE solution for the derivatives is also compared to the difference quotient techniques (DQ) calculated using two simulations M t,p and M t,p+h differing by a small perturbation h for all parameters p ∈ (T 1 , T 2 , M 0 , B 1 ) and time points t: The size of the perturbation is decreased until numerical noise dominates.
To validate the STM approach in the presence of RF pulses, gradients, and relaxation, a slice-selective excitation of a Hamming-windowed sinc-shaped inversion pulse with T RF = 1 ms, BWTP=1, ∆z =10 mm, G z = 10 mT/m is simulated.Relaxation parameters are selected based on typical human white matter values at 3 T: T 1 /T 2 =832/80 ms [35].The STM simulation is computed with the Runga-Kutta solver and the magnetization just before the slice rewinder is compared to a direct simulation exploiting the same solver.For comparison a traditional Bloch simulation technique based on temporal discretization with rotational matrices (ROT) is performed using a discretization rate of 1 MHz.According to the analysis shown in Figure S3, a sampling rate of 1 MHz is required for ROT to accurately model complex spin dynamics that include slice-selective RF pulses.The simulation speed is analyzed for a FLASH sequence with FA=8 • , TR/TE=3.1/1.7 ms and T RF =1 ms simulated for 101 isochromats homogeneously distributed along a slice of 0.02 m width and using a slice-selection gradient of 12 mT/m.The tissue parameters are set to relaxation times of white matter at 3 T.The simulations were executed on a single Intel(R) i7-8565U CPU core at 1.80 GHz.

Validation of Reconstruction
To validate the model-based reconstruction we further perform validations on numerical and experimental phantoms as well as in vivo data for both singleshot IR FLASH and IR bSSFP sequences with tiny golden-angle based radial sampling.The IR bSSFP sequence includes a prior α/2-TR/2 pulse to achieve a smooth signal evolution during the transient state [36,37].To be able to decouple the information of T 1 and T 2 for an IR bSSFP sequence (see section E) a B 1 map is acquired on the same slice using a vendor protocol based on rapid B 1 estimation with preconditioned RF pulses and Turbo-FLASH readout [38].All sequence parameters are shown in Table 1.
Phantom data for an IR FLASH sequence published by Wang et al. [39] was downloaded from Zenodo [40].This data was measured on a 3 T Magnetom Skyra by Siemens Healthcare (Erlangen, Germany) with a 20 channel head-coil.The measured phantom is a commercial reference phantom (Diagnostic Sonar LTD, Scotland, UK, Eurospin II, gel 3, 4, 7, 10, 14, and 16) consisting of six tubes with known T 1 relaxation values surrounded by water.A digital phantom of the dataset is created with the same sequence and acquisition characteristics as the downloaded measurement.The relaxation parameters were set to the estimated reference T 1 and T 2 values from previous studies [13,41].Additional phantom data for a radial single-shot IR bSSFP sequence was acquired on the T 2 spheres of a NIST phantom [42] on the same scanner and the same 20 channel head-coil.For comparison gold-standard maps of T 1 and T 2 are estimated on the same slice of the NIST phantom using fully-sampled single-echo spin-echo sequences.All phantom measurements were performed at constant room temperature (∼21 • C).The estimated quantitative values from gold-standard scans were used to simulate a second digital phantom with the same geometry as the measured T2-spheres of the NIST geometry (model 130).
Radial single-shot IR FLASH and IR bSSFP data for a brain of a healthy volunteer was acquired with a 20 channel head-coil on a Siemens Skyra 3T system (Siemens Healthcare, Erlangen, Germany) after obtaining written informed consent.The IR FLASH data was measured with a TR/TE of 4.1/2.58ms, a flip angle of 6 • , a band-width-time-product (BWTP) of 4 and a RF pulse duration of 1 ms.The field-of-view (FOV) was 220x220 mm 2 measured at 512 samples with two-fold oversampling.The single-shot IR bSSFP data was acquired with two different pulse durations: 1 ms and 2.5 ms respectively.The repetition and echo times were set to 4.88/2.44ms and 10.8/5.4 ms.A FOV of 200x200 mm 2 , a flip angle of 45 • , a BWTP of 4 and a base resolution of 256 was chosen and the same for both measurements.
The phantom and in vivo IR FLASH datasets are reconstructed with both the Look-Locker model and the Bloch model-based technique.As the IR FLASH sequence is insensitive to R 2 , its estimation in the Bloch model-based reconstruction was turned off by setting the scaling in the pre-conditioner to zero.The initial free decay of 15.3 ms during the non-selective hyperbolic secant inversion pulse is corrected in both cases using the correction published by Deichmann et al. [43].
The flexibility of the Bloch model-based reconstruction with its ability to also model more complex sequences is demonstrated using single-shot IR bSSFP.In the analysis of the measured dataset the B 1 map is used in the forward model for scaling the nominal flip angle for each pixel.Additionally, the B 1 estimation is turned off by setting the scaling to zero in the pre-conditioning.The method is numerically validated using a digital phantom of the T 2 sphere of the NIST phantom (model version 130) implemented in BART.The multicoil phantom is simulated in the frequency domain with the same sequence parameters as the measurement.The eight coils are compressed to four virtual coils using a singular-value decomposition (SVD).Complex Gaussian noise is added before coil compression to further avoid an inverse crime.To ensure realistic physical conditions the simulated signal model includes a non-selective hyperbolic secant inversion pulse and slice-selective excitations using multiple isochromats distributed over equally spaced slice-selection gradient positions.While analytical solutions of the Bloch equations require the assumption of perfect inversions and ideal non-selective excitation, the generic simulation of the Bloch model-based reconstruction technique can simulate more realistic signal models.We show how step-wise improvements to the model allow more accurate modelling of actual measurements.This is demonstrated on the numerical and measured NIST phantom datasets by performing reconstructions with various different assumptions about the signal model.The first reconstruction uses a model close to the analytical formula by assuming a perfect inversion and a nonselective excitation.We then add a realistic slice-selective excitation simulated as mean signal of various homogeneously spaced isochromats along the sliceselection gradient.To also model the effect of non-optimal inversion efficiency, the final reconstruction includes an extended model with realistic non-selective hyperbolic secant inversion.
The radial single-shot IR bSSFP in vivo data is reconstructed using the most realistic model.For comparison, the single-shot IR FLASH measurement was acquired on the same slice and reconstructed with the Bloch model-based reconstruction assuming a realistic IR FLASH signal model with non-selective hyperbolic secant inversion and gradient based slice-selective excitation model to estimate a T 1 map.
Details about the measurements can be found in Table 1.

Validation of Bloch Simulation
Figure 2 shows the partial derivatives for an IR bSSFP sequence with respect to R 1 , R 2 and B 1 for the analytical reference, the SAB technique and difference quotient (DQ) techniques with different perturbations h on the left.On the right Figure 2 presents the differences of DQ and SAB to the analytical reference.
As expected, the error of DQ decreases for small perturbations until numerical noise starts to dominate for very small h.
The SAB technique demonstrates a high accuracy and precision of estimating partial derivatives without requiring tuning of the perturbation level.
The error of the STM simulation is dominated by numerical noise due to limited floating point precision.With the parameters used here, the STM technique has substantially lower point-wise errors than ROT.
It demonstrates that STM reproduces the RK54 technique for finding solutions to the Bloch equations extremely well, while ROT is affected by errors due to the discretization with fixed sampling rate and its nature of being a first order method constrained to single floating point precision here.The runtime of RK54, STM and ROT is shown in Figure 3.B.The computational cost of ROT increases linear with higher sampling rates.The STM has higher initial costs than the other techniques which reflect the initial calculation of the state-transition matrices.The other methods are therefore faster for a small number of repetitions.For more repetitions, the STM becomes much faster as it requires only a few matrix multiplications per TR.A detailed comparison of the computational cost and accuracy of the RK54, STM and ROT techniques for various error tolerances and sampling rates can be found in Supplementary Section S3 and Supplementary Figure S3.S4 and Supplementary Figure S6, respectively.

Validation of Reconstruction
In the difference map between the T 1 maps of the two methods only the water background shows areas with minor differences.This probably results from small differences in the Sobolev regularization on the flip angle map in both techniques.At the walls of the inner tubes there is not enough signal and the T 1 maps are not well defined.Results for the radial single-shot in vivo IR FLASH data are shown in Figure 4.C.The parameter maps are visually indistinguishable except for minor artifacts in the areas of the head with flow related effects, which are not modelled by both signal models.The T 1 values in the marked ROIs for representative white and gray matter areas show a very good correspondence.The homogeneity within the white matter is high corresponding to a small standard deviation.
The complex M 0 parameter map in Figure 4.D reconstructed with the Bloch model-based reconstruction is of good quality showing no artifacts and a homogeneous phase.Only in border regions phase changes are present which are most likely caused by fat.The relative flip angle map is globally lower than one.It combines the effect of an imperfect slice-selective excitation and the present B 1 field.The intensity and phase of the estimated SVD-compressed virtual coil sensitivities is comparable to intensity and phase of sensitivities estimated with ESPIRiT [44](results not shown).
Figure 5 shows the reconstructed T 1 and T 2 maps of the digital NIST phantom (5.A) and the measurement (5.B) using the radial single-shot IR bSSFP acquisition for different models compared to the reference values in Bland-Altman plots.The most simplistic model assumes a perfect inversion and an ideal nonselective excitation (Perfect Inversion) and shows inaccuracies in the T 1 and T 2 estimation.By integrating a slice-selective excitation (Slice) the errors in T 2 are significantly reduced leaving an offset in T 1 .Adding a realistic hyperbolic secant inversion pulse to the forward model (Pulse+Slice) corrects for the T 1 offset leading to an accurate estimation of the relaxation parameters.These effects are present in both: the simulation and the measured data reconstructions.Important to note is that the NIST phantoms contains some spheres with extreme T 1 and T 2 parameters [42], which were excluded from the Bland-Altman analysis in Figure 5.A and 5.B for improved visualization.In particular, the three highest T 2 values (1.450 s, 0.388 s and 0.271 s) were removed for the reconstruction from simulated data and the highest and lowest T 2 values (1.450 s, 0.006 s) for the measured data.Especially the simple model has difficulties in finding the correct relaxation parameters in the reconstruction of the measured data, so that the mean value is outside the plotted region.A direct comparison of reference and estimated parameters in a diagonal plot can be found in Supplementary Figure S1.A Bland-Altman plot with all data points is shown in Supplementary Figure S2.
Reconstructed  2. The differences are smaller for the longer TR and longer RF pulse duration T RF compared to the short pulse protocol.This is likely due to the magnetization transfer effect (MT) that affects the IR bSSFP sequence but is not included in the current study.
The reconstruction times depend on the complexity of the forward model.The reconstruction of the IR FLASH datasets took about 80 s on a AMD EPYC 7662 64-Core CPU and a Nvidia A100-SXM-80GB GPU.The reconstruction of the simple forward model of the IR bSSFP NIST phantom dataset took 60 s, while the most complex model reconstruction took 38 min.The longest reconstruction times were required for the in vivo IR bSSFP dataset with short RF pulse.The strong slice-selection gradient prolonged the reconstruction to about 75 min.

Discussion
A nonlinear model-based reconstruction framework can be used in combination with well-crafted sequences and their analytical signal representations to accelerate quantitative MRI.This work presents a generalization of this well-known approach to arbitrary MRI sequences by exploiting the Bloch equations directly as forward model.This method becomes computationally feasible by including a direct sensitivity analysis of the Bloch equations.It allows us to use a generic ODE solver to compute the derivatives required by efficient nonlinear optimization algorithms such as the IRGNM.In comparison to techniques based on difference quotients it produces highly stable and accurate partial derivatives without the need of fine tuning perturbation levels.This was shown by estimating partial derivatives of an IR bSSFP experiment with the mentioned techniques and by comparing their results to the underlying analytical reference.
To further reduce computational demand, we exploit pre-computed STMs.They are used to solve the Bloch equations and the system describing their sensitivities simultaneously for arbitrary initial conditions for a given time span.They reduce the spin dynamics even in the presence of external fields, gradients, and relaxation to single matrix multiplications.It dramatically speeds up the reconstruction whenever the MRI sequence contains repeated patterns as it is often the case.In the presented example of the FLASH sequence with 101 isochromats along a slice and simulated for 1000 repetitions, the run-time of the simulation was reduced by a factor of 10 from 5 s down to 0.5 s in comparison to a regular Runge-Kutta ODE solver.Even in the presence of gradients, RF pulses and relaxation the slice-profile analysis showed the high accuracy of the STM technique in reproducing the ODE solver results.
Experimentally we confirmed that the Bloch model-based reconstruction reproduces the Look-Locker model as a special case.A comparison between both techniques showed only minor differences in the T 1 maps reconstructed from the single-shot phantoms and single-shot in vivo data.The integration of a generic Bloch simulation into the reconstruction adds the flexibility to analyze a broad variety of sequences.As an initial example, we applied the technique to IR bSSFP sequence and validated it using a numerical and measured NIST phantom dataset.By correctly modelling the slice-selective excitation and a non-selective hyperbolic secant inversion pulse highly accurate T 1 and T 2 maps could be obtained.
For a human brain, T 1 maps estimated from an IR bSSFP sequence were compared to Bloch model-based reconstructions of an IR FLASH acquisition of the same slice both including non-selective hyperbolic inversion and a sliceselective excitation.Here, differences could be observed which are likely caused by MT [45].This hypothesis is supported by the fact that prolonging the RF pulse duration and increasing the TR reduced the differences, but preliminary results (Supplementary Section S5) suggest that this does not explain the complete discrepancy and that other effects may also play a role.The NIST phantom measurement is not affected by MT effects, because it is based on water [42].The IR FLASH measurement is assumed to be unaffected by MT because of its small flip angle [46].
At this stage, the most relevant practical limitation is the need to manually tune the scaling factors used for pre-conditioning.For each analyzed sequence the relative scaling between the partial derivatives needs to be balanced manually to ensure smooth convergence.Future work is going to investigate automatic scaling techniques [47,48].

Conclusion
This work developed a generic framework for model-based reconstruction using the Bloch equations.The approach is validated numerically and tested experimentally using phantom and in-vivo scans. 1 Illustration of the operators used in the Bloch model-based reconstruction.The bottom part presents the ODEs for the signal and the derivatives.Figure 2 Left: Temporal evolution of the partial derivatives with respect to R 1 , R 2 and B 1 estimated for an IR bSSFP sequence with the SAB, DQ with varying perturbation levels h and the analytical references.Right: Plot with point-wise errors of the various DQ methods and the SAB with respect to the analytical reference.Note that the errors are presented in ppm for visualization.Figure 3 A: The slice-selection gradient based simulation for a Hammingwindowed sinc-shaped inversion pulse simulated with the RK54 framework is shown (left).The point-wise errors of the RK54 (top), the STM technique (center) and the ROT method (bottom) with sampling rate 1 MHz are plotted for the x-, y-and z-component of the magnetization.Note that the errors are scaled by large factors for visualization.B: The runtime of the STM technique is compared to the reference RK54 method and a ROT simulation performed with a sampling rate of 1 MHz.The simulation is performed for 101 isochromats homogeneously distributed along a slice-selection gradient during a FLASH sequence for various numbers of repetitions.The maximum of 1000 is chosen to cover about 4 s of acquisition, required to measure enough data points for mapping high T 1 values.A more detailed version of this figure has been added to the Supplementary Section S3.    Figure S5 Visualization of the magnetization transfer effect in T 1 maps reconstructed from an IR bSSFP sequence with varying TR and T RF .On the left the T 1 map corresponding to the longest T RF = 2.5 ms is plotted with colored ROIs.The mean value of the relaxation parameter values for the colored areas is plotted with its standard deviation on the right for various length of RFpulses T RF .The dotted lines represent the fitted values following the analysis of section S5. Figure S6 Bloch model-based reconstruction of the short RF pulse in vivo single-shot IR bSSFP dataset.The coil-sensitivities have been estimated with the regular method and have been added as fixed prior knowledge to both reconstruction.The top presents the regularized reconstruction, while the bottom results do not include any regularization on the parameter maps.Table 1 Table listing the sequence parameters for the performed measurements of this work.Table 2 Table listing the single-shot in-vivo IR bSSFP ROI analysis results presented in Figure 6.Table S1 Table listing the fitting parameters estimated for Figure S5.Table S2 Table listing the sequence parameters for the analysis in Figure S5.

A Combining Sensitivity Analysis with a Direct Sensitivity Analysis
The system matrix A(t) in equation 10 can be extended to include the sensitivity analysis for the three partial derivatives R 1 , R 2 and B 1 : with its corresponding parameter vector from equation 9:

B Forward Model Derivatives
The derivative of A in equation 19 follows by exploiting the Jacobi matrix and the product rule similar to [28,13].
The adjoint of the derivative becomes . . .

C Look-Locker Reparameterization
The Look-Locker model represents a special solution of the Bloch equations for an IR FLASH sequence [51].It assumes a perfect inversion, a small flip angle and short repetition times compared to the relaxation effects.Initial relaxation effects between the inversion and the first echo can be compensated analytically [43].The original formulation of the Look-Locker model with the parameters M 0 , M ss and These parameters are related to the underlying physical parameters M 0 , R 1 and α eff .With the assumption of short repetition times [52] the effective relaxation rate splits into R 1 and the readout relaxation rate R 1 determined by the effective flip angle α eff [53].The reparameterized Look-Locker model can be formulated as Here, Equation 29 depends directly on the physical parameters M 0 , R 1 and R 1 = − 1 TR ln cos α eff and allows adding prior knowledge about smooth B1 maps [54,43] to the reconstruction [53].This reparameterized model can be directly compared to a Bloch model-based reconstruction using the same set of parameters.

D Scaling Factors
The output of the forward operator B in equation 19 is the strength of the signal estimated by the Bloch simulation and its partial derivatives.The strength of the signal depends on the sequence and especially the applied flip angle.Therefore, the signals output differs for classical FLASH or bSSFP sequences influencing the weighting between the data fidelity and regularization terms changing the optimization behaviour.
A generic reconstruction aims for robustness against variations of the sequence parameters.It requires a scaling of the signal and its partial derivatives simulated within B. The implemented scaling is motivated by the Look-Locker model assumption that the longitudinal magnetization M z in equation 29 is proportional to the measured signal M xy and scaled to 1. Thus, the signal of a simulation with an initial magnetization of length 1 requires scaling of the simulated signal output M xy by the applied flip angle α and the relaxation effect during the echo time interval ∆t TE : Because of the short echo times of the sequences used in this work, the T 2 relaxation effect can be neglected.This assumption avoids additional T 2 dependencies of the estimated derivatives.The final scaling factor 1 sin α increases the robustness of the forward operator B in equation 19 to the choice of the applied flip angle in FLASH based sequences.For a bSSFP type sequence the flip angle in equation 30 needs to be halved to take its dynamics on the α/2 cone into account [55].

E IR bSSFP Information Encoding
The IR bSSFP signal behaviour is described by a limited exponential growth similar to the IR FLASH sequence [34].A single inversion recovery can encode information for estimating 3 parameters.While the IR FLASH sequence is sensitive to exactly 3 parameters, the IR bSSFP is also sensitive to T 2 , leading to 4 parameters in total.With a single limited exponential growth this additional parameter can not be encoded and two parameters need to be coupled.For bSSFP sequences the relaxation parameters T 1 and T 2 are coupled.Prior knowledge about B 1 can be used to decouple both relaxation parameters [34,56,57].

F Symbolic Derivatives of IR bSSFP
The analytical signal model for an IR bSSFP can be derived from the Bloch equations with the assumptions of hard RF pulses, a perfect inversion and an ideal α/2 − TR/2 magnetization preparation.
The signal is modelled by [34]: with This can be reparameterized to the physical parameters R 1 , R 2 , M 0 and α: with its symbolic derivatives:         For the on-resonance results the difference maps are small.Here, the reconstructed parameter maps for higher error tolerances are very similar to the reference map.For the signal model with simulation of the slice-selection gradient the reconstructed parameter maps are similar up to a tolerance of 0.001, while lower values result in large variations in the T 1 parameter map.The important result is, that the error of the simulation from an STM tolerance value of 0.001 still allows for accurate reconstructions even for complex spin dynamics.

S2 Supporting
In Figure S3.B the error plot from Figure 3 is extend to include more values for tolerance and sampling rate.To translate the required STM tolerance of 0.001 from Figure S3.A for complex dynamics to the required sampling rate for the ROT simulation, it demonstrates the point-wise errors (PWE) for RK54, STM and ROT for the same simulation with slice-selection gradient.The tolerance value for the estimated limit of STM of 0.001 produces an error similar to the 1 MHz sampling rate.This leads to the conclusion that at least a sampling rate of 1 MHz is required for accurate reconstruction of complex spin dynamics involving slice-selection gradients.
The simulation times for the same analysis with slice-selection gradient as presented in Figure 3

The
Bloch model-based reconstruction was compared to the Look-Locker modelbased version for simulated (Figure 4.A) and measured single-shot IR FLASH phantom data (Figure 4.B).Both methods recover high quality T 1 maps with small differences.Values for the same Regions-of-Interests (ROIs) are very simi-lar leading to their position on the diagonal of the Bloch vs. Look-Locker plot on the most right of Figure 4.A and Figure 4.B.The reconstructed tubes are very homogeneous in both reconstructions leading to low standard deviations.Reconstructions using different regularization parameters or no regularization are shown in Supplementary Figure T 1 and T 2 parameter maps for the two single-shot IR bSSFP scans of a human brain are shown in Figure 6.For comparison, a map with the Look-Locker model-based reconstruction of the IR FLASH scan of the same slice is added.Both IR bSSFP reconstructions show large offsets in T 1 compared to the IR FLASH reference.The relaxation values for the two analyzed ROIs are listed in Table

Figure 4
Reconstructed T 1 parameter maps for radial single-shot IR FLASH data acquired from a numerical (A) and measured phantom (B) as well as a human brain (C).The Bloch model-based reconstruction and the differences between the two methods are shown in the middle.The difference map is scaled up by a factor of 20 to improve visualization.On the right the T 1 values of the color-coded ROIs (arrows) of the Bloch reconstruction vs. Look-Locker reconstruction are plotted together with standard deviations.Besides the in vivo T 1 map presented in C, the Bloch model-based technique reconstructs a complex valued M 0 map, a relative flip angle map and complex coil sensitivities shown in D. Figure 5 A: Reconstructed T 1 and T 2 parameter maps and the corresponding ROI values for numerical radial single-shot IR bSSFP data of a digital multicoil reference object simulated in k-space.The left side shows the reconstructed parameter maps and the right the ROI analysis results in Bland-Altman plots relative to the simulated reference values.The analyzed ROIs are marked and numbered in the T 1 map.B: Reconstructed T 1 and T 2 parameter maps and the corresponding ROI values for a radial single-shot IR bSSFP measurement of the T 2 spheres of the NIST system phantom.The left side corresponds to the rightmost ROI analysis.The right side presents the comparison of the analyzed ROIs reconstructed with the Bloch model-based technique for various signal model assumptions compared with gold-standard reference values in Bland-Altman plots.For improved visualization individual outliers are removed from the plot.

Figure 6
The T 1 parameter map reconstructed from a radial single-shot IR FLASH in vivo dataset with a Look-Locker model-based reconstruction is shown on the left.It also shows the reconstructed T 1 parameter maps of a radial singleshot IR bSSFP in vivo dataset acquired on the same brain slice for short RF pulses (D RF : 1 ms, TR: 4.88 ms) on the top and long RF pulses (D RF : 2.5 ms, TR: 10.8 ms) on the bottom reconstructed with the Bloch model-based reconstruction.In the center column the difference maps are shown.The values corresponding to the colored ROIs are listed in Table 2. On the right the T 2 parameter maps for the short and long RF pulse experiments are shown reconstructed from the IR bSSFP sequence.Figure S1 A: Reconstructed T 1 and T 2 parameter maps and the corresponding ROI values for numerical radial single-shot IR bSSFP data of a digital multicoil reference object simulated in k-space.The left side shows the reconstructed parameter maps and the right the ROI analysis results relative to the simulated reference values wit standard derivations.The analyzed ROIs are marked and numbered in the T 1 map.B: Reconstructed T 1 and T 2 parameter maps and the corresponding ROI values for a radial single-shot IR bSSFP measurement of the T 2 spheres of the NIST system phantom.The left side corresponds to the rightmost ROI analysis.The right side presents the comparison of the analyzed ROIs reconstructed with the Bloch model-based technique for various signal model assumptions compared with gold-standard reference values.
Figure S2 A: Reconstructed T 1 and T 2 parameter maps and all corresponding ROI values for simulated radial single-shot IR bSSFP data.Here, the datapoints removed to improve visualization in Figure 5 are included.Reconstructed parameter maps are shown on the left Bland-Altman plots using the simulated reference values are shown on the right.The ROIs are marked and numbered in the T 1 map.B: Reconstructed T 1 and T 2 parameter maps and the corresponding ROI values for a radial single-shot IR bSSFP measurement of the T 2 spheres of the NIST system phantom.The left side corresponds to the rightmost ROI analysis.The right side shows the Bland-Altman plots comparing the values of the ROIs reconstructed with the Bloch model-based technique for various assumptions used in the signal model compared with gold-standard reference values.For improved visualization individual outliers were removed from the plot.FigureS3A: T 1 maps reconstructed from single-shot IR FLASH (same as for Figure4.B) for varying tolerance for computation of STMs.The forward model assumes on-resonant spins.The difference of the reconstructions for tolerances of 1e-5 to 1e-2 compared to the 1e-6 reference are shown in the lower row.The differences are scaled by larges factors for improved visualization.B: Same analysis as for part A, but assuming a complex forward model including a slice-selection gradient.C: Extension of Figure3.A with varying tolerances for RK54 and STM and sampling rates for ROT.D: Extension of Figure3.B with varying tolerances for RK54 and STM as sampling rates for ROT.

Figure S4 Reconstructed T 1
parameter maps similar to Figure 4.B for varying wavelet regularization strengths β 0 for the Bloch model-based reconstruc-tion (upper row).On the left a reference reconstruction estimated with the re-parameterized Look-Locker model-based technique is shown.In the bottom row the differences between Bloch model-based reconstruction and the reference map are shown scaled by a factor of 10 for improved visualization.

Figure 1 :
Figure 1: Illustration of the operators used in the Bloch model-based reconstruction.The bottom part presents the ODEs for the signal and the derivatives.

Figure 3 :
Figure 3: A: The slice-selection gradient based simulation for a Hammingwindowed sinc-shaped inversion pulse simulated with the RK54 framework is shown (left).The point-wise errors of the RK54 (top), the STM technique (center) and the ROT method (bottom) with sampling rate 1 MHz are plotted for the x-, y-and z-component of the magnetization.Note that the errors are scaled by large factors for visualization.B: The runtime of the STM technique is compared to the reference RK54 method and a ROT simulation performed with a sampling rate of 1 MHz.The simulation is performed for 101 isochromats homogeneously distributed along a slice-selection gradient during a FLASH sequence for various numbers of repetitions.The maximum of 1000 is chosen to cover about 4 s of acquisition, required to measure enough data points for mapping high T 1 values.A more detailed version of this figure has been added to the Supplementary Section S3.

Figure 4 :
Figure 4: Reconstructed T 1 parameter maps for radial single-shot IR FLASH data acquired from a numerical (A) and measured phantom (B) as well as a human brain (C).The Bloch model-based reconstruction and the differences between the two methods are shown in the middle.The difference map is scaled up by a factor of 20 to improve visualization.On the right the T 1 values of the color-coded ROIs (arrows) of the Bloch reconstruction vs. Look-Locker reconstruction are plotted together with standard deviations.Besides the in vivo T 1 map presented in C, the Bloch model-based technique reconstructs a complex valued M 0 map, a relative flip angle map and complex coil sensitivities shown in D.

Figure 5 :
Figure 5: A: Reconstructed T 1 and T 2 parameter maps and the corresponding ROI values for numerical radial single-shot IR bSSFP data of a digital multicoil reference object simulated in k-space.The left side shows the reconstructed parameter maps and the right the ROI analysis results in Bland-Altman plots relative to the simulated reference values.The analyzed ROIs are marked and numbered in the T 1 map.B: Reconstructed T 1 and T 2 parameter maps and the corresponding ROI values for a radial single-shot IR bSSFP measurement of the T 2 spheres of the NIST system phantom.The left side corresponds to the rightmost ROI analysis.The right side presents the comparison of the analyzed ROIs reconstructed with the Bloch model-based technique for various signal model assumptions compared with gold-standard reference values in Bland-Altman plots.For improved visualization individual outliers are removed from the plot.

Figure 6 :
Figure 6: The T 1 parameter map reconstructed from a radial single-shot IR FLASH in vivo dataset with a Look-Locker model-based reconstruction is shown on the left.It also shows the reconstructed T 1 parameter maps of a radial singleshot IR bSSFP in vivo dataset acquired on the same brain slice for short RF pulses (D RF : 1 ms, TR: 4.88 ms) on the top and long RF pulses (D RF : 2.5 ms, TR: 10.8 ms) on the bottom reconstructed with the Bloch model-based reconstruction.In the center column the difference maps are shown.The values corresponding to the colored ROIs are listed in Table 2. On the right the T 2 parameter maps for the short and long RF pulse experiments are shown reconstructed from the IR bSSFP sequence.

Figure S2 :
Figure S2: A: Reconstructed T 1 and T 2 parameter maps and all corresponding ROI values for simulated radial single-shot IR bSSFP data.Here, the datapoints removed to improve visualization in Figure 5 are included.Reconstructed parameter maps are shown on the left Bland-Altman plots using the simulated reference values are shown on the right.The ROIs are marked and numbered in the T 1 map.B: Reconstructed T 1 and T 2 parameter maps and the corresponding ROI values for a radial single-shot IR bSSFP measurement of the T 2 spheres of the NIST system phantom.The left side corresponds to the rightmost ROI analysis.The right side shows the Bland-Altman plots comparing the values of the ROIs reconstructed with the Bloch model-based technique for various assumptions used in the signal model compared with gold-standard reference values.For improved visualization individual outliers were removed from the plot.
are shown in S3.C for multiple solver tolerances and sampling rates.

Figure S3 :
Figure S3: A: T 1 maps reconstructed from single-shot IR FLASH (same as for Figure 4.B) for varying tolerance for computation of STMs.The forward model assumes on-resonant spins.The difference of the reconstructions for tolerances of 1e-5 to 1e-2 compared to the 1e-6 reference are shown in the lower row.The differences are scaled by large factors for improved visualization.B: Same analysis as for part A, but assuming a complex forward model including a sliceselection gradient.C: Extension of Figure 3.A with varying tolerances for RK54 and STM and sampling rates for ROT.D: Extension of Figure 3.B with varying tolerances for RK54 and STM as sampling rates for ROT.

Figure S4 :
Figure S4: Reconstructed T 1 parameter maps similar to Figure 4.B for varying wavelet regularization strengths β 0 for the Bloch model-based reconstruction (upper row).On the left a reference reconstruction estimated with the reparameterized Look-Locker model-based technique is shown.In the bottom row the differences between Bloch model-based reconstruction and the reference map are shown scaled by a factor of 10 for improved visualization.

Figure S5 :
Figure S5: Visualization of the magnetization transfer effect in T 1 maps reconstructed from an IR bSSFP sequence with varying TR and T RF .On the left the T 1 map corresponding to the longest T RF = 2.5 ms is plotted with colored ROIs.The mean value of the relaxation parameter values for the colored areas is plotted with its standard deviation on the right for various length of RFpulses T RF .The dotted lines represent the fitted values following the analysis of section S5.

Table 1 :
Table listing the sequence parameters for the performed measurements of this work.