Optimized Quantification of Spin Relaxation Times in the Hybrid State

TL;DR

This paper develops optimized spin ensemble trajectories for relaxometry in the hybrid state, demonstrating improved efficiency over traditional methods through numerical optimization and in vivo validation.

Contribution

It introduces a novel approach to optimize spin trajectories for relaxometry, combining theoretical analysis, numerical optimization, and experimental validation.

Findings

Hybrid-state sequences outperform traditional relaxometry methods.

Optimized trajectories achieve near-optimal signal-to-noise efficiency.

Joint T1 and T2 encoding in the hybrid state is more efficient than sequential methods.

Abstract

Purpose: The analysis of optimized spin ensemble trajectories for relaxometry in the hybrid state. Methods: First, we constructed visual representations to elucidate the differential equation that governs spin dynamics in hybrid state. Subsequently, numerical optimizations were performed to find spin ensemble trajectories that minimize the Cram\'er-Rao bound for $T_1$-encoding, $T_2$-encoding, and their weighted sum, respectively, followed by a comparison of the Cram\'er-Rao bounds obtained with our optimized spin-trajectories, as well as Look-Locker and multi-spin-echo methods. Finally, we experimentally tested our optimized spin trajectories with in vivo scans of the human brain. Results: After a nonrecurring inversion segment on the southern hemisphere of the Bloch sphere, all optimized spin trajectories pursue repetitive loops on the northern half of the sphere in which the beginning of the first and the end of the last loop deviate from the others. The numerical results obtained in this work align well with intuitive insights gleaned directly from the governing equation. Our results suggest that hybrid-state sequences outperform traditional methods. Moreover, hybrid-state sequences that balance $T_1$- and $T_2$-encoding still result in near optimal signal-to-noise efficiency. Thus, the second parameter can be encoded at virtually no extra cost. Conclusion: We provide insights regarding the optimal encoding processes of spin relaxation times in order to guide the design of robust and efficient pulse sequences. We find that joint acquisitions of $T_1$ and $T_2$ in the hybrid state are substantially more efficient than sequential encoding techniques.

Equations

$$\dot{r}(t)=-\left(\frac{{\cos }^2\vartheta (t)}{T_1}+\frac{{\sin }^2\vartheta (t)}{T_2}\right)r(t)+\frac{\cos \vartheta (t)}{T_1}$$

$${\sin }^2\vartheta =\frac{{\sin }^2\frac{\alpha }{2}}{{\sin }^2\frac{\varphi }{2}\cdot {\cos }^2\frac{\alpha }{2}+{\sin }^2\frac{\alpha }{2}},$$

$${\mathbf{b}}_1$$

$${\mathbf{b}}_2$$

$${\mathbf{b}}_3$$

$$rCRB(T_1)$$

$$rCRB(T_2)$$

$${\sigma }_{T_{1,2}}\ge \sqrt{rCRB(T_{1,2})\frac{T_R}{T_{\text{exp}}}}\cdot \sigma \cdot T_{1,2}.$$

Full text

Methods: First, we constructed visual representations to elucidate the differential equation that governs spin dynamics in hybrid state. Subsequently, numerical optimizations were performed to find spin ensemble trajectories that minimize the Cramér-Rao bound for $T_{1}$-encoding, $T_{2}$-encoding, and their weighted sum, respectively, followed by a comparison of the Cramér-Rao bounds obtained with our optimized spin-trajectories, as well as Look-Locker and multi-spin-echo methods. Finally, we experimentally tested our optimized spin trajectories with in vivo scans of the human brain.

Results: After a nonrecurring inversion segment on the southern hemisphere of the Bloch sphere, all optimized spin trajectories pursue repetitive loops on the northern half of the sphere in which the beginning of the first and the end of the last loop deviate from the others. The numerical results obtained in this work align well with intuitive insights gleaned directly from the governing equation. Our results suggest that hybrid-state sequences outperform traditional methods. Moreover, hybrid-state sequences that balance $T_{1}$- and $T_{2}$-encoding still result in near optimal signal-to-noise efficiency. Thus, the second parameter can be encoded at virtually no extra cost.

Conclusion: We provide insights regarding the optimal encoding processes of spin relaxation times in order to guide the design of robust and efficient pulse sequences. We find that joint acquisitions of $T_{1}$ and $T_{2}$ in the hybrid state are substantially more efficient than sequential encoding techniques.

