A unifying view on extended phase graphs and Bloch simulations for quantitative MRI

Quantitative MRI methods and learning-based algorithms require exact forward simulations. One critical factor to correctly describe magnetization dynamics is the effect of slice-selective RF pulses. While contemporary simulation techniques correctly capture their influence, they only provide final magnetization distributions, require to be run for each parameter set separately, and make it hard to derive general theoretical conclusions and to generate a fundamental understanding of echo formation in the presence of slice-profile effects. This work aims to provide a mathematically exact framework, which is equally intuitive as extended phase graphs (EPGs), but also considers slice-profiles through their natural spatial representation. We show, through an analytical, hybrid Bloch-EPG formalism, that the spatially-resolved EPG approach allows to exactly predict the signal dependency on off-resonance, spoiling moment, microscopic dephasing, and echo time. We also demonstrate that our formalism allows to use the same phase graph to simulate both gradient-spoiled and balanced SSFP-based MR sequences. We present a derivation of the formalism and identify the connection to existing methods, i.e. slice-selective Bloch, slice-selective EPG, and the partitioned EPG. As a use case, the proposed hybrid Bloch-EPG framework is applied to MR Fingerprinting.

www.nature.com/scientificreports/ necessitates the introduction of surrogate states to describe the effect of the RF pulse. These states are no longer separated by the spoiling moment, but by an arbitrary discretization moment arising from the time resolution of the RF pulse simulation. The EPG branches at every discretization point leading to a plethora of configuration pathways, which are no longer easily analyzed in view of echo formation. Thus, the elegant yet simplistic picture of configuration states shifted by spoiling gradients and their interference through RF pulses has been greatly complicated.
In 1994, Sobol and Gauntt 26 described in detail the magnetization response to gradient-echo sequences with constant repetition time and gradient spoiling moment. They identified interference of substates (later referred to as configuration states) to be accountable for artifacts in gradient-echo sequences. While they worked with isolated, regular prisms (voxels), they noted that slice profiles are rarely rectangular and thus precise, a priori calculation of their effect on the interference of substates was intractable. Despite this early work, considerable misunderstanding of spoiled states in gradient echo sequences has persisted in the community (e.g. the notion of "killer gradients" destroying transverse magnetization) and effects such as off-resonance artifacts in FISP-MRF came as a surprise (in defense of all authors cited and fellow researchers who fell for this: they also came as a surprise to us) 18,22,23 .
In this work, we present a contemporary solution to Sobol's and Gauntt's approach for sequences employing arbitrary flip angle trains, which is derived from the rotation operator algorithm that forms the basis of current EPG implementations 8,12,18 . We will work as far as possible with the nomenclature by Weigel 8 and Malik et al. 12 . Where applicable, crosslinks to the work of Sobol and Gauntt are included 26 .
Our formalism is a hybrid between a conventional spatial Bloch simulation and an extended phase graph: it represents spoiling through configuration states and slice profiles in their natural spatial basis. This allows us to retain the elegance of the EPG formalism which enables the straight-forward identification of echo formation, while keeping the intuition of a spatial basis for the description of slice profiles. Moreover, the formalism is easily implemented in contemporary EPG codes and provides analytical access to multiple key signal dependencies: We show that the signal s (n) for each time-point n in an arbitrary gradient-echo sequence with constant TR and fixed spoiling moment k sp is given by a weighted sum over spatially-resolved configurations F +(n) k (z) . Here, k is the configuration order and z the spatial position along the slice profile. The signal is given by where C(z) is the coil sensitivity, M 0 (z) is the equilibrium magnetization and W k (ω, R ′ 2 , TE) is an analytical weighting function that captures dependencies on off-resonance ω and microscopic dephasing R ′ 2 for arbitrary echo times TE (see the "Theory" section for the derivation). Through the exponential weighting factor e ikk sp z , the signal's dependency on the spoiling moment k sp is also obtained analytically. An approximation of the spatially resolved configuration states F +(n) k (z) can be readily obtained by piecewise discretization of the slice profile and independent EPG simulations (spatially resolved EPG).
We identify balanced steady-state free-precession (bSSFP or True-FISP) sequences as a special case of spoiled (FISP) SSFP by considering k sp = 0 . For spoiled SSFP, we recover the partitioned EPG approximation by the limit k sp → ∞ . Omitting spatial integration, the hybrid Bloch-EPG yields results equal to a direct Bloch simulation. Additional Fourier transformation allows to recover the slice-selective EPG solution. Thus, the formalism shows the connection between contemporary solution approaches and sequence variants in one unified framework. As an example, we apply the formalism to the raised problem of off-resonance artifacts in gradient-spoiled MR Fingerprinting (FISP-MRF).

