Display content:

In [1]:

```
% **Blog post code introduction**
%
% Congrats on activating the "All cells" option in this interactive blog post =D
%
% Below, several new HTML blocks have appears prior to the figures, displaying the Octave/MATLAB code that was used to generate the figures in this blog post.
%
% If you want to reproduce the data on your own local computer, you simply need to have qMRLab installed in your Octave/MATLAB path and run the "startup.m" file, as is shown below.
%
% If you want to get under the hood and modify the code right now, you can do so in the Jupyter Notebook of this blog post hosted on MyBinder. The link to it is in the introduction above.
```

In [2]:

Widely considered the gold standard for T_{1} mapping, the inversion recovery technique estimates T_{1} values by fitting the signal recovery curve acquired at different delays after an inversion pulse (180°). In a typical inversion recovery experiment (Figure 1), the magnetization at thermal equilibrium is inverted using a 180° RF pulse. After the longitudinal magnetization recovers through spin-lattice relaxation for predetermined delay (“inversion time”, TI), a 90° excitation pulse is applied, followed by a readout imaging sequence (typically a spin-echo or gradient-echo readout) to create a snapshot of the longitudinal magnetization state at that TI.

Inversion recovery was first developed for NMR in the 1940s (Hahn 1949; Drain 1949), and the first T_{1} map was acquired using a saturation-recovery technique (90° as a preparation pulse instead of 180°) by (Pykett and Mansfield 1978). Some distinct advantages of inversion recovery are its large dynamic range of signal change and an insensitivity to pulse sequence parameter imperfections (Stikov et al. 2015). Despite its proven robustness at measuring T_{1}, inversion recovery is scarcely used in practice, because conventional implementations require repetition times (TRs) on the order of 2 to 5 T_{1} (Steen et al. 1994), making it challenging to acquire whole-organ T_{1} maps in a clinically feasible time. Nonetheless, it is continuously used as a reference measurement during the development of new techniques, or when comparing different T_{1} mapping techniques, and several variations of the inversion recovery technique have been developed, making it practical for some applications (Messroghli et al. 2004; Piechnik et al. 2010).

The steady-state longitudinal magnetization of an inversion recovery experiment can be derived from the Bloch equations for the pulse sequence {θ_{180} – TI – θ_{90} – (TR-TI)}, and is given by:

where M_{z} is the longitudinal magnetization prior to the θ_{90} pulse. If the in-phase real signal is desired, it can be calculated by multiplying Eq. 1 by *k*sin(θ_{90})e^{-TE/T2}, where *k* is a constant. This general equation can be simplified by grouping together the constants for each measurements regardless of their values (i.e. at each TI, same TE and θ_{90} are used) and assuming an ideal inversion pulse:

where the first three terms and the denominator of Eq. 1 have been grouped together into the constant C. If the experiment is designed such that TR is long enough to allow for full relaxation of the magnetization (TR > 5T_{1}), we can do an additional approximation by dropping the last term in Eq. 2:

The simplicity of the signal model described by Eq. 3, both in its equation and experimental implementation, has made it the most widely used equation to describe the signal evolution in an inversion recovery T_{1} mapping experiment. The magnetization curves are plotted in Figure 2 for approximate T_{1} values of three different tissues in the brain. Note that in many practical implementations, magnitude-only images are acquired, so the signal measured would be proportional to the absolute value of Eq. 3.

In [3]:

```
%% MATLAB/OCTAVE CODE
% Adds qMRLab to the path of the environment
cd ../qMRLab
startup
```

In [4]:

```
%% MATLAB/OCTAVE CODE
% Code used to generate the data required for Figure 2 of the blog post
clear all
%% Setup parameters
% All times are in seconds
% All flip angles are in degrees
params.TR = 5.0;
params.TI = linspace(0.001, params.TR, 1000);
params.TE = 0.004;
params.T2 = 0.040;
params.EXC_FA = 90; % Excitation flip angle
params.INV_FA = 180; % Inversion flip angle
params.signalConstant = 1;
%% Calculate signals
%
% The option 'GRE-IR' selects the analytical equations for the gradient echo readout inversion recovery experiment
% The option '4' is a flag that selects the long TR approximation of the analytical solution (TR>>T1), Eq. 3 of the blog post.
%
% To see all the options available, run `help inversion_recovery.analytical_solution`
% White matter
params.T1 = 0.900; % in seconds
signal_WM = inversion_recovery.analytical_solution(params, 'GRE-IR', 4);
% Grey matter
params.T1 = 1.500; % in seconds
signal_GM = inversion_recovery.analytical_solution(params, 'GRE-IR', 4);
% CSF
params.T1 = 4.000; % in seconds
signal_CSF = inversion_recovery.analytical_solution(params, 'GRE-IR', 4);
```

In [5]:

In [6]:

_{1} of the tissue-of-interest, which rarely coincides with the longest T_{1} present (e.g. TR may be sufficiently long for white matter, but not for CSF which could also be present in the volume). Equation 3 also assumes ideal inversion pulses, which is rarely the case due to slice profile effects. Figure 3 displays the inversion recovery signal magnitude (complete relaxation normalized to 1) of an experiment with TR = 5 s and T_{1} values ranging between 250 ms to 5 s, calculated using both equations.

In [7]:

```
%% MATLAB/OCTAVE CODE
% Code used to generate the data required for Figure 3 of the blog post
clear all
%% Setup parameters
% All times are in seconds
% All flip angles are in degrees
params.TR = 5.0;
params.TI = linspace(0.001, params.TR, 1000);
params.TE = 0.004;
params.T2 = 0.040;
params.EXC_FA = 90; % Excitation flip angle
params.INV_FA = 180; % Inversion flip angle
params.signalConstant = 1;
T1_range = 0.25:0.25:5; % in seconds
%% Calculate signals
%
% The option 'GRE-IR' selects the analytical equations for the gradient echo readout inversion recovery experiment
% The option '1' is a flag that selects full analytical solution equation (no approximation), Eq. 1 of the blog post.
% The option '4' is a flag that selects the long TR approximation of the analytical solution (TR>>T1), Eq. 3 of the blog post.
%
% To see all the options available, run `help inversion_recovery.analytical_solution`
for ii = 1:length(T1_range)
params.T1 = T1_range(ii);
signal_T1_Eq1{ii} = inversion_recovery.analytical_solution(params, 'GRE-IR', 1);
signal_T1_Eq3{ii} = inversion_recovery.analytical_solution(params, 'GRE-IR', 4);
end
```

In [8]:

In [9]:

Several factors impact the choice of the inversion recovery fitting algorithm. If only magnitude images are available, then a polarity-inversion is often implemented to restore the non-exponential magnitude curves (Figure 3) into the exponential form (Figure 2). This process is sensitive to noise due to the Rician noise creating a non-zero level at the signal null. If phase data is also available, then a phase term must be added to the fitting equation (Barral et al. 2010). Equation 3 must only be used to fit data for the long TR regime (TR > 5T_{1}), which in practice is rarely satisfied for all tissues in subjects.

Early implementations of inversion recovery fitting algorithms were designed around the computational power available at the time. These included the “null method” (Pykett et al. 1983), assuming that each T_{1} value has unique zero-crossings (see Figure 2), and linear fitting of a rearranged version of Eq. 3 on a semi-log plot (Fukushima & Roeder 1981). Nowadays, a non-linear least-squares fitting algorithm (e.g. Levenberg-Marquardt) is more appropriate, and can be applied to either approximate or general forms of the signal model (Eq. 3 or Eq. 1). More recent work (Barral et al. 2010) demonstrated that T_{1} maps can also be fitted much faster (up to 75 times compared to Levenberg-Marquardt) to fit Eq. 1 – without a precision penalty – by using a reduced-dimension non-linear least squares (RD-NLS) algorithm. It was demonstrated that the following simplified 5-parameter equation can be sufficient for accurate T_{1} mapping:

where *a* and *b* are complex values. If magnitude-only data is available, a 3-parameter model can be sufficient by taking the absolute value of Eq. 4. While the RD-NLS algorithms are too complex to be presented here (the reader is referred to the paper, (Barral et al. 2010)), the code for these algorithms was released open-source along with the original publication, and is also available as a qMRLab T_{1} mapping model. One important thing to note about Eq. 4 is that it is general – no assumption is made about TR – and is thus as robust as Eq. 1 as long as all pulse sequence parameters other than TI are kept constant between each measurement. Figure 4 compares simulated data (Eq. 1) using a range of TRs (1.5T_{1} to 5T_{1}) fitted using either RD-NLS & Eq. 4 or a Levenberg-Marquardt fit of Eq. 2.

In [10]:

```
%% MATLAB/OCTAVE CODE
% Code used to generate the data required for Figure 4 of the blog post
clear all
%% Setup parameters
% All times are in milliseconds
% All flip angles are in degrees
params.TI = 50:50:1500;
TR_range = 1500:50:5000;
params.EXC_FA = 90;
params.INV_FA = 180;
params.T1 = 1000;
%% Calculate signals
%
% The option 'GRE-IR' selects the analytical equations for the gradient echo readout inversion recovery experiment
% The option '1' is a flag that selects full analytical solution equation (no approximation), Eq. 1 of the blog post.
%
% To see all the options available, run `help inversion_recovery.analytical_solution`
for ii = 1:length(TR_range)
params.TR = TR_range(ii);
Mz_analytical(ii,:) = inversion_recovery.analytical_solution(params, 'GRE-IR', 1);
end
%% Fit data using Levenberg-Marquardt with the long TR approximation equation
%
% The option '4' is a flag that selects the long TR approximation of the analytical solution (TR>>T1), Eq. 3 of the blog post.
%
% To see all the options available, run `help inversion_recovery.fit_lm`
for ii=1:length(TR_range)
fitOutput_lm{ii} = inversion_recovery.fit_lm(Mz_analytical(ii,:), params, 4);
T1_lm(ii) = fitOutput_lm{ii}.T1;
end
%% Fit data using the RDLS method (Barral), Eq. 4 of the blog post.
%
% Create a qMRLab inversion recovery model object and load protocol values
irObj = inversion_recovery();
irObj.Prot.IRData.Mat = params.TI';
for ii=1:length(TR_range)
data.IRData = Mz_analytical(ii,:);
fitOutput_barral{ii} = irObj.fit(data);
T1_barral(ii) = fitOutput_barral{ii}.T1;
end
```