\jnlcitation\cname

, , , and , (\cyear2018), \ctitleOptimized Quantification of Spin Relaxation Times in the Hybrid State, \cjournalPre-Print, \cvol.

††articletype: Preprint - Submitted to Magnetic Resonance in Medicine††footnotetext: Word Count: 3591

1 Introduction

The dynamics of large spin-1/2 ensembles in a magnetic field are commonly described by the Bloch equations, which capture the macroscopic effects of spin-lattice and spin-spin interactions with the characteristic time constants $T_{1}$ and $T_{2}$, respectively. Robust and rapid quantification of these parameters is important for more objective diagnoses, longitudinal studies, and computer aided diagnosis. The exponential behavior of $T_{1}$ relaxation makes a series of inversion-recovery experiments a natural choice for quantitative $T_{1}$ mapping. Similarly, a series of spin-echo experiments with varying echo times provides an intuitive way to estimate $T_{2}$. Both methods effectively strive to trace exponential relaxation curves following a single perturbation from thermal equilibrium. However, $T_{1}$ and $T_{2}$ can also be estimated by measuring magnetization in various steady-state conditions 1, 2, 3, 4, or even in more complex conditions that depart from the steady state5. Magnetic Resonance Fingerprinting (MRF)5, for example, opened a door for flexible and efficient parameter mapping with higher signal to noise ratio (SNR) efficiency than was possible with more simplistic spin trajectories. The large number of design options for MRF pulse sequences has, however, limited optimization of MRF experiments to heuristic approaches 5, 6, 7, 8, 9, 10. Furthermore, the lack of an intuitive understanding of the spin dynamics in MRF, and the over-simplified nature of the Bloch models employed in typical MRF approaches, creates a risk of introducing biases in the estimated relaxation times 9. Recently, we analyzed biases introduced by $B_{1}$- and $B_{0}$-inhomogeneities, including intra-voxel dephasing, which is also known as inhomogeneous broadening, and found a sequence design space that combines the steady state’s robustness with the efficiency of the transient state 11. During experiments designed in this space, the spin ensemble establishes a so-called hybrid state, whose dynamics can be described by a closed-form solution of the Bloch equations. Here, we provide an in-depth analysis of hybrid-state sequences and discuss how relaxation times can be encoded efficiently with this approach. By comparing hybrid-state sequences that are optimized to encode either $T_{1}$ or $T_{2}$ alone to a joint optimization, we show that the second parameter can be measured at virtually no extra cost.

2 Methods

2.1 Hybrid-State Spin Dynamics

We start by describing the spin dynamics in spherical coordinates, which are here defined by $x=r\sin{\vartheta}\cos\varphi$, $y=r\sin{\vartheta}\sin\varphi$ and $z=r\cos{\vartheta}$, where $r$ is the radius, $\vartheta$ the polar angle or the angle between the magnetization and the $z$-axis and $\varphi$ the azimuth or the angle between the $x$-axis and the projection of the magnetization onto the $x$-$y$-plane. For practical purposes, we use the limits $-1\leq r\leq 1$, $0\leq\vartheta\leq\pi/2$, and $0\leq\varphi<2\pi$ to uniquely identify the polar coordinates. Recapitulating Ref. 11, the entire spin dynamics in hybrid state is captured by the radius $r$ of the magnetization, i.e. its magnitude combined with a sign, which is controlled only by the polar angle $\vartheta$:

[TABLE]

Here, $\dot{r}(t)$ denotes the derivative of $r$ with respect to time.

The flip angle and the phase of the radio frequency (RF) pulses, as well as the repetition time $T_{R}$ have only an indirect and joint effect on the spin dynamics:

[TABLE]

where $\phi=\omega_{z}T_{R}$ describes the phase accumulated during one repetition time and $\alpha$ denotes the flip angle. This equation reduces to $\vartheta=\alpha/2$ for $\phi=\pi$, which we define as the on-resonance condition. In practice, $\phi=\pi$ is assigned to the on-resonant spin isochromat by the common phase increment of $\pi$ in consecutive RF pulses.

2.2 Numerical Optimizations