Theory
The solution of the homogenized form of the Bloch equations can be obtained by consecutive application of four-dimensional matrix operators (Rotation Operator Algorithm, ROA) 7 , which act on magnetization vectors of the form where M ± = M x ± iM y denotes the transverse magnetization, M z the longitudinal magnetization, and M 0 the equilibrium magnetization and r is the spatial position with r := x, y, z T . Following the nomenclature of Weigel 8 , hard pulses will be denoted by T ϕ (α) with flip angle α and phase ϕ . Relaxation and recovery for time t will be written as R(t) implicitly including the dependency on relaxation parameters T 1 and T 2 . Off-resonance ω and gradients G(t) lead to phase accrual at position r of φ = ω · t + γ dtG(t) · r , which will be denoted by the operator S(φ) and where γ is the gyromagnetic ratio ( γ1 H ≈ 2π · 42.577MHz/T).
In the theoretical excursion of this work, we will work with four-dimensional representations of the operators, which are presented in Appendix A. For the numerical implementation, however, it is more efficient to work with the original three-dimensional matrices, since the fourth dimension-taking care of equilibrium magnetization-remains unchanged and only needs to be considered for the recovery of longitudinal magnetization to equilibrium.
Effective RF pulse matrix and the hard pulse approximation. Slice selection is performed using amplitude-modulated (AM) RF pulses (soft pulse) played concurrently to a constant slice-selection gradient G www.nature.com/scientificreports/ ( Fig. 1a) 27 . The through-slice direction is, without loss of generality, defined to be given by z , parallel to the sliceselect gradient. Pulse duration T ex , slice selection gradient strength G , and the slice thickness δ are connected by the time-bandwidth product κ of the pulse through the relationship 27 Each RF pulse is preceded by a pre-winder gradient with moment k pre and followed by a refocusing gradient of moment k post . Since the Bloch equations are linear, an effective matrix operator RF(z) can be defined 28 .
In this work, we restrict the analysis to short pulses ( T ex ≪ T 1 , T 2 , T ′ 2 ) with sufficiently high time-bandwidth product κ ( 2πκ ≫ T ex · ω ), such that relaxation and precession during the pulse can be replaced by relaxation and precession operators acting after the pulse. Consequently, the matrix operator acts instantaneously at the "focus point" t f of the shaped pulse 29 . The focus point can be obtained through a maximization of the freeinduction decay (FID) and is connected to the rewinder area by k post = γ G·t f 30 . The pulse simulation is performed using the hard-pulse approximation (see Fig. 1a,b) 31 . The RF amplitude modulation B + 1 (t) is discretised into hard-pulses of T ϕ (α i ) , each acting instantaneously at equidistant time points The effective RF pulse matrix is obtained by alternating hard pulses and phase accrual Here, the product ( ) denotes an ordered matrix product such that earlier time-points are found to the right of later time-points. As relaxation is neglected, the effective RF pulse matrix RF(z) is a rotation matrix in each z ( detRF(z) = 1) . It can be written in the form of an effective hard pulse T ϕ(z) (α(z)) and an additional effective phase accrual term S φ pre (z) α(z) is the effective flip angle profile and ϕ(z) the effective pulse phase, which together determine the excitation slice profile of the RF pulse. For echo and storage components, the additional phase-term φ pre (z) needs to be considered (Fig. 1c). This representation is effectively decomposing the RF pulse matrix into Euler angles at each z . The resulting intrinsic rotations are first about the z-axis with angle φ pre (z) − ϕ(z) , second, the tipping (5) RF(z) = T ϕ(z) (α(z))S φ pre (z) .

