Why map T₁?

• T₁, the spin–lattice relaxation time: the rate at which excited spins return to equilibrium with their surroundings
• A quantitative biomarker
• Needed to optimize protocols, e.g. the Ernst angle for maximum SNR in spoiled gradient echo (Ernst 1966)
• A correction factor for other qMRI techniques: quantitative MT (Cercignani 2005; Yarnykh 2002) and DCE (Li 2018; Sung 2013)

Three approaches in this chapter

The inversion recovery experiment

• A 180° RF pulse inverts the magnetization at thermal equilibrium
• M_z recovers through spin–lattice relaxation for a set delay, the inversion time (TI)
• A 90° excitation, then a readout (spin echo or gradient echo), takes a snapshot of M_z at that TI
• Repeat for several TIs, then fit the recovery curve to estimate T₁

Where it comes from

• First developed for NMR in the 1940s (Hahn 1949; Drain 1949)
• The first T₁ map used saturation recovery: a 90° preparation pulse instead of 180° (Pykett 1978)

Strengths, and the catch

So IR is scarcely used in practice, but it is the reference measurement for developing and comparing T₁ methods. Several variations make it practical for some applications (Messroghli 2004; Piechnik 2010).

Signal modelling: the general equation

Steady-state M_z from the Bloch equations for {θ₁₈₀ – TI – θ₉₀ – (TR − TI)}:

• M_z is the longitudinal magnetization just before the θ₉₀ pulse
• For the in-phase real signal, multiply by \(k\sin(\theta_{90})\,e^{-TE/T_2}\), where k is a constant

Simplifying: constants and an ideal inversion

Group the constants for each measurement (the same TE and θ₉₀ at every TI) and assume an ideal inversion pulse:

C groups the first three terms and the denominator of Eq. 2.1.

Simplifying further: the long-TR approximation

If TR is long enough for full relaxation (TR > 5 T₁), drop the last term of Eq. 2.2:

• Simple in both equation and experiment: the most widely used IR signal model
• Many implementations acquire magnitude-only images, so the measured signal is proportional to |Eq. 2.3|

When the approximation breaks down

In practice, Eq. 2.1 is the better choice for simulating an IR experiment:

• TR is chosen to exceed 5 T₁ of the tissue of interest, which is rarely the longest T₁ present: long enough for white matter, but not for CSF
• Eq. 2.3 also assumes ideal inversion pulses, rarely the case because of slice profile effects

Next: magnitude signal for TR = 5 s and T₁ from 250 ms to 5 s, computed with both equations.

Data fitting: what to watch for

• Magnitude-only images: polarity inversion restores the exponential form, but it is sensitive to noise: Rician noise gives a non-zero level at the signal null
• Phase data available: a phase term must be added to the fitting equation (Barral 2010)
• Eq. 2.3 must only be used for long-TR data (TR > 5 T₁), which in practice is rarely true for every tissue

Fitting algorithms

Early methods, built for the computers of the time

• Null method: each T₁ has a unique zero crossing (Pykett 1983)
• Linear fit of a rearranged Eq. 2.3 on a semi-log plot (Fukushima 1981)

Today

• Non-linear least squares, e.g. Levenberg–Marquardt, on the approximate (Eq. 2.3) or general (Eq. 2.1) model

Reduced-dimension NLS (RD-NLS)

• Fits Eq. 2.1 up to 75× faster than Levenberg–Marquardt, without a precision penalty (Barral 2010)
• A simplified 5-parameter model is sufficient for accurate T₁ mapping:

• a and b are complex; for magnitude-only data, a 3-parameter fit of |Eq. 2.4| can be sufficient
• General: no assumption about TR, so as robust as Eq. 2.1 when every parameter except TI is kept constant
• Open-source code released with the paper, and available as a qMRLab T₁ model

Why it is the gold standard