Various aspects of the MRF pulse sequence design can be improved. A complete analysis considering all design parameters simultaneously is currently out of reach. Instead, this work focuses on the analysis and optimization of the $T_{1}$- and $T_{2}$-encoding power of the RF-pulse train in hybrid-state pulse sequences with balanced gradient moments, independent from the specific k-space trajectory that may be used in the imaging sequence. To this end, the Cramér-Rao bound 12, 13 is used to provide a universal limit for the noise variance of the estimated parameters. Given an unbiased estimator, i.e. the fitting algorithm used to calculate the proton density ($PD$), $T_{1}$ and $T_{2}$, the noise variance of these parameters is at least as big as the corresponding Cramér-Rao lower bound. This very general and established metric has been utilized for optimizing MR parameter mapping experiments in Refs. 14, 15, 16 amongst others, and for MRF in particular in Ref. 17.

The hybrid state allows us to approximate a voxel’s signal at the echo time $T_{E}=T_{R}/2$ by a single isochromat 11, similar to balanced steady-state sequences 18. With $x(t)=r(t)\sin\vartheta(t)$, we can calculate a discretized signal vector $\mathbf{x}\in\mathbb{R}^{1\times N_{t}}$ from Eqs. (1) and (2). This vector describes the observed signal at $N_{t}$ time points, and is used to calculate the Fisher information matrix $\mathbf{F}\in\mathbb{R}^{3\times 3}$ whose entries $\mathbf{F}_{i,j}=\mathbf{b}_{i}^{T}\mathbf{b}_{j}/\sigma^{2}$ are given by

[TABLE]

Here, $\sigma^{2}$ describes the input variance. The vectors $\mathbf{b}_{i}$ describe the derivatives of the signal evolution with respect to all considered parameters. By normalizing with the input variance $\sigma^{2}$, the duration of the experiment $T_{\text{exp}}=N_{t}T_{R}$, the repetition time $T_{R}$, and the squared relaxation time, we can define the relative Cramér-Rao bounds to be

[TABLE]

The normalization by the variance cancels out the variance in the definition of the Fisher information matrix, the normalization by the squared relaxation time is performed in order to best reflect the $T_{1,2}$-to-noise ratio (defined as $T_{1,2}/\sigma_{T_{1,2}}$), and the normalization with $T_{\text{exp}}/T_{R}$ makes the $rCRB$ invariant to the number of measurements so that it describes the noise efficiency per unit time. In absolute numbers, the noise in the resulting parameter estimations is given by

[TABLE]

Note that the simulations published in Ref. 11 indicate that the actual noise in the estimated parameters is close to this bound.

Given that $\vartheta$ is the effective drive of the spin dynamics 11, we optimized $\vartheta$ directly and picked the sequence parameter ($\alpha,T_{R},\phi$) retrospectively (cf. Eq. (2)). Eq. (1) is an uncoupled first order differential equation and can be solved for different boundary conditions 11. Here, we focus on inversion-recovery balanced hybrid-state free precession (IR-bHSFP) sequences, i.e. hybrid-state sequences that depart from thermal equilibrium by the application of an inversion pulse, which is accounted for by the boundary condition $r(0)=-1$.

We used a Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm 19 with $rCRB(T_{1})$, $rCRB(T_{2})$ and $rCRB(T_{1})+rCRB(T_{2})$ as objective functions. The numerical optimization was based on $\vartheta(nT_{R})$ with the repetition time of $T_{R}=4.5~{}\text{ms}$ and with $n\in\{1,2,\ldots,N_{t}\}$. The derivatives in Eqs. (3)-(5) were calculated analytically, as was the gradient of the objective function with respect to each $\vartheta(nT_{R})$.

Since the $rCRB$ intrinsically compares a fingerprint to its surrounding in the parameter space, only a single set of relaxation times was used for the optimization. In particular, we used the relaxation times $T_{1}={781}~{}\text{ms}$ and $T_{2}={65}~{}\text{ms}$, which correspond to values reported for white matter6. All optimizations were initialized with the pattern provided in Ref. 9.

We performed the optimizations of IR-bHSFP sequences that exploit the full quadrant of the Bloch sphere ($0\leq{\vartheta}\leq\pi/2$), and repeated the same optimizations with $0\leq{\vartheta}\leq\pi/4$ in order to limit the flip angle to $\alpha\leq\pi/2$, ensuring consistent slice profiles by virtue of the linearity in the small tip-angle approximation 20, and aiding compliance with safety considerations by avoiding high power large flip-angle pulses. For comparison, we also performed unconstrained optimizations of Look-Locker 21 and multi-spin-echo experiments.