Figure 1.
Slice-selective RF pulse and its numerical representation. (a) Slice-selection gradient and amplitudemodulated RF pulse for slice-selective excitation shown together with prephasing and refocusing gradients. Same-colored dashed and filled gradient areas are of equal size. (b) Graphical representation of the hard pulse approximation assuming short RF pulses compared to rate of change of the B1 + envelope. (c) Effective RF pulse representation using the effective slice profile α(z) and pulse phase ϕ(z) as well as an effective prephaser φ pre (z). www.nature.com/scientificreports/ of the magnetization by the flip angle α(z) -a rotation about the rotated x ′ axis, and finally, a back rotation of ϕ(z) about the z″-axis.
Sequence description: MR fingerprinting and configuration state imaging. Throughout this work, we assume constant off-resonance ω , repetition time TR and spoiling gradients with moment k sp = e z · k sp in the direction of slice-selection e z and an arbitrary set of RF pulses, here given solely by a flip angle variation ( Fig. 2a) realized by pulse amplitude scaling. The fundamental building block of MR Fingerprinting sequences are either balanced or gradient-spoiled SSFP-type sequence modules (Fig. 2b,c) 2,3 . Each is represented by a sequence of operators: slice-selective RF pulse RF n (z) with nominal flip angle α n , where n enumerates consecutive TR intervals, acquisition AQ , relaxation R(TR) and phase accrual S ω · TR + k sp · z 2,3 . Note, that balanced SSFP is obtained by considering the special case k sp = 0.
The magnetization M (n) (z) directly after the nth RF pulse is then given by where M (0) (z) = M 0 (z)(0, 0, 1, 1) T is the initial magnetization. For simplicity, we will assume acquisition to happen instantaneously at the echo time TE relative to the focus point of the soft pulse and assume readout gradients to be zeroth-moment nulled, i.e. magnetization remains unchanged by acquisition. At TE , magnetization is given by For spoiled SSFP, configurations other than the FID can be refocused, which is the basis for e.g. double-and triple-echo steady-state (DESS and TESS) relaxometry 32,33 , time-reversed FISP (PSIF or T2-weighted FFE), and diffusion-weighted SSFP sequences 34,35 . Acquisition is preceded by a refocusing gradient with a zeroth moment of integer multiples of k sp . To keep the total spoiling moment constant, compensation gradients with inverse gradient sign are inserted after the readout (Fig. 2c). The magnetization of the qth echo is given by Hybrid Bloch-EPG formalism. Magnetization dynamics. The aim is to reformulate the recurrence Eq. (6) to obtain an analytical expression in the phase accrual φ := ω · TR + k sp · z per TR . This can be achieved, as shown in Appendix B, by a Fourier series expansion of the phase accrual operator (6) M (n) (z) = RF n (z) · R(TR) · S ω · TR + k sp · z M (n−1) (z), Through the transformation, the dependency of the magnetization M (n) (z) on the spoiling moment k sp and the off-resonance ω is now fully analytical and given by a Fourier series. The series coefficients x  The is a convenient way to write the shifting operator of the EPG in analytical form. The first element of the diagonal matrix corresponds to up-shifting of F + and the second to down-shifting of F − configurations. k is referred to as the order of the configuration x (n) k (z) . Here, spatial variation in slice-direction is introduced in the SR-EPG only through the effective RF pulse matrix RF n (z) , however, variations could also be introduced through spatially-dependent relaxation parameters and transmitfield inhomogenieties 12 .
Signal reception. Using the coil sensitivity in slice direction C(z) , which collates all dependencies on coil geometry and induction physics into one convenient sensitivity factor, the signal received directly after the n th RF pulse can be written as The solution ansatz for M (n) (z) (Eq. 10) can be inserted and simplified leading to Following Leupold 36 , the signal Eq. (13) can further be extended to include signal decay due to Lorentziandistributed, microscopic inhomogeneities in the magnetic field distribution leading to an effective decay rate R ′ 2 , as well as, signal decay and phase accrual to an arbitrary echo time-point 0 ≤ TE ≤ TR . This leads us to the hybrid Bloch-EPG equation stated in the introduction Thus, the acquired signal s (n) is determined by summation over configuration orders of the Fourier transformation dze ikk sp ·z · · · of weighted spatially-resolved configurations W k ω, R ′ 2 , TE F . The weighting function W k provides the analytical dependency on off-resonance ω , R ′ 2 and the TE , while the Fourier transformation allows for arbitrary spoiling moments.
Note: Substituting k sp → k sp + γ �B 0 (z i ) , the formalism equivalently applies in locally first-order to inhomogeneities in the static magnetic field �B 0 (z i ) , i.e. local field gradients 26 .
We termed the formalism the hybrid Bloch-EPG since it solves the interplay of slice profiles, spoiling gradients, and off-resonance by combining the natural representation of both: (1) extended phase graphs for spoiling and off-resonance and (2) the spatial domain representation of slice profiles as in conventional Bloch simulations.
Discretization and spatially-resolved extended phase graphs. For computer simulations, we discretise the spatial domain in slice direction into N equidistant bins centered at position z i with width z ( i = 1 . . . N ). We define an interpolation kernel function I(z i , z) , which allows us to approximate the RF pulse matrix by For simplicity, we assume the matrix to be piecewise constant within each bin leading to the following kernel function identifying the discrete spatially-resolved configuration states F k,i , and assuming spatially constant coil sensitivity C and equilibrium magnetization M 0 , we arrive at the approximate solution to the hybrid Bloch-EPG which can be directly evaluated from a spatially-resolved extended phase graph simulation. The term zsinc k k sp z 2 in the kernel's Fourier transformation F k,i k sp is equivalent to the "damping factor" D l introduced in Sobol and Gauntt 26 , and accounts for dephasing within each bin (voxel). Note: the equivalent formulation derived by Sobol and Gauntt was stated for a single voxel, i.e. without the summation over the slice profile denoted by i · · · here.
Configuration state imaging. For refocusing of an arbitrary configuration q-in the following we will call it the refocused echo q to distinguish it from configuration state F + k -we can modify Eq. On the other hand, Fourier transformation of the transverse magnetization allows to recover the solution of the slice-selective EPG (ssEPG), which directly yields the magnetization's k-space representation by performing the hard pulse approximation using the EPG (here, to distinguish configuration order k from k-space, k ′ is used to denote spatial frequencies) These relationships are graphically depicted in Fig. 3. A spatially-resolved EPG is shown in the top left corner for the configuration order k = −6 . . . 6 as a function of the position z along the slice profile.
Bloch simulation (Fig. 3, left to right). Analytical weighting with the weighting kernel W k ω, R ′ 2 , TE as well as resolving the spoiling and echo refocusing by applying e i(k−q)k sp ·z to the data is employed to resolve the intermediate representation of the configuration states. These are then superimposed to obtain the slice profile through the hybrid Bloch-EPG formulation (black, right), which is equivalent to direct Bloch simulation (blue). The root-mean-square error of the complex slice profiles normalized to the signal amplitude obtained via Bloch simulation (nRMSE) is 0.71% (see "Methods", Eq. (40)). By integrating over the magnetization profile, signal s (n) is obtained. The predicted MR signal deviates by 3.8 · 10 −4 % between hybrid Bloch-EPG and Bloch simulation.
Note: A change in off-resonance (here, ω = 0 is shown) will lead to different relative phases of the weighted configurations leading to a phase shift of the oscillations in the intermediate representation. This will modify the slice profile by interference and lead to variations in the gradient-echo signal, which has been observed and described for MRF previously 18,22,23 .
Slice-selective EPG (Fig. 3, top to bottom). Similarly, we arrive at the intermediate k-space representation ( k ′ ) through spatial Fourier transformation, where configuration states appear to be shifted by the spoiling moment www.nature.com/scientificreports/ k sp . Here, off-resonance ω will modify the complex phase of each configuration. Summing over the configuration order k leads to the slice-selective EPG representation in the bottom row, where the weighted configurations will interfere, which was previously described by Sobol and Gauntt 26 . The nRMSE of the k-spaces of hybrid Bloch-EPG and ssEPG is 0.15% (see "Methods", Eq. (41)). Here, the MR signal is given by the DC component ( k ′ = 0 ) as known from conventional EPGs 8,12,18 . The predicted signal of the ssEPG deviates from the Bloch simulation by 2 · 10 −2 % . Refocusing of higher configuration orders leads to shifting of the state picture, centering a different peak at k ′ = 0 and the formation of the q th refocused echo.
Special cases of spoiling: balanced SSFP and partitioned EPG. Two special cases of spoiling can be discussed, k sp = 0 equivalent to balanced SSFP; and k sp → ∞ , which will be shown to be equivalent to the partitioned EPG solution.
k sp = 0-balanced SSFP. Balanced SSFP can be considered a special case of spoiled sequences with vanishing spoiling moment ( k sp = 0 ). In this case, the hybrid Bloch-EPG Eq. (20) can be simplified to . Through summation over state orders k , interference between configuration states leads to modification of the slice profile and hence a change in observed signal. A single SR-EPG simulation provides analytical access to four dimensions of the solution space: TE, ω , R ′ 2 , and the spoiling moment k sp . While ssEPG and Bloch simulations can be obtained from each other, the mapping of SR-EPG to Bloch or ssEPG cannot be reversed as the configurations cannot be disentangled from the spatial information alone. The partitioned EPG (pEPG) approximation is given by summation over the F0 state in the SR-EPG (top, left). Solid and broken lines denote imaginary-and real-part, respectively. www.nature.com/scientificreports/ Hence, in balanced SSFP, the configuration states only contain information about off-resonance, which is by definition independent of the spatial coordinate. Thus, slice profiles can directly be integrated over before weighting and summation. k sp → ∞-partitioned EPG. In the limit k sp → ∞ , the sinc-function in Eq. (18) can be simplified to 37 leading Eq. (20) to become Hence in the limit of infinite spoiling, the slice profile integration can be directly performed and only one state, the q-th configuration will contribute. Interference of adjacent configurations, as predicted by the hybrid Bloch-EPG, does not take place. This implies that no magnitude or phase variation other than the trivial phase accrual ω · (qTR + TE) can be observed 18 . This approximation is equivalent to the approach by Lebel and McPhee, which was termed partitioned EPG (pEPG) by Ostenson 13,14,18 .  18 , an offresonance dependency remains, which is sensitive to the spoiling moment. (3) For infinite spoiling, a single state will contribute the signal, which is equivalent to the solution of the partitioned EPG (pEPG). Here, banding is fully suppressed and no off-resonance dependency apart from the trivial phase accrual at TE is observed. For the simulation of the off-resonance dependency, the employed weighting function is marked and the accompanying parameters are given in bold. www.nature.com/scientificreports/ Figure 4 visualizes the signal formation within the hybrid Bloch-EPG from one SR-EPG (left) for the three possible cases of spoiling: no spoiling ( k sp = 0 , top row), finite spoiling ( k sp = 0 , middle row) and strong spoiling ( k sp → ∞ , bottom row). After spatial Fourier transformation, a configuration order k vs. refocused echo q plot can be obtained (left column). For any given q , configurations along the white, horizontal lines are selected. With k sp = 0 , the echo order q is irrelevant as no spoiling moment is applied, i.e. it leads to degeneracy for all echo orders. For k sp → ∞ (bottom row), only one configuration remains, which leads to the pEPG solution. Configurations are weighted by W k (top), which contributes the analytical dependency on TE , R ′ 2 and ω . Summation over the weighted configurations allows to recover, e.g. the signal dependency on off-resonance, here shown for TE = 0.
To further the understanding of the three cases of spoiling in SSFP, Fig. 5 shows the possible gradient echo sequences in a combined configuration k vs. k-space ( k ′ ) graph. The plot on the right is achieved by spatially Fourier transforming the SR-EPG on the left. During signal acquisition, a superposition of weighted configuration states k W k F + k is observed. Each configuration state is evaluated at a different spatial frequency k ′ , which can be found by selecting them along straight lines.
Balanced SSFP corresponds to selecting spatial frequencies of configurations along horizontal lines (green), i.e. all configurations are evaluated at the same spatial frequency. Angulated lines (red) depict spoiled SSFP sequences. The slope is given by the spoiling moment k sp . Asymmetry in the slice profile phase (see Fig. 11) arising from the non-linear response of the Bloch equation's solution with respect to RF fields leads to an asymmetric k′-space, which is seen in Fig. 5 by the asymmetry between the upper and the lower half-space. Hence, negative spoiling moments lead to a different spoiled SSFP signal than positive. Vertical shifting of lines (broken) corresponds to adding gradients before the acquisition object with respective compensation gradients after the readout, i.e. phase encoding gradients. Small shifts ( ≪ k sp ) are attributed to phase encoding as used in Fourier imaging (green, broken). Shifts of integer multiples of k sp refocus different configuration orders, e.g. F −1 or the PSIF-echo (orange, broken).