• A typical protocol uses a long TR and enough inversion times for a stable fit (typically 5 or more), covering [0, TR]
• A wide dynamic range, up to [−kM₀, kM₀]: many TIs with high SNR along the recovery curve (Fukushima 1981)
• Largely insensitive to flip-angle inaccuracies and imperfect spoiling (Stikov 2015): only TI changes, and a single acquisition is made per TR

Protocol design

• Avoid inversion times that null the signal of the tissue of interest: those magnitude images are dominated by Rician noise, which can hurt the fit at low SNR
• Runs on standard pulse sequences available on most scanners, using a customized protocol
• Needs no extra calibration measurements

Example: the interactive preprint of the ISMRM Reproducible Research Study Group 2020 challenge on IR T₁ mapping (Boudreau 2023).

The main drawback: time

• Long TRs are needed: a few T₁ for general models (Eqs. 2.1, 2.4), up to 5 T₁ for the long-TR model (Eq. 2.3)
• One TI per TR (2–5 s) and one phase-encode line per excitation (~100–200 lines): 10–25 minutes for a single-slice T₁ map
• So whole-organ maps are hard to get in clinically feasible times

Still useful as a reference measurement, or for a single-slice T₁ estimate used to optimize other pulse sequences.

Faster variants: Look-Locker

• After the inversion, a periodic train of small pulses θ samples the recovery (Look 1970):

• Samples the inversion recovery curve much more efficiently, at the cost of lower SNR
• Each TI depends on the previous θ pulses, so it is more sensitive to B₁ inhomogeneity and imperfect spoiling (Gai 2013; Stikov 2015)
• Widely used for rapid T₁ mapping; MOLLI and ShMOLLI for cardiac T₁ mapping (Messroghli 2004; Piechnik 2010)

Saturation recovery

• The inversion pulse is replaced with a saturation pulse: θ₉₀ – TI – θ₉₀; used for the very first T₁ map (Pykett 1978)
• No long TR needed: every θ₉₀ resets M_z to the same initial state
• But:
  • TIs still need to reach the order of T₁
  • The dynamic range is halved: [0, M₀]
  • The shortest TIs, the fastest to acquire, have the lowest SNR

The idea

• Acquire several spoiled gradient echo images, each with a different excitation flip angle θ_n (θ_i ≠ θ_j)
• Very short TRs, on the order of 10 ms
• The steady-state signal is very sensitive to T₁ over a wide range of flip angles

Origins and strengths

• Originates from NMR: adopted for its time efficiency and its ability to measure accurate T₁ values simultaneously over a wide range (Christensen 1974; Gupta 1977)
• For imaging, VFA can use a 3D acquisition instead of multislice:
  • higher SNR
  • reduced slice profile effects

The catch: flip angle accuracy

• The signal is very sensitive to inaccuracies in the flip angle
• The nominal flip angle (set at the scanner) differs from the actual one experienced by the spins: up to ±30 % at 3.0 T, increasing with field strength
• VFA typically needs another quantitative map, the transmit RF amplitude B₁⁺ (B₁), to calibrate the flip angle (Sled 1998)
• The extra B₁ map reduces the time savings, and its errors propagate into the T₁ map (Boudreau 2017; Lee 2017)

Signal modelling

Steady-state M_z of an ideal VFA experiment, solved analytically from the Bloch equations for the spoiled gradient echo sequence {θ_n – TR}:

• M₀: magnetization at thermal equilibrium; TR: repetition time; θ_n: excitation flip angle

Next: M_z curves for three brain-tissue T₁ values over a range of θ_n and TR.

The Ernst angle

The flip angle that maximizes the steady-state signal depends on T₁ and TR. It can be solved analytically from Eq. 2.5 (Ernst 1966):

Model assumptions

The closed-form solution (Eq. 2.5) assumes:

• The longitudinal magnetization has reached a steady state after a large number of TRs
• The transverse magnetization is perfectly spoiled at the end of each TR

These may not hold if care is not taken. Bloch simulations, which solve the Bloch equations numerically for a set of spins at each time point, give a more realistic estimate when a steady state is not reached.