2.3 In Vivo Experiments

An asymptomatic volunteer’s brain was imaged following written informed consent, and according to a protocol approved by our institutional review board. Measurements were performed with the IR-bHSFP sequences that minimize $rCRB(T_{1})$, $rCRB(T_{2})$ and $rCRB(T_{1})+rCRB(T_{2})$, limited to $0\leq{\vartheta}\leq\pi/4$. All experiments were performed on a 3T Skyra scanner (Siemens, Erlangen, Germany). The 16 head elements of the manufacturer’s 20 channel head/neck coil were used for signal reception.

At the beginning of the sequence, a secant inversion pulse 22 with a duration of 10.24ms was applied, followed by a spoiler gradient. The other RF-pulses were implemented as slice-selective sinc pulses with a time-bandwidth product of 2 and a duration of 1.10ms.

Spatial encoding was performed with a radial trajectory and a golden angle increment 23. The spatial resolution of the maps is $1~{}\text{mm}\times 1~{}\text{mm}\times 3~{}\text{mm}$ at a FOV of $256~{}\text{mm}\times 256~{}\text{mm}\times 3~{}\text{mm}$. The readout dwell time was set to $2.4~{}\upmu\text{s}$ and an oversampling factor of 2 was applied. The total scan time of each sequence was approximately 3.8s.

The raw data were compressed to 8 virtual receive coils via SVD compression. Thereafter, image reconstruction was performed with the low rank alternating direction method of multipliers (ADMM) approach proposed in Ref. 24. The employed dictionaries were computed with Eqs. (1) and (2) and covered the range between 300ms and 6s in logarithmic steps of 2%. The dictionaries covered the range of $T_{2}$ values between 10ms and 3s in logarithmic steps of 2%. Slice profile correction was incorporated by adopting Ref. 25. Off-resonance and $B_{1}$ correction could be applied on the basis of separately acquired maps. However, we did not correct for these effects here in order to avoid complicating the noise analysis. The dictionary was compressed to include singular vectors corresponding to the 6 largest singular values resulting from a singular value decomposition. The data consistency step of the ADMM alogrithm was performed with 20 conjugate gradient steps. In order to prevent non-linear effects from complicating the noise assessment, only a single ADMM iteration was performed and no spatial regularization was applied.

3 Results

3.1 Visualization of the Differential Equation

In order to provide some intuition for the hybrid-state spin dynamics, we visualize its governing differential equation (Eq. (1)) in Fig. 1a. One can identify the so-called steady-state ellipse 26, 27, 28, 29, which separates areas of the Bloch sphere in which the magnetization grows (red) or shrinks (blue). The absolute value of the derivative of Eq. (1) with respect to $T_{1}$ is in the upper limit twice as high on the southern hemisphere than on the northern hemisphere (note the asymmetric color bar in Fig. 1b), which indicates that $T_{1}$ can be encoded faster on the southern hemisphere. A second benefit of inversion recovery experiments becomes evident when comparing the derivative with respect to $T_{1}$ (b) to the one with respect to $T_{2}$ (c). One changes sign between the hemispheres while the other does not. Thus, encoding data on the southern and northern hemisphere minimizes the correlation between the derivatives, which should yield a superior SNR in the estimated relaxation times by virtue of a lower Cramér-Rao bound (cf. Eqs. 4-5).

When considering $T_{1}$ only, Fig. 1b suggests that there is benefit in keeping the magnetization close to the $z$-axis on the southern hemisphere, and close to the origin on the northern hemisphere, since the absolute value of the derivative is largest in those areas. The derivative with respect to $T_{2}$, on the other hand, reaches its maximum far away from the $z$-axis (c). Starting from thermal equilibrium, this suggests that it makes sense to acquire signal while the magnetization relaxes along the equator. After the magnetization is brought to zero, the behavior of the derivative induces the magnetization to grow to a comparably large $r$ in order to drive the magnetization back into areas with a large derivative with respect to $T_{2}$.

3.2 Spin Dynamics on the Bloch Sphere