Application of the hybrid Bloch-EPG to MR fingerprinting and configuration state imaging.
Off-resonance. In Fig. 6, the signal is plotted as a function of time index and for two different relaxation parameters for both balanced and spoiled SSFP sequences. The respective off-resonance value is color coded, where positive and negative off-resonance share the same color bar as they are indistinguishable in the magni- For spoiled SSFP, the signal phase is shown relative to the phase accrual observed within the pEPG formalism, i.e. accounting for sign change due to inversion recovery and the trivial accrual of ω(qTR + TE) due to the refocusing of different configuration orders q . Spoiled SSFP ( k sp = 4 · 2π/δ , δ : slice thickness) shows greatly reduced off-resonance dependency compared to bSSFP. Higher configuration orders experience larger magnitude and phase variation, which can be attributed to neighbouring, lower configurations being of larger amplitude, thus leading to stronger interference.  www.nature.com/scientificreports/ tude plots. The spoiling moment of k sp = 4 · 2π/δ reduces the off-resonance dependency greatly compared to bSSFP for all displayed configuration orders. While magnitude variation persists equally for all configurations, phase variation is increased for higher-order configurations. Phase variation is predominantly found in the region of low signal magnitude as well as for high flip angles. Fig. 7, the signal for time-point 90 in the MR fingerprinting sequence is shown for both balanced and spoiled SSFP sequences and four configuration orders. As in conventional balanced SSFP, the bSSFP-MRF signal behaves like a gradient-echo for short TE (1), and spin-echo like for long TE (2) 38 . Configurations of higher order ( F 1 and F −2 ) cannot be observed due to the rapid transverse decay for R2′ = 1/10 ms. Constructive (3) and destructive (4) interference of configurations leads to the banding structure observed in balanced SSFP.