Spoiling

• Sufficient spoiling is likely the most challenging parameter to control in a VFA experiment
• A combination of gradient spoiling and RF phase spoiling is typically recommended (Handbook 2004; Zur 1991)
• RF phase increment after each TR: \(\Phi_{n}=\Phi_{n-1}+n\Phi_{0}\), with \(\Phi_0 = 117°\) (Zur 1991)
• Very strong gradients introduce diffusion effects that further improve spoiling (Yarnykh 2010)

Data fitting: linearizing the signal

• Non-linear least squares on Eq. 2.5 (two free variables: T₁ and M₀) is valid, but rarely done in practice
• A simple refactoring of Eq. 2.5, without any approximation, gives the form y = mx + b (Gupta 1977):

The last term is the same for every θ_n, so it can be grouped into the constant:

T₁ from the slope

A linear regression of S_n/sin(θ_n) against S_n/tan(θ_n) gives T₁ from its slope:

• With only two flip angles, a very common protocol, the slope is simply the elementary slope between the two points

Next: both forms of the signal (Eq. 2.5 and Eq. 2.8) for a noisy dataset.

Protocol design: which flip angles?

• Most VFA experiments use the minimum number of flip angles, two, to minimize acquisition time
• For a single T₁, the best precision comes from the two flip angles that give 71 % of the maximum signal (at the Ernst angle) (Deoni 2003; Schabel 2008)
• Time allowing, extra flip angles are acquired at higher values and in between, where differences between tissue T₁ values are larger
• With more than two angles, the points no longer share the same noise weighting, which may bias linear fits at low SNR: use non-linear least squares or a weighted linear fit (Chang 2008)

Correcting the flip angle with a B₁ map

• The excitation RF field (B₁⁺) of a loaded coil varies in space, unless RF shimming is available, which is not common at clinical field strengths (Sled 1998)
• Correct the flip angle voxelwise with a B₁ map acquired in the same session:

• B₁ is normalized: unitless, and equal to 1 where the RF field has the expected amplitude

Next: fitted T₁ from Monte Carlo data simulated with biased flip angles, fitted without and with B₁ correction.

Pitfalls

• Accuracy is very sensitive to the pulse sequence implementation, so VFA is less robust than inversion recovery (Baudrexel 2017; Lutti 2013; Stikov 2015)
• Insufficient spoiling can bias T₁ by up to 30 % relative to inversion recovery (Stikov 2015)
• Accuracy and precision depend strongly on the quality of the B₁ map (Lee 2017), which varies between implementations (Boudreau 2017)
• Omitting B₁ correction is not recommended: most scanners have default sequences that, with careful settings, give B₁ maps of sufficient quality very rapidly (Boudreau 2017; Samson 2006; Wang 2005)

Why VFA is still widely used

• One of the most widely used T₁ mapping methods in research
• Rapid acquisition and rapid image processing
• Widespread availability
• A great candidate for other quantitative protocols: quantitative MT (Cercignani 2005; Yarnykh 2002) and DCE (Li 2018; Sung 2013)

Dictionary-based techniques

• Growing in popularity: available on clinical scanners, rapid scan times, fast post-processing
• They use numerical dictionaries: databases of pre-calculated signal values simulated for a wide range of tissue and protocol combinations
• Examples for T₁: MR Fingerprinting (Ma 2013), some compressed sensing methods (Doneva 2010; Li 2012), and MP2RAGE (Marques 2010)
• Two categories:
  • dictionaries that exploit redundancy in parametric data to accelerate imaging (CS, MRF)
  • dictionaries that estimate quantitative maps from the reconstructed images (MP2RAGE)

MP2RAGE: a hybrid of IR and VFA

• An extension of the conventional MPRAGE sequence widely used in clinical studies (Haase 1989; Mugler 1990)
• A 180° inversion pulse prepares the magnetization for T₁ sensitivity at the start of each TR_MP2RAGE
• Then two images at two inversion times, using gradient echo (GRE) blocks with low flip angles and short TRs