Fig. 2 depicts the hybrid-state spin ensemble trajectories and the corresponding driving functions $\vartheta(t)$ that were numerically optimized for $T_{1}$ and/or $T_{2}$ encoding. Some features are shared by all patterns, while others differ depending on the figure of merit. In the paragraphs that follow, the common features are described first, followed by the optimization-specific features.

All patterns are comparatively smooth, which is not explicitly enforced by the optimization. They start with a nonrecurring inversion segment on the southern hemisphere, followed by repetitive loops on the northern hemisphere. The beginning of the first loop differs slightly from that of the other loops, which can be explained by the initial conditions of the loop. At the time when the magnetization crosses the origin (turquoise segments in Fig. 2), the derivatives with respect to the relaxation times are comparatively large. Furthermore, the derivative with respect to $T_{2}$ has not yet changed its sign. These two properties make this segment a valuable asset to help minimize the correlation between the derivatives with respect to the two relaxation times. As a consequence, the optimized trajectories have a comparatively large polar angle in this segment. In contrast, the subsequent loops start with a segment in which the magnetization recovers along the $z$-axis, reflecting small and correlated derivatives.

Another common feature of all trajectories is that the polar angle decreases rapidly when the shrinking magnetization reaches the steady-state ellipse (magnifications in Fig. 2a,d,g,j,m,p). The arrival at the steady-state ellipse (blue ellipse in Fig. 2) also concludes the last loop in the $T_{2}$-specific and in the jointly optimized patterns so that they take advantage of the $T_{2}$-dominated shrinkage of the magnetization (d,g,m,p). By contrast, the $T_{1}$-specific trajectories (a,j) conclude by maximizing $r$ and $dr/dT_{1}$ and then briefly bringing the prepared magnetization close to the equator in order to maximize the signal.

Taking a closer look at the $T_{1}$-optimized hybrid-state pattern, we observe that the spin trajectory stays close to the $z$-axis on the southern hemisphere and close the origin on the northern hemisphere (Fig. 2a), which is exactly the behavior we expect from intuitive understanding of the derivatives (Section 3.1). In fact, the numerical optimization suggests that it is beneficial to forgo any signal at the beginning of the sequence (dark blue segment) in favor of an information-rich signal thereafter. Comparing the unconstrained optimization ($0\leq\vartheta\leq\pi/2$) with the constrained one ($0\leq\vartheta\leq\pi/4$), we find that the spin trajectory remains virtually unaffected on the southern hemisphere (Fig. 2j-l). On the northern hemisphere, however, the unconstrained pattern exploits the maximum polar angle of $\pi/2$ to quickly kill the magnetization and create a sort of saturation-recovery loop. When the polar angle is constrained, on the other hand, the magnetization shrinks more slowly and, consequently, fewer loops can be performed in the same amount of time.

When a hybrid-state sequence is optimized purely for $T_{2}$ encoding without constraining the polar angle, the magnetization first decays along the equator (Fig. 2d), which is in agreement with the intuitive understanding described in Section 3.1. As $r$ approaches zero, its derivatives with respect to the relaxation times become small as well and the spins follow a trajectory along the $z$-axis (bright blue segment). Thereafter, the spins follows comparatively large loops as predicted. Surprisingly, the optimization does not result in further segments with $\vartheta=\pi/2$. Instead, the trajectory approaches this value only as the magnetization approaches the origin. Limiting the polar angle to $0\leq\vartheta\leq\pi/4$ distorts the spin trajectory and the maximum value $\vartheta=\pi/4$ is exploited during several segments (m,n).

After the inversion pulse, the trajectory of the combined optimization oscillates between large and small $\vartheta$-values. The corresponding derivatives of the magnetization have peaks that are shifted with respect to one another (Fig. 2g-i). The size of the loops lies in between the ones of the $T_{1}$- and $T_{2}$-optimized trajectories. When the polar angle is limited to $0\leq{\vartheta}\leq\pi/4$, the spin trajectory is distorted in a similar way as for the $T_{2}$-optimized case.