R2′, off-resonance and echo-time dependency. In
Spoiling moment. In Fig. 8, the signal response for time-point 90 is plotted against off-resonance and spoiling moment for the four configuration orders. At k sp = 0 , all configurations show the same balanced SSFP response. For small deviations of k sp from zero, the different configuration states can be thought of as phase encodings of the balanced SSFP slice profiles. Asymmetry in the effective slice profile leads to different response for positive and negative spoiling moments. For sufficiently large spoiling moments, both magnitude and non-linear phase variations can be effectively suppressed.

Suppression of magnitude variation in MRF.
Signal magnitude and phase vary with spoiling moment and fingerprinting index. Figure 9 shows the standard deviation of the signal magnitude with respect to off-resonance normalized to the average MRF signal as a function of spoiling moment k sp and time n , where �· · · � ω denotes the average over ω . Periodic minima are observed, which are encountered at approximately 2π/δ . In the time-point region between 15 and 50 (flip angle ramp between 15° and 50°), the location of the minima varies. Thus, no spoiling moment can be found that reduces the off-resonance signal variation below 1% (black) for all time points simultaneously. For gradientrecalled echo (GRE, F 0 state of gradient-spoiled SSFP) and time-reversed GRE ( F −1 state of gradient-spoiled SSFP), a spoiling moment of at least 8 · 2π/δ is necessary to reduce magnitude variation below 1%. in Eq. (18) in the hybrid Bloch-EPG (yellow, dotted). This demonstrates the importance of the damping factor to resolve the effect of spoiling on a sub-discretization level. Aliasing occurs as soon as the damping factor becomes zero for the highest occupied state k max i.e. a state k max is populated, which can no longer be resolved by the slice profile discretization alone. The respective configuration state at which aliasing occurs is given by and evaluates to 250, 25 and 12.6 for the shown case, which is in agreement with the findings shown in Fig. 10. Figure 10. Comparison of signal-time behavior for different slice profile discretizations and simulation methods. The hybrid Bloch-EPG with reduced number of slice profile discretization points N sp (solid, gray) approximates well the signal obtained with the Hybrid Bloch-EPG and N sp = 3001 (black). Omitting the damping factor in the hybrid Bloch-EPG signal Eq. (19) leads to aliasing (yellow, dotted), which is equivalent to the results of a Bloch simulation (violet, broken).

Figure 9.
Off-resonance signal variation as a function of time point and spoiling moment. standard deviation of the signal magnitude with respect to off-resonance normalized to the mean signal amplitude plotted against spoiling moment and temporal index. Between index 15 and 50, minima of the signal variation do not follow a horizontal line, thus, no single spoiling moment can be found that reduces slice profile induced off-resonance variation to below 1% in this exemplary case. Only with strong spoiling ( ≥ 8 · 2π/δ ) can the variation consistently be suppressed for the F0 configuration. This demonstrates that the pEPG is ill-suited to predict MRF signal time-courses for standard spoiled SSFP-based MRF and even more so in the context of configuration state imaging. www.nature.com/scientificreports/