Encoding, and why two images

• Within a GRE block, each excitation is followed by a constant in-plane ("y") phase encode, varied for each TR_MP2RAGE, and a different 3D ("z") phase encode, varied at each TR
• The centre of k-space in the 3D direction is acquired at each block's TI
• Motivation: an MPRAGE-like image with self-bias correction of B₀ and B₁⁻, and a first-order correction of B₁⁺
• Because two TIs are acquired (MPRAGE acquires one), T₁ can be inferred, giving quantitative T₁ maps

The MP2RAGE signal

For complex data (magnitude and phase, or real and imaginary), S_MP2RAGE combines the images at the two TIs (Marques 2010):

• Bounded between [−0.5, 0.5]
• The phase data help reduce some B₀ inhomogeneity effects

Real or polarity-restored data

For real data, or magnitude data with polarity restoration:

Signal equations: assumptions

• MP2RAGE is a hybrid of IR and VFA, so its signal equations are more complex
• A steady state is not reached during the short GRE train: the signal at the centre of k-space, which defines the contrast, depends on the number of phase-encoding steps
• The equations that follow assume the 3D phase-encoding dimension is fully sampled: no partial Fourier or parallel imaging

Derivation: appendix of Marques et al. 2010.

Signal at the two inversion times

Definitions

Steady state before the inversion

What S_MP2RAGE removes, and what remains

Data fitting with a dictionary

• No conventional minimization (e.g. Levenberg–Marquardt): MP2RAGE uses pre-calculated signal values over a wide range of T₁, then interpolates to find the T₁ that matches the observed signal
• Dictionaries can be generated before scanning, and interpolating is much faster than most fitting algorithms, so maps can be shown directly on the scanner console
• Signal curves must be unique for each parameter value: S_MP2RAGE is unique for each T₁ for a given protocol, so a [T₁, S_MP2RAGE] dictionary is enough (Marques 2010)
• B₁ sensitivity is reduced with protocol optimization, and further by adding B₁ as a dictionary dimension [T₁, B₁, S_MP2RAGE], at the cost of an extra B₁ map (Marques 2013)

Software

• MP2RAGE is increasingly distributed by MRI vendors, typically with a fitting package that reconstructs T₁ maps online
• Several open-source packages are also available (Marques 2017; de Hollander 2017)
• For new users these are recommended over programming one from scratch: there are many potential pitfalls, e.g. adjusting the equations for partial Fourier or parallel imaging

Benefits

• Widespread availability and turnkey acquisition and fitting: a main reason for growing interest in quantitative T₁ in clinical and neuroscience studies
• Highly reproducible brain T₁ in an inter- and intra-site study at eight sites (same hardware and software, 7T) (Voelker 2016)
• Among the fastest acquisition and post-processing times of quantitative T₁ methods
• Very high resolution: 1 mm isotropic at 3T and submillimeter at 7T, both in under 10 min (Fujimoto 2014), opening the door to cortical studies (Waehnert 2014; Beck 2018; Haast 2018)

Pitfalls

• Reproducibility depends on using one pre-calculated dictionary: different dictionaries (e.g. across sites and vendors) give minor differences in T₁ for the same data
• Residual B₁ sensitivity can be large enough to need correction with a measured B₁ map (Marques 2013; Haast 2018), but B₁ mapping brings its own errors: choose the method and post-processing (e.g. filtering) carefully (Boudreau 2017)
• The equations, and so the dictionaries, assume monoexponential longitudinal relaxation: suboptimal estimates of the long T₁ component under a biexponential model (Rioux 2016), more important at higher fields

Summary

Exercises and discussion

• The chapter ends with nine exercises that use the interactive figures, the notebooks (on Binder) and qMRLab
• Questions and answers: github.com/qMRLab/mooc/discussions

Open the exercises in the book →