The hybrid-state equation of spin motion (Eq. (1)) can be written in a dimensionless form (Eq. (LABEL:eq:Bloch_r_dimensionless) in the supporting material), which highlights the fact that only the ratios $T_{\text{exp}}/T_{2}$ and $T_{1}/T_{2}$ of the involved time constants affect the hybrid-state spin trajectories. Varying $T_{\text{exp}}/T_{2}$ in the joint optimization results in different trade-offs between the inversion phase and the repetitive loops (see supporting Fig. LABEL:fig:Bloch_sphere_vary_Texp). Short sequences ($T_{\text{exp}}\approx T_{2}$) require a $T_{2}$-optimized inversion phase similar to Fig. 2m since the temporal constraints allow only for small, $T_{1}$-dominated loops. Increasing the sequence duration allows for larger loops so that the optimized inversion phase increasingly trades $T_{2}$ for $T_{1}$ encoding (cf. Fig. 2p). Long sequences ($T_{\text{exp}}\gg T_{2}$) use a $T_{1}$-optimized inversion-recovery segment (cf. Fig. 2j) followed by a combination of large loops that predominantly encode $T_{2}$ and smaller loops, which provide additional $T_{1}$ encoding.

Supporting Fig. LABEL:fig:Bloch_sphere_vary_T1 analyzes the effect of variations of $T_{1}/T_{2}$ on the spin ensemble trajectory and demonstrates that the optimized loop size decreases with an increasing ratio. Large $T_{1}$ values limit how much the magnetization can grow within a given time. Therefore, the optimized $\vartheta$-pattern spends most of the time for this recovery and the resulting loop size is small. One can further observe that the optimization of $rCRB(T_{1})+rCRB(T_{2})$ results in different trade-offs between $T_{1}$ and $T_{2}$ encoding at different $T_{1}/T_{2}$ ratios (Fig. LABEL:fig:Bloch_sphere_vary_T1c). At intermediate ratios, the algorithm slightly favors the $T_{1}$-encoding, which is reflected both in the Cramér-Rao bound and in the corresponding spin ensemble trajectory, as evident by the inversion-recovery segment. Both at small and at large ratios, the individual Cramér-Rao bounds are similar.

3.3 The Cramér-Rao Bound

This section analyzes the SNR-efficiency of different optimized spin ensemble trajectories. Examining $rCRB(T_{1})$, we can observe that the purely $T_{1}$-optimized trajectories shine most at short experiment durations $T_{\text{exp}}$ (Fig. 3a). As the duration increases, the advantage of the $T_{1}$-optimized trajectories over the joint optimizations starts to vanish. At $T_{\text{exp}}=3.8$s, which is used for the in vivo experiments, we expect only a minor advantage of the $T_{1}$-optimized sequence in comparison to the joint optimization. In contrast, the $T_{2}$-optimized sequence shows a poor $T_{1}$-encoding power. Furthermore, we can see that the limit $\vartheta\leq\pi/4$ has only minor effects on the $T_{1}$ and joint optimizations (transparent vs. solid marks in Fig. 3). For comparison, we also show the Cramér-Rao bound of an optimized Look-Locker sequence 21 in which we also allowed the flip angle to vary over time. One can observe a substantially worse performance compared to the hybrid state, even if one is interested only in $T_{1}$.

Examining $rCRB(T_{2})$, we also find the difference of the purely $T_{2}$-optimized and the jointly optimized sequences to be rather small (Fig. 3b), especially when limiting the polar angle to $\pi/4$. The $T_{1}$-optimized patterns, on the other hand, show a poor $T_{2}$-encoding power. For comparison, we also depict the Cramér-Rao bound of a multi-spin-echo sequence. For short experiments, this sequence has the same $rCRB(T_{2})$ as the hybrid-state sequence that was optimized for $T_{2}$ without constraining $\vartheta$. In fact, both optimizations result in the same sequences in this limit (not shown here). For long experiments, the performance of the multi-spin-echo sequence does, however, degrade and the more general hybrid state allows for a more efficient sequence design.

So far, the Cramér-Rao bound was analyzed only at the specific $T_{1}=781~{}\text{ms}$ and $T_{2}=65~{}\text{ms}$ that were used for the optimization. Fig. 4 analyzes the performance of sequences optimized for these particular relaxation times over a larger parameter space, i.e. over the range of $T_{1}$ and $T_{2}$ values commonly found in biological tissue. We found that optimizing for a single set of relaxation times does not merely result in good performance in the proximity of this set of parameters, which is visualized by the red dot. Instead, we found good performance over a large area in the $T_{1}$-$T_{2}$-space with the lowest Cramér-Rao bound not being located at the relaxation times used for the optimizations.