Discussion
We have presented a formalism that employs the spatially-resolved EPG as an input to obtain the MR signal of spoiled and balanced SSFP sequences with variable, slice-selective RF pulses. The formalism allows to capture the dependency on spoiling moment k sp , echo time TE , microscopic dephasing R2 ′ , and off-resonance ω analytically. We recovered three known solution methods to the Bloch equations, i.e. spatial-domain simulations (Bloch simulation), slice-selective extended phase graphs (ssEPG), and the partitioned EPG (pEPG) and presented them in one unified picture.
Interpretation of the hybrid Bloch-EPG. Considering gradient-spoiled SSFP, spoiling gradients store magnetization history in high spatial frequencies of the magnetization distribution. The present formalism suggests that off-resonance artifacts in spoiled SSFP sequences arise from overlapping past and current magnetization in k-space leading to interference of configuration states. This occurs whenever a state's k-space representation has finite intensity at multiples of the spoiling moment. Since each configuration order k shows different off-resonance phase ω(kTR + TE) , interference of the magnetization states leads to signal magnitude and phase variations. Changes in the spoiling moment can either reduce or amplify interference. Generally, the overlap of the configurations is reduced with increasing spoiling moment 26 .
Through the analytical dependency on the spoiling moment, the special case k sp = 0 can be considered, which corresponds to balanced SSFP (True-FISP). To the knowledge of the authors, this is the first time that a common formalism for balanced and spoiled SSFP-based MRF sequences is presented that allows to seamlessly and analytically sweep between these sequences. This also applies when fully neglecting any spatial variation in the magnetization apart from spoiling (pEPG). In this case, interference of configuration states does not occur for k sp = 0 . However, setting k sp = 0 , the balanced SSFP solution can again be obtained 36 . Since SR-EPG calculations are readily implemented and already widely used, adoption of the hybrid Bloch-EPG formula is straightforward.
Fundamental to the slice-selective EPG formalism is a constant spoiling moment increment k , which is given by the time-resolution t employed in the hard pulse approximation (Eq. 4) In order to ensure linear phase graph growth, spoiling moment k sp , refocusing gradient moment k post and prewinder k pre need to be integer multiples of k . This can either be enforced through rounding, rendering the solution inexact or by designing the pulse sequence respectively. Contrary, in Bloch simulations, the minimal number of isochromats is determined by the spoiling moment and number of sequence repetitions simulated to avoid aliasing artifacts. In the hybrid Bloch-EPG, however, the mathematical separation of slice profile and spoiling means that aliasing cannot occur, and the slice profile resolution can be chosen according to accuracy needs (Fig. 10). Here, the damping factor plays a crucial role by analytically resolving the intra-voxel phase dispersion on the sub-discretization level. Likewise, the number of configuration states solely depends on the number of time points simulated for the MRF sequence and is fully independent of the time-resolution of the RF pulse. This flexibility and intrinsic numerical stability of the method comes at the expense of increased memory requirements, as the full SR-EPG needs to be retained in order to evaluate the analytical dependencies. The implementations of Bloch, ssEPG and hybrid Bloch-EPG presented here were not designed with computational efficiency in mind. An exact computational complexity analysis and performance comparison of the techniques in the realm of MRF dictionary generation is beyond the scope of this work.
Extensions to the hybrid Bloch-EPG. Off-resonance distribution. In Eq. (14), we followed Leupold 36 , to extend our model to include dephasing due to microscopic field inhomogeneities effectively assuming a Lorentzian off-resonance distribution. Using the work of Ganter 39 , other distributions can equally be considered by integrating the signal s (n) over an off-resonance distribution p(ω) . Since ω only enters into s (n) via the exponential term e iω(k·TR+TE) contained in the weighting function W k , we can readily define a generalized weighting function where p(t) is the Fourier transformation of the off-resonance distribution function p(ω) For example, by substituting the Fourier transformations provided by Ganter and keeping a global offset frequency ω , we arrive at the following two expressions for the weighting function for Gaussian-distributed inhomogeneities and spherically distributed inhomogeneities with J 1 being the Bessel function of the first kind 37 . www.nature.com/scientificreports/ Spectral information. Analogously, the model can be extended to include multi-peak spectra analytically, e.g. for accurate modelling of fat using a 7-peak model 40 . Assuming peak frequencies of ω j , line-width R ′ 2,j , and peak weights w j , the weighting function can be replaced by assuming equal relaxation properties (T1 & T2) for each peak. If peaks have different relaxation properties, a separate phase graph calculation must be performed for each peak before summation.
Variable echo time and configuration state imaging. Our derivation does not impose restrictions on the type of readout employed in each TR. Thus, sequences with variable echo times and multi-echo sequences can equally be discussed 2,41 . This also includes sequences refocusing multiple configuration states, e.g. double-and tripleecho steady-state (DESS & TESS) 32,33 or MRF employing configuration state imaging.
Arbitrary spoiling and k-space readout. We restricted spoiling to the through-slice direction, which allowed us to discuss slice profiles as a source for off-resonance artifacts in spoiled SSFP-based MRF. Using the same approach, however, we can also discuss arbitrary spoiling directions. In this case, three factors need to be considered separately: (1) the spatial distribution of the configuration states F + k (r) , which is a consequence of RF pulses, B1 + inhomogeneity, and variation of relaxation properties; (2) the coil sensitivity-weighted equilibrium magnetization distribution C(r)M 0 (r) ; and (3) Fourier encoding of the spatial frequency k ′ with the associated time point of encoding τ (k ′ ), which depends on the chosen acquisition trajectory. The signal Eq. (14) can then be generalized to Employing the convolution theorem, the signal is given by Thus, interference of configurations can be equally formulated as a question of how well separated states are in k-space 42 , an insight already formulated and discussed in detail by Sobol and Gauntt 26 . In addition, effects of non-instantaneous acquisition, such as blurring through off-resonance and k-space filtering as a result of signal decay are equally captured through the dependency of the model on the sampling time-point τ (k ′ ).
Variable TR sequences. Central to the derivation is a constant phase increment φ = ω · TR + k sp · z . This directly enforces a constant repetition-time TR, off-resonance ω , and spoiling moment k sp . Only in this case, the phase graphs for off-resonance and spoiling are equal. In addition, this leads to a constant phase advance angle φ per TR and the phase graph will grow linearly with each iteration. Dropping this requirement, exponential growth ( ∼ 3 n ) of pathways quickly renders the SR-EPG intractable 8 . Dropping the constant TR requirement, a phase graph needs to be simulated for each off-resonance frequency separately.
Equilibrium magnetization and coil sensitivity. Contrary to RF pulses, equilibrium magnetization and coil sensitivities enter as factors into the generalized signal equation (Eq. 35). From a state-mixing perspective, they can be represented as convolutions in k-space (Eq. 36). Hence, they generally lead to broadening of the configuration's k-space representation and lead to increased mixing. Considering coil sensitivities as the only spatially varying component, they are expected to show minor state interference since they typically vary slowly in space. Contrary, equilibrium magnetization can vary strongly, e.g. at tissue interfaces or due to contrast agent uptake, leading to increased interference of configurations and thus off-resonance artifacts.
Continuous broadening of configuration states. Considering a single RF pulse within the small tip-angle approximation ( α ≪ 1rad ), Pauly's k-space interpretation of slice-selective excitation can be employed to estimate the spatial frequencies excited 43 . The maximal frequency component in the excitation slice profile is then given by half the gradient moment of the slice selection gradient γ GT ex /2 . Thus, the minimal spoiling moment to fully dephase the transverse magnetization after a single RF pulse is given by Since RF pulses can be represented by independent matrix products in each spatial position, they also transform to convolutions in k-space. Hence, consecutive application of multiple RF pulses will generally lead to a broadening of the configuration state's k-space populating at most the spatial frequency nπκ/δ after the n th RF pulse. This hinders complete spoiling unless signal cancellation due to interference is achieved (Figs. 7 and 8).
B1 + inhomogeneity. Transmit field inhomogeneity ( B + 1 ) affects the RF pulse matrices and thus is equally expected to lead to continuous configuration state broadening despite being slowly varying in space.
Implications for MR fingerprinting and configuration state imaging. A detailed discussion of the effect of slice profile and spoiling in the context of off-resonance artifacts has already been presented by Ostenson et al. 18 We will therefore only discuss the specific findings in the context of configuration interference and configuration state imaging. There are two regions in the test example MRF sequence, where magnitude and phase variation due to off-resonance are strongest: (1) the time at which longitudinal magnetization passes through zero during the recovery process (for F(0) and T1 = 800 ms, T2 = 50 ms (red curve), time index 25…40 in Fig. 6) and (2) high flip angle regions (time index 50…100). Following the argument of Sobol and our findings regarding interference of configuration states, an increase in off-resonance variation is related to increased contributions from unwanted states F k =q relative to the refocused state F q .
In the case of configuration state imaging, a time-reversed spoiled SSFP-based fingerprinting sequence reading out the F −1 configuration shows similar behavior to the F 0 imaging, albeit with generally reduced signal magnitude. For higher configurations ( F +1 and F −2 ), configuration state interference leads to pronounced signal magnitude variations (see Figs. 6 and 9). We attribute this to the generally low intensity of the F +1 and F −2 states compared to their neighboring F 0 and F −1 states. Of note, while favorable spoiling moments exist, the valleys of low magnitude variations become steeper in Fig. 9, rendering them more susceptible to changes in the spoiling moment, e.g. due to local field inhomogeneities.
To the knowledge of the authors, MRF combined with configuration state imaging has so far not been published apart from time-reversed gradient-spoiled MRF (PSIF) 44 . While the combination of the two seems straightforward, the increase in off-resonance variation for higher order states as well as the difficulty of finding a common spoiling moment that works equally well across configurations might be one obstacle hindering successful implementation. Especially time-reversed spoiled-SSFP might be interesting in the realm of diffusionsensitive MR Fingerprinting 34,35 .
To the knowledge of the authors, current MRF implementations employ RF pulse amplitude scaling to implement variable flip angles. This typically leads to an amplification of side lobes at high flip angles ( Fig. 11(1)) and non-linear phase accrual over the slice profile ( Fig. 11 (2) and (3)), which in turn result in increased off-resonance variation in high flip angle regions. Optimal design of RF pulses, pre-winder, and refocusing gradients might allow control over configuration interference and thus off-resonance artifacts in gradient-spoiled MRF.
The effective pulse representation (Eq. 5) allows to fully describe slice selective RF pulses through a positiondependent phase accrual φ pre (z) and a hard pulse. Hence, together with the hybrid Bloch-EPG, we can formulate a model for slice profile corrections in DESS & TESS 45 sequences using the analytical solutions for the transverse configuration states F + k in steady state by Hänicke (use Eq. (21) in Hänicke, where S n denotes the steady state solutions of the signal magnitude of the configuration orders n) 46 . By defining we can evaluate the signal of the q th refocused echo in the presence of an arbitrary slice profile by Finally, the hybrid Bloch-EPG allows to recover four dimensions analytically: TE , ω , R ′ 2 , and the spoiling moment k sp . Optimal parameter combinations of R ′ 2 and ω , could be determined by conventional fitting routines rather than being resolved through discrete atoms in a dictionary. This could help to reduce the dimensionality of the dictionaries and improve model-data consistency. Analytical access to k sp further allows for fast evaluation of optimal spoiling moments to improve sequence design and mitigate slice profile effects. www.nature.com/scientificreports/