In the supplementary material, we bridge the gap between the Cramér-Rao bound and the correlation coefficients, which are used to reconstruct the parameter maps in MRF 5, 24 (Fig. LABEL:fig:correlation_WM and LABEL:fig:correlation_max).

3.4 In Vivo Experiments

Fig. 5 depicts quantitative in vivo maps acquired with the patterns shown in Fig. 2k,n,q and with one radial k-space spoke acquired for each time frame. The total scan time was approximately $3.8$s for each of the three experiments. Systematic deviations can be noted in the measured relaxation times (Fig. 5). These variations are most likely caused by magnetization transfer 30, which is not addressed in the present work. Focusing on the noise properties, one can observe a good agreement between the Cramér-Rao bound predictions and the resulting noise level in the parameter maps. As anticipated, the purely $T_{1}$-optimized pattern and the jointly-optimized pattern achieve a very similar SNR level in the $T_{1}$-maps (cf. Fig. 3a), while the $T_{2}$-optimized pattern is slightly worse. The gap is substantially larger in the $T_{2}$-maps, where the $T_{2}$-optimized and the jointly-optimized patterns achieve similar performance, but the $T_{1}$-optimized pattern has a substantially larger standard deviation (cf. the ROI analysis displayed in the corner of Fig. 5).

4 Discussion

The purpose of this paper is to shed some light on the optimal encoding of spin relaxation times. In hybrid state, the spin dynamics is trapped in a single dimension and can be described by a single, uncoupled differential equation (Eq. (1)), which may be conveniently visualized (Fig. 1). Insights gleaned from intuitive understanding of the governing equation were found to align well with numerical optimizations (Fig. 2). Numerically searching for the optimal spin trajectories is a non-convex problem so we can only speculate about the optimality of our results. However, the notable correspondence to the intuitive understanding of the governing differential equation, together with the number of reproducible features in the optimized trajectories, gives us confidence that the obtained numerical results are closely related to the truly optimal trajectories.

The Cramér-Rao bound was used as figure of merit since it provides a universal lower bound for the noise transfer from a time series to the estimated parameters. It also provides the flexibility to tailor the acquisition to, for example, encode $T_{1}$ while accounting for correlations in $T_{2}$-space. Within the limits of this analysis—most notably our neglect of the incomplete spatial encoding and the non-convexity of the optimization problems—we showed that the hybrid state has superior SNR properties compared to traditional methods, namely Look-Locker for $T_{1}$-encoding and multi-spin-echo sequences for $T_{2}$-encoding. Moreover, we found that the joint optimizations have nearly the same performance in encoding each relaxation time compared to hybrid-state sequences that are specialized to encode a single parameter. Thus, the second relaxation time can be measured at virtually no extra cost.

Here, we optimized the sequences only for a single set of relaxation times, which correspond to white matter. However, we found that such optimized sequences perform well over a large range of relaxation times. If one wanted to extend this area further, one could optimize for the sum of the Cramér-Rao bounds at different relaxation times.

Our least-constrained optimizations used the limits $-1\leq r(t)\leq 1$ and $0\leq{\vartheta}(t)\leq\pi/2$. This unusual definition of spherical coordinates was chosen to deliberately prohibit magnetization from revisiting the southern hemisphere after passing through the origin. This allows us to implement the RF-pulses without major violations of the small tip-angle approximation. However, we demonstrated the benefits of the southern hemisphere and for long experiments as might be required, e.g., for 3D imaging, it might be desirable to revisit the southern hemisphere. This could be achieved by simply changing the limits to $0\leq r(t)\leq 1$ and $0\leq{\vartheta}(t)\leq\pi$. However, this would change the optimization landscape, and the convergence behavior might be impaired. Alternatively, one could use anti-periodic boundary conditions, as demonstrated in Ref. 11.

The present work focuses on signal-to-noise properties and the employed model is derived from the standard Bloch equation. The in vivo data shown in Fig. 5 are corrupted by systematic errors such as $B_{0}$- and $B_{1}$-inhomogeneities, as well as magnetization transfer effects. The first two problems can be corrected during the reconstruction process 5, 31, 7, 25, and it may be possible to also address magnetization transfer effects using a multi pool model 30. Future work will include numerical optimizations that incorporate these effects in order to minimize correlations between the corresponding parameters and the spin relaxation times.