Conclusion
The hybrid Bloch-EPG formalism proposed in this work utilizes a numerically calculated spatially-resolved EPG to analytically recover the MR signal dependency on echo time, microscopic dephasing, off-resonance, and spoiling moment. It is applicable to sequences with varying RF pulses but constant TR and spoiling moment, as they are commonly used in MR Fingerprinting. The hybrid Bloch-EPG formalism allows to seamlessly tune between spoiled and balanced SSFP-based MRF sequences. It was shown that off-resonance artifacts arise from interference of configuration states, which is a direct effect of incomplete spoiling due to non-rectangular slice profiles. The formalism retains a clear separation of configuration states and slice profiles, which allows for a more fundamental understanding of echo formation and artifact generation than alternative Bloch or sliceselective EPG calculations.

Methods
The hybrid Bloch-EPG formalism was implemented in MATLAB 2020b (Mathworks, Natik, MA, USA) and compared to implementations of both Bloch simulation and slice-selective EPG. Hann-apodized sinc pulses with flip angles α n according to Fig. 2a were generated with κ = 10 and nominal slice thickness of 5 mm by amplitude scaling. Refocussing gradient area k post was calculated once for the maximal flip angle of 70° and the prewinder was defined as k pre = − γ GT ex + k post . For each n = 1 . . . 100 , effective RF pulse matrices RF n (z) were obtained using the hard pulse approximation. An adiabatic inversion pulse preceded the test sequence, Figure 11. Effective RF pulse decomposition. Two exemplary slice profiles for (a) 30° and (b) 70° RF pulse decomposed into effective flip angle α(z) (black) and rf pulse phase profiles ϕ(z) (blue, solid). The broken line corresponds to the effective pre-phaser φ pre (z) , which needs to be considered for storage and echo components of the pulse. (1) At high flip angle, slide lobes are amplified due to violation of the small tip angle approximation.
(2) At low flip angles, the effective pulse phase is not flat as the gradient waveforms for the sinc pulses were optimized for the maximum flip angle (3) as is performed on contemporary scanner systems. For large flip angles, pronounced non-linearity in both phase terms arises at the boundary of the slice profiles. www.nature.com/scientificreports/ which was simulated as an ideal pulse inverting the longitudinal magnetization only. Repetition time was assumed to be constant with TR = 15 ms. Bloch simulation was performed using the rotation operator algorithm and propagation of the magnetization through the hard pulse approximation (Eq. 4) in each repetition. The slice-selective EPG was implemented according to Ostenson et al. by performing the hard pulse approximation in the Fourier domain 18 . For comparison of the techniques in Fig. 3, rounding of spoiling, refocusing and prephasing moments to the state increment of the slice-selective EPG ( �k = γ G�t ) was performed as well as choosing 5119 isochromats and 145 temporal steps to ensure equivalence of the methods. For visualization purposes only, RF pulse phase alternation was not simulated.
For all other results and if not otherwise stated, plots depict phase graphs and signals after 91 RF pulses, i.e. time-point 90 within the fingerprinting train (Fig. 2a) with simulation parameters of T1 = 800 ms , T2 = 50 ms , R ′ 2 = 0 , TE = 0 , and ω = 0 . The spatially-resolved EPG was simulated with 3001 isochromats, and a field-of-view of three times the nominal slice thickness.
For the comparison of hybrid Bloch-EPG magnetization profiles and the Bloch simulation (see Fig. 3), the normalized root-mean-square error (nRMSE) was calculated as follows and the nRMSE for the ssEPG according to Code to reproduce the figures in this paper can be found in the supporting material.

Data availability
All data presented was generated through custom code, which is shared as part of this publication.